v0.16.0
Loading...
Searching...
No Matches
Hdiv.cpp
Go to the documentation of this file.
1/** \file Hdiv.cpp
2
3 \brief Implementation of H-curl base
4
5 Based on Hierarchic Finite Element Bases on Unstructured Tetrahedral
6 Meshes, by Mark Ainsworth and Joe Coyle
7 Shape functions for MBTRI/MBTET and HCurl space
8
9*/
10
11using namespace MoFEM;
12
13boost::function<int(int)> AinsworthOrderHooks::broken_nbfacetri_edge_hdiv =
14 [](int p) { return p; };
15boost::function<int(int)> AinsworthOrderHooks::broken_nbfacetri_face_hdiv =
16 [](int p) { return p; };
17boost::function<int(int)> AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv =
18 [](int p) { return p; };
19boost::function<int(int)> AinsworthOrderHooks::broken_nbvolumetet_face_hdiv =
20 [](int p) { return p; };
21boost::function<int(int)> AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv =
22 [](int p) { return p; };
23
24MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(
25 int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3],
26 double *diff_phi_f_e[4][3], int gdim,
27 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
28 double *L, double *diffL,
29 const int dim)) {
30
31 MoFEMFunctionBeginHot;
32#ifndef NDEBUG
33 if (!diff_phi_f_e)
34 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
35 "expected to return derivatives");
36#endif
37
38 for (int ff = 0; ff < 4; ff++) {
39 CHKERR Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
40 &faces_nodes[3 * ff], p[ff], N, diffN, phi_f_e[ff], diff_phi_f_e[ff],
41 gdim, 4, base_polynomials);
42 }
43
44 MoFEMFunctionReturnHot(0);
45}
46
47static inline double get_ksi(const double n0, const double n1) {
48 static_assert(
49 AinsworthHdivOptions::edge_coordinate ==
50 AinsworthHdivOptions::EdgeCoordinate::PAPER ||
51 AinsworthHdivOptions::edge_coordinate ==
52 AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_O ||
53 AinsworthHdivOptions::edge_coordinate ==
54 AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_I,
55 "Unsupported Ainsworth H(div) edge coordinate");
56
57 if constexpr (AinsworthHdivOptions::edge_coordinate ==
58 AinsworthHdivOptions::EdgeCoordinate::PAPER) {
59 return n1 - n0;
60 } else if constexpr (AinsworthHdivOptions::edge_coordinate ==
61 AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_O) {
62 return 1 - 2 * n0;
63 } else if constexpr (AinsworthHdivOptions::edge_coordinate ==
64 AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_I) {
65 return 2 * n1 - 1;
66 }
67 return 0;
68}
69
70static inline FTensor::Tensor1<double, 3> get_diff_ksi(
71 const double *diffN, const int n0, const int n1) {
72 static_assert(
73 AinsworthHdivOptions::edge_coordinate ==
74 AinsworthHdivOptions::EdgeCoordinate::PAPER ||
75 AinsworthHdivOptions::edge_coordinate ==
76 AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_O ||
77 AinsworthHdivOptions::edge_coordinate ==
78 AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_I,
79 "Unsupported Ainsworth H(div) edge coordinate");
80
81 FTensor::Tensor1<double, 3> t_diff_ksi;
82 for (int dd = 0; dd != 3; ++dd) {
83 if constexpr (AinsworthHdivOptions::edge_coordinate ==
84 AinsworthHdivOptions::EdgeCoordinate::PAPER) {
85 t_diff_ksi(dd) = diffN[3 * n1 + dd] - diffN[3 * n0 + dd];
86 } else if constexpr (AinsworthHdivOptions::edge_coordinate ==
87 AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_O) {
88 t_diff_ksi(dd) = -2 * diffN[3 * n0 + dd];
89 } else if constexpr (AinsworthHdivOptions::edge_coordinate ==
90 AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_I) {
91 t_diff_ksi(dd) = 2 * diffN[3 * n1 + dd];
92 }
93 }
94 return t_diff_ksi;
95}
96
97MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
98 int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3],
99 double *diff_phi_f_e[3], int gdim, int nb,
100 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
101 double *L, double *diffL,
102 const int dim)) {
103
104 constexpr int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
105 constexpr int face_opposite_edges_node[] = {2, 0, 1};
106 FTENSOR_INDEX(3, i);
107 FTENSOR_INDEX(3, j);
108 FTENSOR_INDEX(3, k);
109
110 MoFEMFunctionBeginHot;
111 if (p < 1)
112 MoFEMFunctionReturnHot(0);
113
114 FTensor::Tensor1<double, 3> t_edge_cross[3];
115 FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
116 FTensor::Tensor1<double, 3> t_diff_ksi[3];
117 if (diffN) {
118 t_node_diff_ksi[0] =
119 FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
120 t_node_diff_ksi[1] =
121 FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
122 t_node_diff_ksi[2] =
123 FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
124 t_node_diff_ksi[3] =
125 FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
126 for (int ee = 0; ee < 3; ee++) {
127 const int n0 = faces_nodes[face_edges_nodes[ee][0]];
128 const int n1 = faces_nodes[face_edges_nodes[ee][1]];
129 t_diff_ksi[ee] = get_diff_ksi(diffN, n0, n1);
130 t_edge_cross[ee](i) = levi_civita(i, j, k) * t_node_diff_ksi[n0](j) *
131 t_node_diff_ksi[n1](k);
132 }
133 } else {
134 for (int ee = 0; ee < 3; ee++) {
135 t_edge_cross[ee](0) = 1;
136 t_edge_cross[ee](1) = 0;
137 t_edge_cross[ee](2) = 0;
138 }
139 }
140 double psi_l[p + 1];
141
142 for (int ee = 0; ee != 3; ee++) {
143 const int i0 = faces_nodes[face_edges_nodes[ee][0]];
144 const int i1 = faces_nodes[face_edges_nodes[ee][1]];
145 const int iO = faces_nodes[face_opposite_edges_node[ee]];
146 auto t_psi_f_e = getFTensor1FromPtr<3>(phi_f_e[ee]);
147 for (int ii = 0; ii != nb * gdim; ii += nb) {
148 const double n0 = N[ii + i0];
149 const double n1 = N[ii + i1];
150 const double lambda = N[ii + iO];
151 double ksi = get_ksi(n0, n1);
152 base_polynomials(p, ksi, NULL, psi_l, NULL, 3);
153 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
154 int jj = 0;
155 for (int l = 0; l <= p - 1; l++) {
156 t_psi_f_e(i) = (lambda * t_psi_l) * t_edge_cross[ee](i);
157 ++t_psi_f_e;
158 ++t_psi_l;
159 ++jj;
160 }
161#ifndef NDEBUG
162 if (jj != NBFACETRI_AINSWORTH_EDGE_HDIV(p)) {
163 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
164 "wrong order %d != %d", jj, NBFACETRI_AINSWORTH_EDGE_HDIV(p));
165 }
166#endif
167 }
168 }
169
170 if (diff_phi_f_e) {
171 double diff_psi_l[3 * (p + 1)];
172 for (int ee = 0; ee != 3; ee++) {
173 const int i0 = faces_nodes[face_edges_nodes[ee][0]];
174 const int i1 = faces_nodes[face_edges_nodes[ee][1]];
175 const int iO = faces_nodes[face_opposite_edges_node[ee]];
176 auto t_diff_phi_f_e = getFTensor2HVecFromPtr<3, 3>(diff_phi_f_e[ee]);
177 for (int ii = 0; ii != nb * gdim; ii += nb) {
178 const double n0 = N[ii + i0];
179 const double n1 = N[ii + i1];
180 const double lambda = N[ii + iO];
181 double ksi = get_ksi(n0, n1);
182 base_polynomials(p, ksi, &t_diff_ksi[ee](0), psi_l, diff_psi_l, 3);
183 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
184 auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>{
185 &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]};
186 for (int l = 0; l <= p - 1; l++) {
187 t_diff_phi_f_e(i, j) =
188 (t_node_diff_ksi[iO](j) * t_psi_l + lambda * t_diff_psi_l(j)) *
189 t_edge_cross[ee](i);
190 ++t_diff_psi_l;
191 ++t_diff_phi_f_e;
192 ++t_psi_l;
193 }
194 }
195 }
196 }
197
198 MoFEMFunctionReturnHot(0);
199}
200
201MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions(
202 int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[],
203 double *diff_phi_f[], int gdim,
204 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
205 double *L, double *diffL,
206 const int dim)) {
207
208 MoFEMFunctionBeginHot;
209
210#ifndef NDEBUG
211 if (!diff_phi_f)
212 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
213 "expected to return derivatives");
214#endif
215
216 for (int ff = 0; ff < 4; ff++) {
217 CHKERR Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
218 &faces_nodes[3 * ff], p[ff], N, diffN, phi_f[ff], diff_phi_f[ff], gdim,
219 4, base_polynomials);
220 }
221 MoFEMFunctionReturnHot(0);
222}
223
224MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
225 int *face_nodes, int p, double *N, double *diffN, double *phi_f,
226 double *diff_phi_f, int gdim, int nb,
227 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
228 double *L, double *diffL,
229 const int dim)) {
230 FTENSOR_INDEX(3, i);
231 FTENSOR_INDEX(3, j);
232 FTENSOR_INDEX(3, k);
233
234 MoFEMFunctionBeginHot;
235 if (p < 3)
236 MoFEMFunctionReturnHot(0);
237
238 const int vert_i = face_nodes[1];
239 const int vert_j = face_nodes[2];
240 const int i0 = face_nodes[0];
241 FTensor::Tensor1<double, 3> t_cross;
242 FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
243 FTensor::Tensor1<double, 3> t_diff_ksi0i;
244 FTensor::Tensor1<double, 3> t_diff_ksi0j;
245
246 if (diffN) {
247 t_node_diff_ksi[0] =
248 FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
249 t_node_diff_ksi[1] =
250 FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
251 t_node_diff_ksi[2] =
252 FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
253 t_node_diff_ksi[3] =
254 FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
255 t_diff_ksi0i(i) = get_diff_ksi(diffN, i0, vert_i)(i);
256 t_diff_ksi0j(i) = get_diff_ksi(diffN, i0, vert_j)(i);
257 t_cross(i) = levi_civita(i, j, k) * t_node_diff_ksi[vert_i](j) *
258 t_node_diff_ksi[vert_j](k);
259 } else {
260 t_cross(0) = 1;
261 t_cross(1) = 0;
262 t_cross(2) = 0;
263 }
264
265 double psi_l[p + 1], diff_psi_l[3 * (p + 1)];
266 double psi_m[p + 1], diff_psi_m[3 * (p + 1)];
267 auto t_psi_f = getFTensor1FromPtr<3>(phi_f);
268
269 for (int ii = 0; ii < gdim; ii++) {
270
271 const int node_shift = ii * nb;
272 const double ni = N[node_shift + vert_i];
273 const double nj = N[node_shift + vert_j];
274 const double n0 = N[node_shift + i0];
275 const double ksi0i = get_ksi(n0, ni);
276 const double ksi0j = get_ksi(n0, nj);
277 double beta_0ij = n0 * ni * nj;
278 base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
279 base_polynomials(p, ksi0j, NULL, psi_m, NULL, 3);
280
281 int jj = 0;
282 int oo = 0;
283 for (; oo <= p - 3; oo++) {
284 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
285 for (int l = 0; l <= oo; l++) {
286 int m = oo - l;
287 if (m >= 0) {
288 t_psi_f(i) = (beta_0ij * t_psi_l * psi_m[m]) * t_cross(i);
289 ++t_psi_f;
290 }
291 ++t_psi_l;
292 ++jj;
293 }
294 }
295#ifndef NDEBUG
296 if (jj != NBFACETRI_AINSWORTH_FACE_HDIV(p)) {
297 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
298 "wrong order %d != %d", jj, NBFACETRI_AINSWORTH_FACE_HDIV(p));
299 }
300#endif
301 }
302
303 if (diff_phi_f) {
304 auto t_diff_phi_f = getFTensor2HVecFromPtr<3, 3>(diff_phi_f);
305
306 for (int ii = 0; ii < gdim; ii++) {
307
308 const int node_shift = ii * nb;
309 const double ni = N[node_shift + vert_i];
310 const double nj = N[node_shift + vert_j];
311 const double n0 = N[node_shift + i0];
312 const double ksi0i = get_ksi(n0, ni);
313 const double ksi0j = get_ksi(n0, nj);
314 double beta_0ij = n0 * ni * nj;
315 FTensor::Tensor1<double, 3> t_diff_beta_0ij;
316 t_diff_beta_0ij(i) = (ni * nj) * t_node_diff_ksi[i0](i) +
317 (n0 * nj) * t_node_diff_ksi[vert_i](i) +
318 (n0 * ni) * t_node_diff_ksi[vert_j](i);
319 base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
320 base_polynomials(p, ksi0j, &t_diff_ksi0j(0), psi_m, diff_psi_m, 3);
321
322 int jj = 0;
323 int oo = 0;
324 for (; oo <= p - 3; oo++) {
325 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
326 auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>{
327 &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]};
328 for (int l = 0; l <= oo; l++) {
329 int m = oo - l;
330 if (m >= 0) {
331 auto t_diff_psi_m =
332 FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>{
333 &diff_psi_m[m], &diff_psi_m[p + 1 + m],
334 &diff_psi_m[2 * p + 2 + m]};
335 t_diff_phi_f(i, j) = ((t_psi_l * psi_m[m]) * t_diff_beta_0ij(j) +
336 (beta_0ij * psi_m[m]) * t_diff_psi_l(j) +
337 (beta_0ij * t_psi_l) * t_diff_psi_m(j)) *
338 t_cross(i);
339 ++t_diff_phi_f;
340 }
341 ++t_psi_l;
342 ++t_diff_psi_l;
343 ++jj;
344 }
345 }
346#ifndef NDEBUG
347 if (jj != NBFACETRI_AINSWORTH_FACE_HDIV(p)) {
348 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
349 "wrong order %d != %d", jj, NBFACETRI_AINSWORTH_FACE_HDIV(p));
350 }
351#endif
352 }
353 }
354 MoFEMFunctionReturnHot(0);
355}
356
357MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(
358 int p, double *N, double *diffN, double *phi_v_e[6],
359 double *diff_phi_v_e[6], int gdim,
360 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
361 double *L, double *diffL,
362 const int dim)) {
363
364 constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
365 {0, 3}, {1, 3}, {2, 3}};
366
367 MoFEMFunctionBeginHot;
368 if (p < 2)
369 MoFEMFunctionReturnHot(0);
370 if (diffN == NULL) {
371 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
372 }
373
374 FTensor::Tensor1<double, 3> t_coords[4] = {
375 FTensor::Tensor1<double, 3>(0., 0., 0.),
376 FTensor::Tensor1<double, 3>(1., 0., 0.),
377 FTensor::Tensor1<double, 3>(0., 1., 0.),
378 FTensor::Tensor1<double, 3>(0., 0., 1.)};
379 FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
380 FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
381 FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
382 FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
383 FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
384
385 FTENSOR_INDEX(3, i);
386 FTENSOR_INDEX(3, j);
387
388 FTensor::Tensor1<double, 3> t_tou_e;
389 FTensor::Tensor1<double, 3> t_diff_ksi0i;
390 FTensor::Tensor1<double, 3> t_diff_beta_e;
391
392 double psi_l[p + 1];
393 double diff_psi_l[3 * (p + 1)];
394
395 for (int ee = 0; ee != 6; ee++) {
396 t_tou_e(i) =
397 t_coords[edges_nodes[ee][1]](i) - t_coords[edges_nodes[ee][0]](i);
398 t_diff_ksi0i = get_diff_ksi(diffN, edges_nodes[ee][0], edges_nodes[ee][1]);
399 auto t_phi_v_e = getFTensor1FromPtr<3>(phi_v_e[ee]);
400 for (int ii = 0; ii != gdim; ii++) {
401 const int node_shift = ii * 4;
402 const double ni = N[node_shift + edges_nodes[ee][1]];
403 const double n0 = N[node_shift + edges_nodes[ee][0]];
404 const double beta_e = ni * n0;
405 double ksi0i = get_ksi(n0, ni);
406 base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
407 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
408 for (int l = 0; l <= p - 2; l++) {
409 t_phi_v_e(i) = (beta_e * t_psi_l) * t_tou_e(i);
410 ++t_phi_v_e;
411 ++t_psi_l;
412 }
413 }
414 }
415
416 if (diff_phi_v_e) {
417 for (int ee = 0; ee != 6; ee++) {
418 t_tou_e(i) =
419 t_coords[edges_nodes[ee][1]](i) - t_coords[edges_nodes[ee][0]](i);
420 t_diff_ksi0i =
421 get_diff_ksi(diffN, edges_nodes[ee][0], edges_nodes[ee][1]);
422 auto t_diff_phi_v_e = getFTensor2HVecFromPtr<3, 3>(diff_phi_v_e[ee]);
423 for (int ii = 0; ii != gdim; ii++) {
424 const int node_shift = ii * 4;
425 const double ni = N[node_shift + edges_nodes[ee][1]];
426 const double n0 = N[node_shift + edges_nodes[ee][0]];
427 const double beta_e = ni * n0;
428 const double ksi0i = get_ksi(n0, ni);
429 t_diff_beta_e(i) = ni * t_node_diff_ksi[edges_nodes[ee][0]](i) +
430 t_node_diff_ksi[edges_nodes[ee][1]](i) * n0;
431 base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
432 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
433 auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
434 diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
435 for (int l = 0; l <= p - 2; l++) {
436 t_diff_phi_v_e(i, j) =
437 (t_diff_beta_e(j) * t_psi_l + beta_e * t_diff_psi_l(j)) *
438 t_tou_e(i);
439 ++t_diff_phi_v_e;
440 ++t_diff_psi_l;
441 ++t_psi_l;
442 }
443 }
444 }
445 }
446
447 MoFEMFunctionReturnHot(0);
448}
449
450MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(
451 int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[],
452 int gdim,
453 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
454 double *L, double *diffL,
455 const int dim)) {
456
457 constexpr int faces_nodes[4][3] = {
458 {0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
459
460 MoFEMFunctionBeginHot;
461 if (p < 3)
462 MoFEMFunctionReturnHot(0);
463
464 FTensor::Tensor1<double, 3> t_coords[4] = {
465 FTensor::Tensor1<double, 3>(0., 0., 0.),
466 FTensor::Tensor1<double, 3>(1., 0., 0.),
467 FTensor::Tensor1<double, 3>(0., 1., 0.),
468 FTensor::Tensor1<double, 3>(0., 0., 1.)};
469
470 FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
471 FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
472 FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
473 FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
474 FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
475
476 FTENSOR_INDEX(3, i);
477 FTENSOR_INDEX(3, j);
478
479 FTensor::Tensor1<double, 3> t_tau0i[4], t_tau0j[4];
480 FTensor::Tensor1<double, 3> t_diff_ksi0i[4], t_diff_ksi0j[4];
481 for (int ff = 0; ff != 4; ff++) {
482 const int v0 = faces_nodes[ff][0];
483 const int vi = faces_nodes[ff][1];
484 const int vj = faces_nodes[ff][2];
485 t_tau0i[ff](i) = t_coords[vi](i) - t_coords[v0](i);
486 t_tau0j[ff](i) = t_coords[vj](i) - t_coords[v0](i);
487 t_diff_ksi0i[ff] = get_diff_ksi(diffN, v0, vi);
488 t_diff_ksi0j[ff] = get_diff_ksi(diffN, v0, vj);
489 }
490
491 double psi_l[p + 1], psi_m[p + 1];
492 double diff_psi_l[3 * (p + 1)], diff_psi_m[3 * (p + 1)];
493 for (int ff = 0; ff != 4; ff++) {
494 const int v0 = faces_nodes[ff][0];
495 const int vi = faces_nodes[ff][1];
496 const int vj = faces_nodes[ff][2];
497 auto t_phi_v_f = getFTensor1FromPtr<3>(phi_v_f[ff]);
498 auto t_diff_phi_v_f = getFTensor2HVecFromPtr<3, 3>(diff_phi_v_f[ff]);
499 for (int ii = 0; ii < gdim; ii++) {
500 const int node_shift = 4 * ii;
501 const double n0 = N[node_shift + v0];
502 const double ni = N[node_shift + vi];
503 const double nj = N[node_shift + vj];
504 const double beta_f = n0 * ni * nj;
505 FTensor::Tensor1<double, 3> t_diff_beta_f;
506 t_diff_beta_f(i) = (ni * nj) * t_node_diff_ksi[v0](i) +
507 (n0 * nj) * t_node_diff_ksi[vi](i) +
508 (n0 * ni) * t_node_diff_ksi[vj](i);
509 const double ksi0i = get_ksi(n0, ni);
510 const double ksi0j = get_ksi(n0, nj);
511 base_polynomials(p, ksi0i, &t_diff_ksi0i[ff](0), psi_l, diff_psi_l, 3);
512 base_polynomials(p, ksi0j, &t_diff_ksi0j[ff](0), psi_m, diff_psi_m, 3);
513 FTensor::Tensor1<double, 3> t_diff_a;
514 int jj = 0;
515 for (int oo = 0; oo <= p - 3; oo++) {
516 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
517 auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
518 diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
519 for (int l = 0; l <= oo; l++) {
520 int m = oo - l;
521 if (m >= 0) {
522 auto t_diff_psi_m = FTensor::Tensor1<double, 3>{
523 diff_psi_m[m], diff_psi_m[p + 1 + m],
524 diff_psi_m[2 * p + 2 + m]};
525 const double a = beta_f * t_psi_l * psi_m[m];
526 t_phi_v_f(i) = a * t_tau0i[ff](i);
527 ++t_phi_v_f;
528 ++jj;
529 t_phi_v_f(i) = a * t_tau0j[ff](i);
530 ++t_phi_v_f;
531 ++jj;
532
533 t_diff_a(j) = (t_psi_l * psi_m[m]) * t_diff_beta_f(j) +
534 (beta_f * psi_m[m]) * t_diff_psi_l(j) +
535 (beta_f * t_psi_l) * t_diff_psi_m(j);
536 t_diff_phi_v_f(i, j) = t_diff_a(j) * t_tau0i[ff](i);
537 ++t_diff_phi_v_f;
538 t_diff_phi_v_f(i, j) = t_diff_a(j) * t_tau0j[ff](i);
539 ++t_diff_phi_v_f;
540
541 ++t_psi_l;
542 ++t_diff_psi_l;
543 }
544 }
545 }
546 if (jj != NBVOLUMETET_AINSWORTH_FACE_HDIV(p)) {
547 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
548 "wrong order %d != %d", jj,
549 NBVOLUMETET_AINSWORTH_FACE_HDIV(p));
550 }
551 }
552 }
553 MoFEMFunctionReturnHot(0);
554}
555
556MoFEMErrorCode MoFEM::Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(
557 int p, double *N, double *diffN, double *phi_v, double *diff_phi_v,
558 int gdim,
559 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
560 double *L, double *diffL,
561 const int dim)) {
562
563 MoFEMFunctionBeginHot;
564 if (p < 4)
565 MoFEMFunctionReturnHot(0);
566
567 FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
568 FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
569 FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
570 FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
571 FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
572
573 FTensor::Index<'i', 3> i;
574 FTensor::Index<'j', 3> j;
575 FTensor::Number<0> N0;
576 FTensor::Number<1> N1;
577 FTensor::Number<2> N2;
578
579 FTensor::Tensor1<double, 3> t_diff_ksi0i;
580 FTensor::Tensor1<double, 3> t_diff_ksi0j;
581 FTensor::Tensor1<double, 3> t_diff_ksi0k;
582
583 t_diff_ksi0i = get_diff_ksi(diffN, 0, 1);
584 t_diff_ksi0j = get_diff_ksi(diffN, 0, 2);
585 t_diff_ksi0k = get_diff_ksi(diffN, 0, 3);
586
587 double psi_l[p + 1];
588 double diff_psi_l[3 * (p + 1)];
589 double psi_m[p + 1];
590 double diff_psi_m[3 * (p + 1)];
591 double psi_n[p + 1];
592 double diff_psi_n[3 * (p + 1)];
593
594 auto t_phi_v = getFTensor1FromPtr<3>(phi_v);
595 auto t_diff_phi_v = getFTensor2HVecFromPtr<3, 3>(diff_phi_v);
596
597 FTensor::Tensor1<double, 3> t_diff_beta_v;
598 for (int ii = 0; ii < gdim; ii++) {
599 const int node_shift = ii * 4;
600 const double n0 = N[node_shift + 0];
601 const double ni = N[node_shift + 1];
602 const double nj = N[node_shift + 2];
603 const double nk = N[node_shift + 3];
604 const double ksi0i = get_ksi(n0, ni);
605 const double ksi0j = get_ksi(n0, nj);
606 const double ksi0k = get_ksi(n0, nk);
607 const double beta_v = n0 * ni * nj * nk;
608 t_diff_beta_v(i) = (ni * nj * nk) * t_node_diff_ksi[0](i) +
609 (n0 * nj * nk) * t_node_diff_ksi[1](i) +
610 (n0 * ni * nk) * t_node_diff_ksi[2](i) +
611 (n0 * ni * nj) * t_node_diff_ksi[3](i);
612 base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
613 base_polynomials(p, ksi0j, &t_diff_ksi0j(0), psi_m, diff_psi_m, 3);
614 base_polynomials(p, ksi0k, &t_diff_ksi0k(0), psi_n, diff_psi_n, 3);
615
616 FTensor::Tensor1<double, 3> t_diff_a;
617
618 int jj = 0;
619 for (int oo = 0; oo <= p - 4; oo++) {
620 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
621 auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
622 diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
623 for (int l = 0; l <= oo; l++) {
624 auto t_psi_m = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_m);
625 auto t_diff_psi_m = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
626 diff_psi_m, &diff_psi_m[p + 1], &diff_psi_m[2 * p + 2]);
627 for (int m = 0; (l + m) <= oo; m++) {
628 int n = oo - l - m;
629 if (n >= 0) {
630 FTensor::Tensor1<double, 3> t_diff_psi_n(diff_psi_n[n],
631 diff_psi_n[p + 1 + n],
632 diff_psi_n[2 * p + 2 + n]);
633 const double a = beta_v * t_psi_l * t_psi_m * psi_n[n];
634 t_phi_v(0) = a;
635 t_phi_v(1) = 0;
636 t_phi_v(2) = 0;
637 ++t_phi_v;
638 t_phi_v(0) = 0;
639 t_phi_v(1) = a;
640 t_phi_v(2) = 0;
641 ++t_phi_v;
642 t_phi_v(0) = 0;
643 t_phi_v(1) = 0;
644 t_phi_v(2) = a;
645 ++t_phi_v;
646 t_diff_a(j) = (t_psi_l * t_psi_m * psi_n[n]) * t_diff_beta_v(j) +
647 (beta_v * t_psi_m * psi_n[n]) * t_diff_psi_l(j) +
648 (beta_v * t_psi_l * psi_n[n]) * t_diff_psi_m(j) +
649 (beta_v * t_psi_l * t_psi_m) * t_diff_psi_n(j);
650 t_diff_phi_v(N0, j) = t_diff_a(j);
651 t_diff_phi_v(N1, j) = 0;
652 t_diff_phi_v(N2, j) = 0;
653 ++t_diff_phi_v;
654 t_diff_phi_v(N0, j) = 0;
655 t_diff_phi_v(N1, j) = t_diff_a(j);
656 t_diff_phi_v(N2, j) = 0;
657 ++t_diff_phi_v;
658 t_diff_phi_v(N0, j) = 0;
659 t_diff_phi_v(N1, j) = 0;
660 t_diff_phi_v(N2, j) = t_diff_a(j);
661 ++t_diff_phi_v;
662 ++jj;
663 }
664 ++t_psi_m;
665 ++t_diff_psi_m;
666 }
667 ++t_psi_l;
668 ++t_diff_psi_l;
669 }
670 }
671
672 if (3 * jj != NBVOLUMETET_AINSWORTH_VOLUME_HDIV(p)) {
673 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
674 "wrong order %d != %d", jj,
675 NBVOLUMETET_AINSWORTH_VOLUME_HDIV(p));
676 }
677 }
678
679 MoFEMFunctionReturnHot(0);
680}
681
682MoFEMErrorCode
683MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N,
684 double *diffN, double *phi_f,
685 double *diff_phi_f, int gdim, int nb) {
686 const int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
687 const int face_opposite_edges_node[] = {2, 0, 1};
688
689 MoFEMFunctionBeginHot;
690
691 FTensor::Index<'i', 3> i;
692 FTensor::Index<'j', 3> j;
693
694 FTensor::Tensor1<double, 3> t_cross[3];
695 FTensor::Tensor2<double, 3, 3> t_diff_cross(0., 0., 0., 0., 0., 0., 0., 0.,
696 0.);
697 FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
698 FTensor::Tensor1<double, 3> t_node_diff_sum_n0_n1;
699
700 const int i0 = faces_nodes[0];
701 const int i1 = faces_nodes[1];
702 const int i2 = faces_nodes[2];
703 const int o[] = {faces_nodes[face_opposite_edges_node[0]],
704 faces_nodes[face_opposite_edges_node[1]],
705 faces_nodes[face_opposite_edges_node[2]]};
706
707 FTensor::Tensor1<double, 3> t_diff_n0_p_n1;
708 FTensor::Tensor1<double, 3> t_diff_n0_p_n1_p_n2;
709
710 if (diff_phi_f) {
711 t_node_diff_ksi[0] =
712 FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
713 t_node_diff_ksi[1] =
714 FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
715 t_node_diff_ksi[2] =
716 FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
717 t_node_diff_ksi[3] =
718 FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
719 t_diff_cross(i, j) = 0;
720 for (int ee = 0; ee != 3; ee++) {
721 int ei0 = faces_nodes[face_edges_nodes[ee][0]];
722 int ei1 = faces_nodes[face_edges_nodes[ee][1]];
723 t_cross[ee](0) = t_node_diff_ksi[ei0](1) * t_node_diff_ksi[ei1](2) -
724 t_node_diff_ksi[ei0](2) * t_node_diff_ksi[ei1](1);
725 t_cross[ee](1) = t_node_diff_ksi[ei0](2) * t_node_diff_ksi[ei1](0) -
726 t_node_diff_ksi[ei0](0) * t_node_diff_ksi[ei1](2);
727 t_cross[ee](2) = t_node_diff_ksi[ei0](0) * t_node_diff_ksi[ei1](1) -
728 t_node_diff_ksi[ei0](1) * t_node_diff_ksi[ei1](0);
729 FTensor::Tensor1<double, 3> t_diff_o(
730 diffN[3 * o[ee] + 0], diffN[3 * o[ee] + 1], diffN[3 * o[ee] + 2]);
731 t_diff_cross(i, j) += t_cross[ee](i) * t_diff_o(j);
732 // cerr << t_cross[ee](0) << " " << t_cross[ee](1) << " " <<
733 // t_cross[ee](2) << endl;
734 }
735 // cerr << endl << endl;
736 t_diff_n0_p_n1(i) = t_node_diff_ksi[i0](i) + t_node_diff_ksi[i1](i);
737 t_diff_n0_p_n1_p_n2(i) = t_diff_n0_p_n1(i) + t_node_diff_ksi[i2](i);
738 } else {
739 for (int ee = 0; ee != 3; ee++) {
740 t_cross[ee](0) = 1;
741 t_cross[ee](1) = 0;
742 t_cross[ee](2) = 0;
743 }
744 }
745
746 FTensor::Tensor1<double *, 3> t_phi(&phi_f[HVEC0], &phi_f[HVEC1],
747 &phi_f[HVEC2], 3);
748 boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_phi_ptr;
749 if (diff_phi_f) {
750 t_diff_phi_ptr = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
751 new FTensor::Tensor2<double *, 3, 3>(
752 &diff_phi_f[HVEC0_0], &diff_phi_f[HVEC0_1], &diff_phi_f[HVEC0_2],
753 &diff_phi_f[HVEC1_0], &diff_phi_f[HVEC1_1], &diff_phi_f[HVEC1_2],
754 &diff_phi_f[HVEC2_0], &diff_phi_f[HVEC2_1], &diff_phi_f[HVEC2_2],
755 9));
756 }
757
758 double fi[p + 1], diff_fi[3 * p + 3];
759 double fj[p + 1], diff_fj[3 * p + 3];
760 double tmp_fj[p + 1], tmp_diff_fj[3 * p + 3];
761 for (int ii = 0; ii != gdim; ii++) {
762 const int shift = ii * nb;
763 double n0 = N[shift + i0];
764 double n1 = N[shift + i1];
765 double n2 = N[shift + i2];
766 double *diff_n1 = (diff_phi_f) ? &t_node_diff_ksi[i1](0) : NULL;
767 double *diff_n0_p_n1 = (diff_phi_f) ? &t_diff_n0_p_n1(0) : NULL;
768 ierr = Jacobi_polynomials(p, 0, n1, n0 + n1, diff_n1, diff_n0_p_n1, fi,
769 diff_phi_f ? diff_fi : NULL, 3);
770 CHKERRG(ierr);
771 for (int pp = 0; pp <= p; pp++) {
772 double *diff_n2 = (diff_phi_f) ? &t_node_diff_ksi[i2](0) : NULL;
773 double *diff_n0_p_n1_p_n2 = (diff_phi_f) ? &t_diff_n0_p_n1_p_n2(0) : NULL;
774 ierr = Jacobi_polynomials(pp, 2 * pp + 1, n2, n0 + n1 + n2, diff_n2,
775 diff_n0_p_n1_p_n2, tmp_fj,
776 diff_phi_f ? tmp_diff_fj : NULL, 3);
777 CHKERRG(ierr);
778 fj[pp] = tmp_fj[pp];
779 if (diff_phi_f) {
780 diff_fj[0 * (p + 1) + pp] = tmp_diff_fj[0 * (pp + 1) + pp];
781 diff_fj[1 * (p + 1) + pp] = tmp_diff_fj[1 * (pp + 1) + pp];
782 diff_fj[2 * (p + 1) + pp] = tmp_diff_fj[2 * (pp + 1) + pp];
783 }
784 }
785 double no0 = N[shift + o[0]];
786 double no1 = N[shift + o[1]];
787 double no2 = N[shift + o[2]];
788 FTensor::Tensor1<double, 3> base0;
789 base0(i) = no0 * t_cross[0](i) + no1 * t_cross[1](i) + no2 * t_cross[2](i);
790 int jj = 0;
791 for (int oo = 0; oo < p; oo++) {
792 FTensor::Tensor0<double *> t_fi(fi);
793 FTensor::Tensor1<double *, 3> t_diff_fi(&diff_fi[0], &diff_fi[p + 1],
794 &diff_fi[2 * p + 2], 1);
795 for (int ll = 0; ll <= oo; ll++) {
796 const int mm = oo - ll;
797 if (mm >= 0) {
798 const double a = t_fi * fj[mm];
799 // cerr << ll << " " << mm << " " << a << endl;
800 t_phi(i) = a * base0(i);
801 if (diff_phi_f) {
802 FTensor::Tensor1<double *, 3> t_diff_fj(
803 &diff_fj[0 + mm], &diff_fj[p + 1 + mm],
804 &diff_fj[2 * p + 2 + mm], 1);
805 (*t_diff_phi_ptr)(i, j) =
806 a * t_diff_cross(i, j) +
807 (t_diff_fi(j) * fj[mm] + t_fi * t_diff_fj(j)) * base0(i);
808 ++(*t_diff_phi_ptr);
809 ++t_diff_fi;
810 }
811 ++t_fi;
812 ++t_phi;
813 ++jj;
814 }
815 }
816 }
817 if (jj != NBFACETRI_DEMKOWICZ_HDIV(p)) {
818 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
819 "wrong number of base functions "
820 "jj!=NBFACETRI_DEMKOWICZ_HDIV(p) "
821 "%d!=%d",
822 jj, NBFACETRI_DEMKOWICZ_HDIV(p));
823 }
824 }
825
826 MoFEMFunctionReturnHot(0);
827}
828
829MoFEMErrorCode MoFEM::Hdiv_Demkowicz_Interior_MBTET(
830 int p, double *N, double *diffN, int p_f[], double *phi_f[4],
831 double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim) {
832
833 const int opposite_face_node[4] = {2, 0, 1, 3};
834 // list of zero node faces
835 const int znf[] = {0, 2, 3};
836 MoFEMFunctionBegin;
837
838 FTensor::Index<'i', 3> i;
839 FTensor::Index<'j', 3> j;
840
841 FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
842 t_node_diff_ksi[0] =
843 FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
844 t_node_diff_ksi[1] =
845 FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
846 t_node_diff_ksi[2] =
847 FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
848 t_node_diff_ksi[3] =
849 FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
850 FTensor::Tensor1<double, 3> t_m_node_diff_ksi[4];
851 for (int ff = 0; ff != 4; ++ff) {
852 t_m_node_diff_ksi[ff](i) = -t_node_diff_ksi[ff](i);
853 }
854
855 FTensor::Tensor1<double *, 3> t_phi_v(&phi_v[HVEC0], &phi_v[HVEC1],
856 &phi_v[HVEC2], 3);
857 FTensor::Tensor2<double *, 3, 3> t_diff_phi_v(
858 &diff_phi_v[HVEC0_0], &diff_phi_v[HVEC0_1], &diff_phi_v[HVEC0_2],
859 &diff_phi_v[HVEC1_0], &diff_phi_v[HVEC1_1], &diff_phi_v[HVEC1_2],
860 &diff_phi_v[HVEC2_0], &diff_phi_v[HVEC2_1], &diff_phi_v[HVEC2_2], 9);
861
862 MatrixDouble fk(3, p + 1), diff_fk(3, 3 * p + 3);
863
864 for (int ii = 0; ii != gdim; ii++) {
865 const int shift = 4 * ii;
866
867 for (int ff = 0; ff != 3; ff++) {
868 const int fff = znf[ff];
869 const int iO = opposite_face_node[fff];
870 const double nO = N[shift + iO];
871 for (int pp = 1; pp <= p; pp++) {
872 CHKERR IntegratedJacobi_polynomials(
873 pp, 2 * pp + 2, nO, 1 - nO, &t_node_diff_ksi[iO](0),
874 &t_m_node_diff_ksi[iO](0), &fk(ff, 0), &diff_fk(ff, 0), 3);
875 }
876 }
877
878 int jj = 0;
879 for (int oo = 2; oo <= p; oo++) {
880 for (int k = 1; k != oo; k++) {
881 int OO = oo - k;
882 if (OO >= 0) {
883 int s = NBFACETRI_DEMKOWICZ_HDIV(OO - 1);
884 // Note that we do faces 0,2,3, skipping 1. All the faces which have
885 // zero node in it.
886 int nb_dofs = NBFACETRI_DEMKOWICZ_HDIV(p_f[znf[0]]);
887 int sp[] = {ii * 3 * nb_dofs + 3 * s, ii * 3 * nb_dofs + 3 * s,
888 ii * 3 * nb_dofs + 3 * s};
889 FTensor::Tensor1<double *, 3> t_phi_f[] = {
890 FTensor::Tensor1<double *, 3>(&phi_f[znf[0]][sp[0] + HVEC0],
891 &phi_f[znf[0]][sp[0] + HVEC1],
892 &phi_f[znf[0]][sp[0] + HVEC2], 3),
893 FTensor::Tensor1<double *, 3>(&phi_f[znf[1]][sp[1] + HVEC0],
894 &phi_f[znf[1]][sp[1] + HVEC1],
895 &phi_f[znf[1]][sp[1] + HVEC2], 3),
896 FTensor::Tensor1<double *, 3>(&phi_f[znf[2]][sp[2] + HVEC0],
897 &phi_f[znf[2]][sp[2] + HVEC1],
898 &phi_f[znf[2]][sp[2] + HVEC2], 3)};
899 int sdp[] = {ii * 9 * nb_dofs + 9 * s, ii * 9 * nb_dofs + 9 * s,
900 ii * 9 * nb_dofs + 9 * s};
901 FTensor::Tensor2<double *, 3, 3> t_diff_phi_f[] = {
902 FTensor::Tensor2<double *, 3, 3>(
903 &diff_phi_f[znf[0]][sdp[0] + HVEC0_0],
904 &diff_phi_f[znf[0]][sdp[0] + HVEC0_1],
905 &diff_phi_f[znf[0]][sdp[0] + HVEC0_2],
906 &diff_phi_f[znf[0]][sdp[0] + HVEC1_0],
907 &diff_phi_f[znf[0]][sdp[0] + HVEC1_1],
908 &diff_phi_f[znf[0]][sdp[0] + HVEC1_2],
909 &diff_phi_f[znf[0]][sdp[0] + HVEC2_0],
910 &diff_phi_f[znf[0]][sdp[0] + HVEC2_1],
911 &diff_phi_f[znf[0]][sdp[0] + HVEC2_2], 9),
912 FTensor::Tensor2<double *, 3, 3>(
913 &diff_phi_f[znf[1]][sdp[1] + HVEC0_0],
914 &diff_phi_f[znf[1]][sdp[1] + HVEC0_1],
915 &diff_phi_f[znf[1]][sdp[1] + HVEC0_2],
916 &diff_phi_f[znf[1]][sdp[1] + HVEC1_0],
917 &diff_phi_f[znf[1]][sdp[1] + HVEC1_1],
918 &diff_phi_f[znf[1]][sdp[1] + HVEC1_2],
919 &diff_phi_f[znf[1]][sdp[1] + HVEC2_0],
920 &diff_phi_f[znf[1]][sdp[1] + HVEC2_1],
921 &diff_phi_f[znf[1]][sdp[1] + HVEC2_2], 9),
922 FTensor::Tensor2<double *, 3, 3>(
923 &diff_phi_f[znf[2]][sdp[2] + HVEC0_0],
924 &diff_phi_f[znf[2]][sdp[2] + HVEC0_1],
925 &diff_phi_f[znf[2]][sdp[2] + HVEC0_2],
926 &diff_phi_f[znf[2]][sdp[2] + HVEC1_0],
927 &diff_phi_f[znf[2]][sdp[2] + HVEC1_1],
928 &diff_phi_f[znf[2]][sdp[2] + HVEC1_2],
929 &diff_phi_f[znf[2]][sdp[2] + HVEC2_0],
930 &diff_phi_f[znf[2]][sdp[2] + HVEC2_1],
931 &diff_phi_f[znf[2]][sdp[2] + HVEC2_2], 9)};
932 for (int ij = s; ij != NBFACETRI_DEMKOWICZ_HDIV(OO); ij++) {
933 for (int ff = 0; ff != 3; ff++) {
934 FTensor::Tensor1<double, 3> t_diff_fk(diff_fk(ff, 0 * p + k - 1),
935 diff_fk(ff, 1 * p + k - 1),
936 diff_fk(ff, 2 * p + k - 1));
937 t_phi_v(i) = fk(ff, k - 1) * t_phi_f[ff](i);
938 t_diff_phi_v(i, j) = t_diff_fk(j) * t_phi_f[ff](i) +
939 fk(ff, k - 1) * t_diff_phi_f[ff](i, j);
940 ++t_phi_v;
941 ++t_diff_phi_v;
942 ++t_phi_f[ff];
943 ++t_diff_phi_f[ff];
944 ++jj;
945 }
946 }
947 }
948 }
949 }
950 if (jj != NBVOLUMETET_DEMKOWICZ_HDIV(p)) {
951 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
952 "wrong number of base functions "
953 "jj!=NBVOLUMETET_DEMKOWICZ_HDIV(p) "
954 "%d!=%d",
955 jj, NBVOLUMETET_DEMKOWICZ_HDIV(p));
956 }
957 }
958
959 MoFEMFunctionReturn(0);
960}
get_ksi
static double get_ksi(const double n0, const double n1)
Definition Hdiv.cpp:47
get_diff_ksi
static FTensor::Tensor1< double, 3 > get_diff_ksi(const double *diffN, const int n0, const int n1)
Definition Hdiv.cpp:70
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition Templates.hpp:2561
a
constexpr double a
Definition approx_sphere.cpp:30
IntegratedJacobi_polynomials
PetscErrorCode IntegratedJacobi_polynomials(int p, double alpha, double x, double t, double *diff_x, double *diff_t, double *L, double *diffL, const int dim)
Calculate integrated Jacobi approximation basis.
Definition base_functions.c:84
Jacobi_polynomials
PetscErrorCode Jacobi_polynomials(int p, double alpha, double x, double t, double *diff_x, double *diff_t, double *L, double *diffL, const int dim)
Calculate Jacobi approximation basis.
Definition base_functions.c:17
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition definitions.h:499
HVEC0
@ HVEC0
Definition definitions.h:199
HVEC1
@ HVEC1
Definition definitions.h:199
HVEC2
@ HVEC2
Definition definitions.h:199
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
HVEC1_1
@ HVEC1_1
Definition definitions.h:209
HVEC0_1
@ HVEC0_1
Definition definitions.h:208
HVEC1_0
@ HVEC1_0
Definition definitions.h:206
HVEC2_1
@ HVEC2_1
Definition definitions.h:210
HVEC1_2
@ HVEC1_2
Definition definitions.h:212
HVEC2_2
@ HVEC2_2
Definition definitions.h:213
HVEC2_0
@ HVEC2_0
Definition definitions.h:207
HVEC0_2
@ HVEC0_2
Definition definitions.h:211
HVEC0_0
@ HVEC0_0
Definition definitions.h:205
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
NBVOLUMETET_DEMKOWICZ_HDIV
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:140
NBVOLUMETET_AINSWORTH_FACE_HDIV
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:134
NBFACETRI_DEMKOWICZ_HDIV
#define NBFACETRI_DEMKOWICZ_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:139
NBFACETRI_AINSWORTH_FACE_HDIV
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:131
NBVOLUMETET_AINSWORTH_VOLUME_HDIV
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:135
NBFACETRI_AINSWORTH_EDGE_HDIV
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:130
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
lambda
static double lambda
Definition incompressible_elasticity.cpp:181
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:23
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition Exceptions.hpp:76
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[], double *diff_phi_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
Definition Hdiv.cpp:201
MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition Hdiv.cpp:683
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3], double *diff_phi_f_e[4][3], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
Definition Hdiv.cpp:24
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3], double *diff_phi_f_e[3], int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
Definition Hdiv.cpp:97
MoFEM::getFTensor2HVecFromPtr< 3, 3 >
auto getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
Definition Templates.hpp:641
MoFEM::Hdiv_Demkowicz_Interior_MBTET
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
Definition Hdiv.cpp:829
MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition Hdiv.cpp:450
MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
Definition Hdiv.cpp:224
MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_e[6], double *diff_phi_v_e[6], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base function, Edge-based interior (volume) functions by Ainsworth .
Definition Hdiv.cpp:357
MoFEM::Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Interior bubble functions by Ainsworth .
Definition Hdiv.cpp:556
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
N
const int N
Definition speed_test.cpp:3
MoFEM::AinsworthHdivOptions::edge_coordinate
static constexpr EdgeCoordinate edge_coordinate
Selected edge coordinate variant.
Definition Hdiv.hpp:50
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv
static boost::function< int(int)> broken_nbvolumetet_edge_hdiv
Definition Hdiv.hpp:27
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_face_hdiv
static boost::function< int(int)> broken_nbvolumetet_face_hdiv
Definition Hdiv.hpp:28
MoFEM::AinsworthOrderHooks::broken_nbfacetri_face_hdiv
static boost::function< int(int)> broken_nbfacetri_face_hdiv
Definition Hdiv.hpp:26
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv
static boost::function< int(int)> broken_nbvolumetet_volume_hdiv
Definition Hdiv.hpp:29
MoFEM::AinsworthOrderHooks::broken_nbfacetri_edge_hdiv
static boost::function< int(int)> broken_nbfacetri_edge_hdiv
Definition Hdiv.hpp:25
Generated by Doxygen 1.9.8 and hosted at University of Glasgow