v0.14.0
Loading...
Searching...
No Matches
EdgeQuadHexPolynomials.cpp
Go to the documentation of this file.
1/** \file EdgeQuadHexPolynomials.cpp
2
3 \brief Implementation of hierarchical Edge, Quad, and Hex shape functions of
4 type H1, Hcurl, Hdiv
5*/
6
7using namespace MoFEM;
8
9namespace DemkowiczHexAndQuad {
10
11MoFEMErrorCode monom_ordering(int *perm, int p, int q, int r = 0) {
12 MoFEMFunctionBeginHot;
13
14 if (r > 0) {
15
16 for (int m = 0; m != std::max(std::max(p, q), r) + 1; ++m) {
17
18 const int i = std::min(m, p);
19 const int j = std::min(m, q);
20 const int k = std::min(m, r);
21
22 if (i == m)
23 for (int jj = 0; jj != j; ++jj) {
24 for (int kk = 0; kk != k; ++kk) {
25 *(perm++) = i;
26 *(perm++) = jj;
27 *(perm++) = kk;
28 }
29 }
30
31 if (j == m)
32 for (int ii = 0; ii != i; ++ii) {
33 for (int kk = 0; kk != k; ++kk) {
34 *(perm++) = ii;
35 *(perm++) = j;
36 *(perm++) = kk;
37 }
38 }
39
40 if (k == m)
41 for (int ii = 0; ii != i; ++ii) {
42 for (int jj = 0; jj != j; ++jj) {
43 *(perm++) = ii;
44 *(perm++) = jj;
45 *(perm++) = k;
46 }
47 }
48
49 if (j == m || k == m)
50 for (int ii = 0; ii != i; ++ii) {
51 *(perm++) = ii;
52 *(perm++) = j;
53 *(perm++) = k;
54 }
55
56 if (i == m || k == m)
57 for (int jj = 0; jj != j; ++jj) {
58 *(perm++) = i;
59 *(perm++) = jj;
60 *(perm++) = k;
61 }
62
63 if (i == m || j == m)
64 for (int kk = 0; kk != k; ++kk) {
65 *(perm++) = i;
66 *(perm++) = j;
67 *(perm++) = kk;
68 }
69
70 *(perm++) = i;
71 *(perm++) = j;
72 *(perm++) = k;
73 }
74 } else {
75
76 for (int m = 0; m != std::max(p, q) + 1; ++m) {
77
78 const int i = std::min(m, p);
79 const int j = std::min(m, q);
80
81 if (j == m)
82 for (int ii = 0; ii != i; ++ii) {
83 *(perm++) = ii;
84 *(perm++) = j;
85 *(perm++) = 0;
86 }
87
88 if (i == m)
89 for (int jj = 0; jj != j; ++jj) {
90 *(perm++) = i;
91 *(perm++) = jj;
92 *(perm++) = 0;
93 }
94
95 *(perm++) = i;
96 *(perm++) = j;
97 *(perm++) = 0;
98 }
99 }
100 MoFEMFunctionReturnHot(0);
101}
102
103static inline void get_ksi_hex(int shift, double *N, double *N_diff,
104 double ksi[3], double diff_ksi[3][3]) {
105
106 constexpr std::array<size_t, 4> ksi_nodes[2][3] = {
107
108 {{1, 2, 6, 5}, {3, 2, 6, 7}, {4, 5, 6, 7}},
109
110 {{0, 3, 7, 4}, {0, 1, 5, 4}, {0, 1, 2, 3}}
111
112 };
113
114 for (size_t i = 0; i != 3; ++i) {
115 for (auto n : ksi_nodes[0][i])
116 ksi[i] += N[shift + n];
117 for (auto n : ksi_nodes[1][i])
118 ksi[i] -= N[shift + n];
119 }
120
121 for (size_t i = 0; i != 3; ++i) {
122 for (auto d = 0; d != 3; ++d) {
123 for (auto n : ksi_nodes[0][i]) {
124
125 diff_ksi[i][d] += N_diff[3 * (shift + n) + d];
126 }
127 }
128 for (auto n : ksi_nodes[1][i]) {
129 for (auto d = 0; d != 3; ++d) {
130 diff_ksi[i][d] -= N_diff[3 * (shift + n) + d];
131 }
132 }
133 }
134}
135
136} // namespace DemkowiczHexAndQuad
137
138/*
139 0 *------------* 1
140*/
141MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::H1_BubbleShapeFunctions_ONSEGMENT(
142 int p, double *N, double *diffN, double *bubbleN, double *diff_bubbleN,
143 int nb_integration_pts) {
144 MoFEMFunctionBeginHot;
145 const int nb_dofs = (p - 1);
146 if (nb_dofs > 0) {
147 constexpr int n0 = 0;
148 constexpr int n1 = 1;
149 double diff_mu = diffN[n1] - diffN[n0];
150 for (int q = 0; q != nb_integration_pts; q++) {
151 int shift = 2 * q;
152 double mu = N[shift + n1] - N[shift + n0];
153 double L[p + 2];
154 double diffL[p + 2];
155 CHKERR Lobatto_polynomials(p + 1, mu, &diff_mu, L, diffL, 1);
156 int qd_shift = nb_dofs * q;
157 for (int n = 0; n != nb_dofs; n++) {
158 bubbleN[qd_shift + n] = L[n + 2];
159 diff_bubbleN[qd_shift + n] = diffL[n + 2];
160 }
161 }
162 }
163 MoFEMFunctionReturnHot(0);
164}
165
166MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::L2_ShapeFunctions_ONSEGMENT(
167 int p, double *N, double *diffN, double *funN, double *funDiffN,
168 int nb_integration_pts) {
169 MoFEMFunctionBeginHot;
170 const int nb_dofs = p + 1;
171 if (nb_dofs > 0) {
172 constexpr int n0 = 0;
173 constexpr int n1 = 1;
174 double diff_mu = diffN[n1] - diffN[n0];
175 for (int q = 0; q != nb_integration_pts; q++) {
176 int shift = 2 * q;
177 double mu = N[shift + n1] - N[shift + n0];
178 double L[p + 1];
179 double diffL[p + 1];
180 CHKERR Legendre_polynomials(p, mu, &diff_mu, L, diffL, 1);
181 int qd_shift = (p + 1) * q;
182 for (int n = 0; n <= p; n++) {
183 funN[qd_shift + n] = L[n];
184 funDiffN[qd_shift + n] = diffL[n];
185 }
186 }
187 }
188 MoFEMFunctionReturnHot(0);
189}
190
191/*
192 Quads
193 3-------2------2
194 | | eta
195 | | ^
196 3 1 |
197 | | |
198 | | 0----- > ksi
199 0-------0------1
200*/
201
202MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::H1_EdgeShapeFunctions_ONQUAD(
203 int *sense, int *p, double *N, double *diffN, double *edgeN[4],
204 double *diff_edgeN[4], int nb_integration_pts) {
205 MoFEMFunctionBeginHot;
206
207 constexpr int n0 = 0;
208 constexpr int n1 = 1;
209 constexpr int n2 = 2;
210 constexpr int n3 = 3;
211
212 for (int q = 0; q != nb_integration_pts; q++) {
213 const int shift = 4 * q;
214 const double shape0 = N[shift + n0];
215 const double shape1 = N[shift + n1];
216 const double shape2 = N[shift + n2];
217 const double shape3 = N[shift + n3];
218 const double ksi01 = (shape1 + shape2 - shape0 - shape3) * sense[n0];
219 const double ksi12 = (shape2 + shape3 - shape1 - shape0) * sense[n1];
220 const double ksi23 = (shape3 + shape0 - shape2 - shape1) * sense[n2];
221 const double ksi30 = (shape0 + shape1 - shape3 - shape2) * sense[n3];
222
223 const double mu[] = {ksi01, ksi12, ksi23, ksi30};
224 const double mu_const[] = {shape0 + shape1, shape1 + shape2,
225 shape2 + shape3, shape3 + shape0};
226
227 const int diff_shift = 2 * shift;
228 double diff_mu_const[4][2];
229 double diff_mu[4][2];
230 for (int d = 0; d != 2; d++) {
231 const double diff_shape0 = diffN[diff_shift + 2 * n0 + d];
232 const double diff_shape1 = diffN[diff_shift + 2 * n1 + d];
233 const double diff_shape2 = diffN[diff_shift + 2 * n2 + d];
234 const double diff_shape3 = diffN[diff_shift + 2 * n3 + d];
235 diff_mu_const[0][d] = diff_shape0 + diff_shape1;
236 diff_mu_const[1][d] = diff_shape1 + diff_shape2;
237 diff_mu_const[2][d] = diff_shape2 + diff_shape3;
238 diff_mu_const[3][d] = diff_shape3 + diff_shape0;
239
240 const double diff_ksi01 =
241 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) * sense[0];
242 const double diff_ksi12 =
243 (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) * sense[1];
244 const double diff_ksi23 =
245 (diff_shape3 + diff_shape0 - diff_shape2 - diff_shape1) * sense[2];
246 const double diff_ksi30 =
247 (diff_shape0 + diff_shape1 - diff_shape3 - diff_shape2) * sense[3];
248 diff_mu[0][d] = diff_ksi01;
249 diff_mu[1][d] = diff_ksi12;
250 diff_mu[2][d] = diff_ksi23;
251 diff_mu[3][d] = diff_ksi30;
252 }
253
254 for (int e = 0; e != 4; e++) {
255 const int nb_dofs = p[e] - 1;
256 if (nb_dofs > 0) {
257 double L[p[e] + 2];
258 double diffL[2 * (p[e] + 2)];
259 CHKERR Lobatto_polynomials(p[e] + 1, mu[e], diff_mu[e], L, diffL, 2);
260 int qd_shift = (p[e] - 1) * q;
261 for (int n = 0; n != p[e] - 1; ++n) {
262 edgeN[e][qd_shift + n] = mu_const[e] * L[n + 2];
263 for (int d = 0; d != 2; ++d) {
264 diff_edgeN[e][2 * qd_shift + 2 * n + d] =
265 mu_const[e] * diffL[d * (p[e] + 2) + n + 2] +
266 diff_mu_const[e][d] * L[n + 2];
267 }
268 }
269 }
270 }
271 }
272 MoFEMFunctionReturnHot(0);
273}
274
275MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::H1_FaceShapeFunctions_ONQUAD(
276 int *faces_nodes, int *p, double *N, double *diffN, double *faceN,
277 double *diff_faceN, int nb_integration_pts) {
278
279 if (p[0] > 1 && p[1] > 1) {
280 const int nb_dofs = (p[0] - 1) * (p[1] - 1);
281
282 const int n0 = faces_nodes[0];
283 const int n1 = faces_nodes[1];
284 const int n2 = faces_nodes[2];
285 const int n3 = faces_nodes[3];
286
287 MoFEMFunctionBeginHot;
288 int permute[(p[0] - 1) * (p[1] - 1)][3];
289 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0], p[0] - 2,
290 p[1] - 2);
291 for (int q = 0; q != nb_integration_pts; q++) {
292
293 const int shift = 4 * q;
294 const double shape0 = N[shift + n0];
295 const double shape1 = N[shift + n1];
296 const double shape2 = N[shift + n2];
297 const double shape3 = N[shift + n3];
298 const double ksi01 = (shape1 + shape2 - shape0 - shape3);
299 const double ksi12 = (shape2 + shape3 - shape1 - shape0);
300
301 const int diff_shift = 2 * shift;
302 double diff_ksi01[2], diff_ksi12[2];
303 for (int d = 0; d != 2; d++) {
304 const double diff_shape0 = diffN[diff_shift + 2 * n0 + d];
305 const double diff_shape1 = diffN[diff_shift + 2 * n1 + d];
306 const double diff_shape2 = diffN[diff_shift + 2 * n2 + d];
307 const double diff_shape3 = diffN[diff_shift + 2 * n3 + d];
308 diff_ksi01[d] = (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3);
309 diff_ksi12[d] = (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0);
310 }
311
312 double L01[p[0] + 2];
313 double diffL01[2 * (p[0] + 2)];
314 CHKERR Lobatto_polynomials(p[0] + 1, ksi01, diff_ksi01, L01, diffL01, 2);
315 double L12[p[1] + 2];
316 double diffL12[2 * (p[1] + 2)];
317 CHKERR Lobatto_polynomials(p[1] + 1, ksi12, diff_ksi12, L12, diffL12, 2);
318
319 int qd_shift = nb_dofs * q;
320 for (int n = 0; n != nb_dofs; ++n) {
321 int s1 = permute[n][0];
322 int s2 = permute[n][1];
323 faceN[qd_shift + n] = L01[s1 + 2] * L12[s2 + 2];
324 for (int d = 0; d != 2; ++d) {
325 diff_faceN[2 * (qd_shift + n) + d] =
326 diffL01[d * (p[0] + 2) + s1 + 2] * L12[s2 + 2] +
327 L01[s1 + 2] * diffL12[d * (p[1] + 2) + s2 + 2];
328 }
329 }
330 }
331 }
332 MoFEMFunctionReturnHot(0);
333}
334
335MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::L2_FaceShapeFunctions_ONQUAD(
336 int *p, double *N, double *diffN, double *faceN, double *diff_faceN,
337 int nb_integration_pts) {
338 MoFEMFunctionBeginHot;
339 const int nb_dofs = (p[0] + 1) * (p[1] + 1);
340 if (nb_dofs > 0) {
341
342 int permute[nb_dofs][3];
343 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0], p[0], p[1]);
344
345 constexpr int n0 = 0;
346 constexpr int n1 = 1;
347 constexpr int n2 = 2;
348 constexpr int n3 = 3;
349
350 for (int q = 0; q != nb_integration_pts; q++) {
351 const int shift = 4 * q;
352 const double shape0 = N[shift + n0];
353 const double shape1 = N[shift + n1];
354 const double shape2 = N[shift + n2];
355 const double shape3 = N[shift + n3];
356 const double ksi01 = (shape1 + shape2 - shape0 - shape3);
357 const double ksi12 = (shape2 + shape3 - shape1 - shape0);
358
359 const int diff_shift = 2 * shift;
360 double diff_ksi01[2], diff_ksi12[2];
361 for (int d = 0; d != 2; d++) {
362 const double diff_shape0 = diffN[diff_shift + 2 * n0 + d];
363 const double diff_shape1 = diffN[diff_shift + 2 * n1 + d];
364 const double diff_shape2 = diffN[diff_shift + 2 * n2 + d];
365 const double diff_shape3 = diffN[diff_shift + 2 * n3 + d];
366 diff_ksi01[d] = (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3);
367 diff_ksi12[d] = (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0);
368 }
369
370 double L01[p[0] + 2];
371 double diffL01[2 * (p[0] + 2)];
372 CHKERR Legendre_polynomials(p[0] + 1, ksi01, diff_ksi01, L01, diffL01, 2);
373 double L12[p[1] + 2];
374 double diffL12[2 * (p[1] + 2)];
375 CHKERR Legendre_polynomials(p[1] + 1, ksi12, diff_ksi12, L12, diffL12, 2);
376
377 int qd_shift = nb_dofs * q;
378 for (int n = 0; n != nb_dofs; ++n) {
379 int s1 = permute[n][0];
380 int s2 = permute[n][1];
381 faceN[qd_shift + n] = L01[s1] * L12[s2];
382 for (int d = 0; d != 2; ++d) {
383 diff_faceN[2 * (qd_shift + n) + d] =
384 diffL01[d * (p[0] + 2) + s1] * L12[s2] +
385 L01[s1] * diffL12[d * (p[1] + 2) + s2];
386 }
387 }
388 }
389 }
390 MoFEMFunctionReturnHot(0);
391}
392
393MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hcurl_EdgeShapeFunctions_ONQUAD(
394 int *sense, int *p, double *N, double *diffN, double *edgeN[4],
395 double *diff_edgeN[4], int nb_integration_pts) {
396 MoFEMFunctionBeginHot;
397
398 constexpr int n0 = 0;
399 constexpr int n1 = 1;
400 constexpr int n2 = 2;
401 constexpr int n3 = 3;
402
403 for (int q = 0; q != nb_integration_pts; q++) {
404
405 const int shift = 4 * q;
406 const double shape0 = N[shift + n0];
407 const double shape1 = N[shift + n1];
408 const double shape2 = N[shift + n2];
409 const double shape3 = N[shift + n3];
410 const double ksi01 = (shape1 + shape2 - shape0 - shape3) * sense[n0];
411 const double ksi12 = (shape2 + shape3 - shape1 - shape0) * sense[n1];
412 const double ksi23 = (shape3 + shape0 - shape2 - shape1) * sense[n2];
413 const double ksi30 = (shape0 + shape1 - shape3 - shape2) * sense[n3];
414
415 const double mu[] = {ksi01, ksi12, ksi23, ksi30};
416 const double mu_const[] = {shape0 + shape1, shape1 + shape2,
417 shape2 + shape3, shape3 + shape0};
418
419 const int diff_shift = 2 * shift;
420 double diff_mu_const[4][2];
421 double diff_mu[4][2];
422 for (int d = 0; d != 2; d++) {
423 const double diff_shape0 = diffN[diff_shift + 2 * n0 + d];
424 const double diff_shape1 = diffN[diff_shift + 2 * n1 + d];
425 const double diff_shape2 = diffN[diff_shift + 2 * n2 + d];
426 const double diff_shape3 = diffN[diff_shift + 2 * n3 + d];
427 diff_mu_const[0][d] = diff_shape0 + diff_shape1;
428 diff_mu_const[1][d] = diff_shape1 + diff_shape2;
429 diff_mu_const[2][d] = diff_shape2 + diff_shape3;
430 diff_mu_const[3][d] = diff_shape3 + diff_shape0;
431
432 const double diff_ksi01 =
433 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) * sense[n0];
434 const double diff_ksi12 =
435 (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) * sense[n1];
436 const double diff_ksi23 =
437 (diff_shape3 + diff_shape0 - diff_shape2 - diff_shape1) * sense[n2];
438 const double diff_ksi30 =
439 (diff_shape0 + diff_shape1 - diff_shape3 - diff_shape2) * sense[n3];
440 diff_mu[0][d] = diff_ksi01;
441 diff_mu[1][d] = diff_ksi12;
442 diff_mu[2][d] = diff_ksi23;
443 diff_mu[3][d] = diff_ksi30;
444 }
445
446 for (int e = 0; e != 4; e++) {
447
448 if (p[e] > 0) {
449 double L[p[e]];
450 double diffL[2 * p[e]];
451
452 CHKERR Legendre_polynomials(p[e] - 1, mu[e], diff_mu[e], L, diffL, 2);
453
454 int qd_shift = p[e] * q;
455 double *t_n_ptr = &edgeN[e][3 * qd_shift];
456 double *t_diff_n_ptr = &diff_edgeN[e][3 * 2 * qd_shift];
457 auto t_n = getFTensor1FromPtr<3>(t_n_ptr);
458 auto t_diff_n = getFTensor2HVecFromPtr<3, 2>(t_diff_n_ptr);
459
460 for (int n = 0; n != p[e]; ++n) {
461 const double a = mu_const[e] * L[n];
462 const double d_a[] = {
463 diff_mu_const[e][0] * L[n] + mu_const[e] * diffL[0 * p[e] + n],
464 diff_mu_const[e][1] * L[n] + mu_const[e] * diffL[1 * p[e] + n]};
465
466 for (int d = 0; d != 2; ++d) {
467 t_n(d) = (a / 2) * diff_mu[e][d];
468 for (int j = 0; j != 2; ++j) {
469 t_diff_n(d, j) = (d_a[j] / 2) * diff_mu[e][d];
470 }
471 }
472 t_n(2) = 0;
473 t_diff_n(2, 0) = 0;
474 t_diff_n(2, 1) = 0;
475 ++t_n;
476 ++t_diff_n;
477 }
478 }
479 }
480 }
481 MoFEMFunctionReturnHot(0);
482}
483
484MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hcurl_FaceShapeFunctions_ONQUAD(
485 int *face_nodes, int *p, double *N, double *diffN, double *faceN[],
486 double *diff_faceN[], int nb_integration_pts) {
487 MoFEMFunctionBeginHot;
488
489 const int pq[2] = {p[0], p[1]};
490 const int qp[2] = {p[1], p[0]};
491
492 const int nb_dofs_fm0 = pq[0] * (qp[1] - 1);
493 const int nb_dofs_fm1 = (pq[0] - 1) * qp[1];
494 int permute_fm0[3 * nb_dofs_fm0];
495 int permute_fm1[3 * nb_dofs_fm1];
496
497 std::array<int *, 2> permute = {permute_fm0, permute_fm1};
498 for (int fm = 0; fm != 2; ++fm) {
499 const int pp = pq[fm];
500 const int qq = qp[fm];
501 CHKERR ::DemkowiczHexAndQuad::monom_ordering(permute[fm], qq - 1, pp - 2);
502 }
503
504 const int n0 = face_nodes[0];
505 const int n1 = face_nodes[1];
506 const int n2 = face_nodes[2];
507 const int n3 = face_nodes[3];
508
509 for (int q = 0; q != nb_integration_pts; q++) {
510
511 const int shift = 4 * q;
512 const double shape0 = N[shift + n0];
513 const double shape1 = N[shift + n1];
514 const double shape2 = N[shift + n2];
515 const double shape3 = N[shift + n3];
516 const double ksi01 = (shape1 + shape2 - shape0 - shape3);
517 const double ksi12 = (shape2 + shape3 - shape1 - shape0);
518
519 const int diff_shift = 2 * shift;
520 double diff_ksi01[2], diff_ksi12[2];
521 for (int d = 0; d != 2; d++) {
522 const double diff_shape0 = diffN[diff_shift + 2 * n0 + d];
523 const double diff_shape1 = diffN[diff_shift + 2 * n1 + d];
524 const double diff_shape2 = diffN[diff_shift + 2 * n2 + d];
525 const double diff_shape3 = diffN[diff_shift + 2 * n3 + d];
526 diff_ksi01[d] = (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3);
527 diff_ksi12[d] = (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0);
528 }
529
530 double mu_ksi_eta[2] = {ksi01, ksi12};
531 double mu_eta_ksi[2] = {ksi12, ksi01};
532
533 double *diff_ksi_eta[2] = {diff_ksi01, diff_ksi12};
534 double *diff_eta_ksi[2] = {diff_ksi12, diff_ksi01};
535
536 for (int family = 0; family != 2; family++) {
537 const int pp = pq[family];
538 const int qq = qp[family];
539 const int nb_dofs = (pp - 1) * qq;
540
541 if (nb_dofs > 0) {
542
543 double zeta[pp + 2];
544 double diff_zeta[2 * (pp + 2)];
545 CHKERR Lobatto_polynomials(pp + 1, mu_ksi_eta[family],
546 diff_ksi_eta[family], zeta, diff_zeta, 2);
547
548 double eta[qq];
549 double diff_eta[2 * qq];
550 CHKERR Legendre_polynomials(qq - 1, mu_eta_ksi[family],
551 diff_eta_ksi[family], eta, diff_eta, 2);
552
553 const int qd_shift = nb_dofs * q;
554 double *t_n_ptr = &faceN[family][3 * qd_shift];
555 double *t_diff_n_ptr = &diff_faceN[family][6 * qd_shift];
556 auto t_n = getFTensor1FromPtr<3>(t_n_ptr);
557 auto t_diff_n = getFTensor2HVecFromPtr<3, 2>(t_diff_n_ptr);
558
559 for (int n = 0; n != nb_dofs; n++) {
560 int i = permute[family][3 * n + 0];
561 int j = permute[family][3 * n + 1];
562
563 const double z = zeta[j + 2];
564 const double e = eta[i];
565 const double a = z * e;
566 const double d_a[] = {
567
568 diff_zeta[0 * (pp + 2) + j + 2] * e + z * diff_eta[0 * qq + i],
569
570 diff_zeta[1 * (pp + 2) + j + 2] * e + z * diff_eta[1 * qq + i]};
571
572 for (int d = 0; d != 2; ++d) {
573 t_n(d) = a * diff_eta_ksi[family][d];
574 for (int m = 0; m != 2; ++m) {
575 t_diff_n(d, m) = d_a[m] * diff_eta_ksi[family][d];
576 }
577 }
578 t_n(2) = 0;
579 t_diff_n(2, 0) = 0;
580 t_diff_n(2, 1) = 0;
581 ++t_n;
582 ++t_diff_n;
583 }
584 }
585 }
586 }
587 MoFEMFunctionReturnHot(0);
588}
589
590MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hdiv_FaceShapeFunctions_ONQUAD(
591 int *face_nodes, int *p, double *N, double *diffN, double *faceN,
592 double *diff_faceN, int nb_integration_pts) {
593 MoFEMFunctionBeginHot;
594
595 const int nb_dofs = (p[0] * p[1]);
596
597 if (nb_dofs > 0) {
598
599 int permute[nb_dofs][3];
600 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0], p[0] - 1,
601 p[1] - 1);
602
603 const int n0 = face_nodes[0];
604 const int n1 = face_nodes[1];
605 const int n2 = face_nodes[2];
606 const int n3 = face_nodes[3];
607
608 FTensor::Index<'i', 3> i;
609 FTensor::Index<'j', 3> j;
610 FTensor::Index<'k', 3> k;
611
612 auto t_n = getFTensor1FromPtr<3>(faceN);
613 auto t_diff_n = getFTensor2HVecFromPtr<3, 2>(diff_faceN);
614
615 for (int q = 0; q != nb_integration_pts; q++) {
616
617 const int shift = 4 * q;
618 const double shape0 = N[shift + n0];
619 const double shape1 = N[shift + n1];
620 const double shape2 = N[shift + n2];
621 const double shape3 = N[shift + n3];
622 const double ksi01 = (shape1 + shape2 - shape0 - shape3);
623 const double ksi12 = (shape2 + shape3 - shape1 - shape0);
624
625 const int diff_shift = 2 * shift;
626 FTensor::Tensor1<double, 3> t_diff_ksi01;
627 FTensor::Tensor1<double, 3> t_diff_ksi12;
628 for (int d = 0; d != 2; d++) {
629 const double diff_shape0 = diffN[diff_shift + 2 * n0 + d];
630 const double diff_shape1 = diffN[diff_shift + 2 * n1 + d];
631 const double diff_shape2 = diffN[diff_shift + 2 * n2 + d];
632 const double diff_shape3 = diffN[diff_shift + 2 * n3 + d];
633 t_diff_ksi01(d + 1) =
634 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3);
635 t_diff_ksi12(d + 1) =
636 (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0);
637 }
638 t_diff_ksi01(0) = 0;
639 t_diff_ksi12(0) = 0;
640
641 FTensor::Tensor1<double, 3> t_cross;
642 t_cross(i) =
643 FTensor::levi_civita(i, j, k) * t_diff_ksi01(j) * t_diff_ksi12(k);
644
645 double zeta[p[0] + 1];
646 double diff_zeta[2 * (p[0] + 1)];
647 CHKERR Legendre_polynomials(p[0], ksi01, &t_diff_ksi01(0), zeta,
648 diff_zeta, 2);
649
650 double eta[p[1] + 1];
651 double diff_eta[2 * (p[1] + 1)];
652 CHKERR Legendre_polynomials(p[1], ksi12, &t_diff_ksi12(0), eta, diff_eta,
653 2);
654
655 for (int n = 0; n != nb_dofs; ++n) {
656 int ii = permute[n][0];
657 int jj = permute[n][1];
658
659 const double z = zeta[ii];
660 const double e = eta[jj];
661 const double ez = e * z;
662
663 auto t_diff_zeta = FTensor::Tensor1<double, 2>{
664 diff_zeta[ii], diff_zeta[(p[0] + 1) + ii]};
665 auto t_diff_eta = FTensor::Tensor1<double, 2>{
666 diff_eta[jj], diff_eta[(p[1] + 1) + jj]};
667
668 FTensor::Index<'J', 2> J;
669 t_n(i) = t_cross(i) * ez;
670 t_diff_n(i, J) = t_cross(i) * (t_diff_zeta(J) * e + t_diff_eta(J) * z);
671
672 ++t_n;
673 ++t_diff_n;
674 }
675 }
676 }
677 MoFEMFunctionReturnHot(0);
678}
679
680/* Reference Hex and its cannonical vertex and edge numbering
681 7 ---------10--------- 6
682 /| /|
683 / | / |
684 11 | 9 | x3
685 / 7 / | |
686 / | / 6 | x2
687 4 ----------8--------- 5 | | /
688 | | | | | /
689 | 3 ----------2---------- 2 | /
690 4 / | / o -------- x1
691 | / 5 /
692 | 3 | 1
693 | / | /
694 |/ |/
695 0 ---------0---------- 1
696
697 Hex Face Cannonical numbering
698
699 1. 0 1 2 3
700 2. 0 1 5 4
701 3. 1 2 6 5
702 4. 3 2 6 7
703 5. 0 3 7 4
704 6. 4 5 6 7
705*/
706
707MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::H1_EdgeShapeFunctions_ONHEX(
708 int *sense, int *p, double *N, double *diff_N, double *edgeN[12],
709 double *diff_edgeN[12], int nb_integration_pts) {
710 MoFEMFunctionBeginHot;
711
712 EntityType sub_type;
713 int num_sub_ent_vertices;
714 const short int *edges_conn[12];
715 for (int e = 0; e != 12; ++e)
716 edges_conn[e] =
717 CN::SubEntityVertexIndices(MBHEX, 1, e, sub_type, num_sub_ent_vertices);
718
719 const short int *face_conn[6];
720 for (int f = 0; f != 6; ++f)
721 face_conn[f] =
722 CN::SubEntityVertexIndices(MBHEX, 2, f, sub_type, num_sub_ent_vertices);
723
724 constexpr int quad_edge[12][2] = {{3, 1}, {0, 2}, {1, 3}, {2, 0},
725 {4, 5}, {4, 5}, {4, 5}, {4, 5},
726 {3, 1}, {0, 2}, {1, 3}, {2, 0}};
727
728 for (int qq = 0; qq != nb_integration_pts; ++qq) {
729
730 const int shift = 8 * qq;
731
732 double ksi[12];
733 double diff_ksi[12][3];
734 for (int e = 0; e != 12; ++e) {
735
736 ksi[e] = 0;
737 for (int d = 0; d != 3; ++d)
738 diff_ksi[e][d] = 0;
739 for (int n = 0; n != 4; ++n) {
740 const auto n1 = shift + face_conn[quad_edge[e][1]][n];
741 const auto n0 = shift + face_conn[quad_edge[e][0]][n];
742 ksi[e] += N[n1] - N[n0];
743 const auto n03 = 3 * n0;
744 const auto n13 = 3 * n1;
745 for (int d = 0; d != 3; ++d)
746 diff_ksi[e][d] += diff_N[n13 + d] - diff_N[n03 + d];
747 }
748
749 ksi[e] *= sense[e];
750 for (int d = 0; d != 3; ++d)
751 diff_ksi[e][d] *= sense[e];
752 }
753
754 double mu[12];
755 double diff_mu[12][3];
756 for (int e = 0; e != 12; ++e) {
757 const auto n0 = shift + edges_conn[e][0];
758 const auto n1 = shift + edges_conn[e][1];
759 mu[e] = N[n0] + N[n1];
760 const auto n03 = 3 * n0;
761 const auto n13 = 3 * n1;
762 for (int d = 0; d != 3; ++d) {
763 diff_mu[e][d] = diff_N[n03 + d] + diff_N[n13 + d];
764 }
765 }
766
767 for (int e = 0; e != 12; e++) {
768
769 const int nb_dofs = (p[e] - 1);
770 if (nb_dofs > 0) {
771
772 double L[p[e] + 2];
773 double diffL[3 * (p[e] + 2)];
774 CHKERR Lobatto_polynomials(p[e] + 1, ksi[e], diff_ksi[e], L, diffL, 3);
775
776 const int qd_shift = nb_dofs * qq;
777 for (int n = 0; n != nb_dofs; n++) {
778 edgeN[e][qd_shift + n] = mu[e] * L[n + 2];
779 for (int d = 0; d != 3; ++d) {
780 diff_edgeN[e][3 * (qd_shift + n) + d] =
781
782 diff_mu[e][d] * L[n + 2]
783
784 +
785
786 mu[e] * diffL[d * (p[e] + 2) + n + 2];
787 }
788 }
789 }
790 }
791 }
792 MoFEMFunctionReturnHot(0);
793}
794
795MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::H1_FaceShapeFunctions_ONHEX(
796 int *face_nodes, int *face_nodes_order, int *p, double *N, double *diffN,
797 double *faceN[6], double *diff_faceN[6], int nb_integration_pts) {
798 MoFEMFunctionBeginHot;
799
800 constexpr int opposite_face_node[6][4] = {
801
802 {3, 2, 6, 7},
803
804 {0, 3, 7, 4},
805
806 {1, 0, 4, 5},
807
808 {2, 1, 5, 6},
809
810 {4, 7, 6, 5},
811
812 {0, 1, 2, 3}
813
814 };
815
816 for (int face = 0; face != 6; face++) {
817
818 if (p[face] > 1) {
819 const int nb_dofs = (p[face] - 1) * (p[face] - 1);
820
821 const int n0 = face_nodes[4 * face + 0];
822 const int n1 = face_nodes[4 * face + 1];
823 const int n2 = face_nodes[4 * face + 2];
824 const int n3 = face_nodes[4 * face + 3];
825
826 const int o0 = opposite_face_node[face][face_nodes_order[4 * face + 0]];
827 const int o1 = opposite_face_node[face][face_nodes_order[4 * face + 1]];
828 const int o2 = opposite_face_node[face][face_nodes_order[4 * face + 2]];
829 const int o3 = opposite_face_node[face][face_nodes_order[4 * face + 3]];
830
831 int permute[nb_dofs][3];
832 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0], p[face] - 2,
833 p[face] - 2);
834
835 for (int qq = 0; qq != nb_integration_pts; qq++) {
836
837 const int shift = 8 * qq;
838
839 const double shape0 = N[shift + n0];
840 const double shape1 = N[shift + n1];
841 const double shape2 = N[shift + n2];
842 const double shape3 = N[shift + n3];
843
844 const double o_shape0 = N[shift + o0];
845 const double o_shape1 = N[shift + o1];
846 const double o_shape2 = N[shift + o2];
847 const double o_shape3 = N[shift + o3];
848
849 const double ksi01 = (shape1 + shape2 - shape0 - shape3) +
850 (o_shape1 + o_shape2 - o_shape0 - o_shape3);
851 const double ksi12 = (shape2 + shape3 - shape1 - shape0) +
852 (o_shape2 + o_shape3 - o_shape1 - o_shape0);
853 const double mu = shape1 + shape2 + shape0 + shape3;
854
855 const int diff_shift = 3 * shift;
856 double diff_ksi01[3], diff_ksi12[3], diff_mu[3];
857 for (int d = 0; d != 3; d++) {
858 const double diff_shape0 = diffN[diff_shift + 3 * n0 + d];
859 const double diff_shape1 = diffN[diff_shift + 3 * n1 + d];
860 const double diff_shape2 = diffN[diff_shift + 3 * n2 + d];
861 const double diff_shape3 = diffN[diff_shift + 3 * n3 + d];
862 const double o_diff_shape0 = diffN[diff_shift + 3 * o0 + d];
863 const double o_diff_shape1 = diffN[diff_shift + 3 * o1 + d];
864 const double o_diff_shape2 = diffN[diff_shift + 3 * o2 + d];
865 const double o_diff_shape3 = diffN[diff_shift + 3 * o3 + d];
866
867 diff_ksi01[d] =
868 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) +
869 (o_diff_shape1 + o_diff_shape2 - o_diff_shape0 - o_diff_shape3);
870 diff_ksi12[d] =
871 (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) +
872 (o_diff_shape2 + o_diff_shape3 - o_diff_shape1 - o_diff_shape0);
873 diff_mu[d] = (diff_shape1 + diff_shape2 + diff_shape0 + diff_shape3);
874 }
875
876 double L01[p[face] + 2];
877 double diffL01[3 * (p[face] + 2)];
878 CHKERR Lobatto_polynomials(p[face] + 1, ksi01, diff_ksi01, L01, diffL01,
879 3);
880 double L12[p[face] + 2];
881 double diffL12[3 * (p[face] + 2)];
882 CHKERR Lobatto_polynomials(p[face] + 1, ksi12, diff_ksi12, L12, diffL12,
883 3);
884
885 int qd_shift = nb_dofs * qq;
886 for (int n = 0; n != nb_dofs; ++n) {
887 const int s1 = permute[n][0];
888 const int s2 = permute[n][1];
889 const double vol = L01[s1 + 2] * L12[s2 + 2];
890 faceN[face][qd_shift + n] = vol * mu;
891 for (int d = 0; d != 3; ++d) {
892 diff_faceN[face][3 * (qd_shift + n) + d] =
893 (diffL01[d * (p[face] + 2) + s1 + 2] * L12[s2 + 2] +
894 diffL12[d * (p[face] + 2) + s2 + 2] * L01[s1 + 2]) *
895 mu +
896 vol * diff_mu[d];
897 }
898 }
899 }
900 }
901 }
902 MoFEMFunctionReturnHot(0);
903}
904
905MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::H1_InteriorShapeFunctions_ONHEX(
906 const int *p, double *N, double *N_diff, double *faceN, double *diff_faceN,
907 int nb_integration_pts) {
908 MoFEMFunctionBeginHot;
909
910 const int nb_bases = (p[0] - 1) * (p[1] - 1) * (p[2] - 1);
911
912 if (nb_bases > 0) {
913
914 int permute[nb_bases][3];
915 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0], p[0] - 2,
916 p[1] - 2, p[2] - 2);
917
918 // double P0[p[0] + 2];
919 // double diffL0[3 * (p[0] + 2)];
920 // double P1[p[1] + 2];
921 // double diffL1[3 * (p[1] + 2)];
922 // double P2[p[2] + 2];
923 // double diffL2[3 * (p[2] + 2)];
924
925 for (int qq = 0; qq != nb_integration_pts; ++qq) {
926
927 const int shift = 8 * qq;
928 double ksi[3] = {0, 0, 0};
929 double diff_ksi[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
930 ::DemkowiczHexAndQuad::get_ksi_hex(shift, N, N_diff, ksi, diff_ksi);
931
932 double L0[p[0] + 2];
933 double diffL0[3 * (p[0] + 2)];
934 double L1[p[1] + 2];
935 double diffL1[3 * (p[1] + 2)];
936 double L2[p[2] + 2];
937 double diffL2[3 * (p[2] + 2)];
938
939 CHKERR Lobatto_polynomials(p[0] + 1, ksi[0], diff_ksi[0], L0, diffL0, 3);
940 CHKERR Lobatto_polynomials(p[1] + 1, ksi[1], diff_ksi[1], L1, diffL1, 3);
941 CHKERR Lobatto_polynomials(p[2] + 1, ksi[2], diff_ksi[2], L2, diffL2, 3);
942
943 const int qd_shift = nb_bases * qq;
944 for (int n = 0; n != nb_bases; ++n) {
945 const int s1 = permute[n][0];
946 const int s2 = permute[n][1];
947 const int s3 = permute[n][2];
948
949 const double l0l1 = L0[s1 + 2] * L1[s2 + 2];
950 const double l0l2 = L0[s1 + 2] * L2[s3 + 2];
951 const double l1l2 = L1[s2 + 2] * L2[s3 + 2];
952
953 faceN[qd_shift + n] = l0l1 * L2[s3 + 2];
954 for (int d = 0; d != 3; ++d) {
955 diff_faceN[3 * (qd_shift + n) + d] =
956 (diffL0[d * (p[0] + 2) + s1 + 2] * l1l2 +
957 diffL1[d * (p[1] + 2) + s2 + 2] * l0l2 +
958 diffL2[d * (p[2] + 2) + s3 + 2] * l0l1);
959 }
960 }
961 }
962 }
963 MoFEMFunctionReturnHot(0);
964}
965
966MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::L2_InteriorShapeFunctions_ONHEX(
967 const int *p, double *N, double *N_diff, double *volN, double *diff_volN,
968 int nb_integration_pts) {
969 MoFEMFunctionBeginHot;
970
971 const int nb_bases = (p[0] + 1) * (p[1] + 1) * (p[2] + 1);
972 if (nb_bases > 0) {
973
974 int permute[nb_bases][3];
975 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0], p[0], p[1],
976 p[2]);
977
978 double P0[p[0] + 2];
979 double diffL0[3 * (p[0] + 2)];
980 double P1[p[1] + 2];
981 double diffL1[3 * (p[1] + 2)];
982 double P2[p[2] + 2];
983 double diffL2[3 * (p[2] + 2)];
984
985 for (int qq = 0; qq != nb_integration_pts; qq++) {
986
987 const int shift = 8 * qq;
988 double ksi[3] = {0, 0, 0};
989 double diff_ksi[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
990 ::DemkowiczHexAndQuad::get_ksi_hex(shift, N, N_diff, ksi, diff_ksi);
991
992 CHKERR Legendre_polynomials(p[0] + 1, ksi[0], diff_ksi[0], P0, diffL0, 3);
993 CHKERR Legendre_polynomials(p[1] + 1, ksi[1], diff_ksi[1], P1, diffL1, 3);
994 CHKERR Legendre_polynomials(p[2] + 1, ksi[2], diff_ksi[2], P2, diffL2, 3);
995
996 const int qd_shift = qq * nb_bases;
997 for (int n = 0; n != nb_bases; ++n) {
998 const int ii = permute[n][0];
999 const int jj = permute[n][1];
1000 const int kk = permute[n][2];
1001
1002 const double p1p2 = P1[jj] * P2[kk];
1003 const double p0p1 = P0[ii] * P1[jj];
1004 const double p0p2 = P0[ii] * P2[kk];
1005
1006 volN[qd_shift + n] = p0p1 * P2[kk];
1007 for (int d = 0; d != 3; ++d) {
1008 diff_volN[3 * (qd_shift + n) + d] =
1009 p1p2 * diffL0[d * (p[0] + 2) + ii] +
1010 p0p2 * diffL1[d * (p[1] + 2) + jj] +
1011 p0p1 * diffL2[d * (p[2] + 2) + kk];
1012 }
1013 }
1014 }
1015 }
1016 MoFEMFunctionReturnHot(0);
1017}
1018
1019MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hcurl_EdgeShapeFunctions_ONHEX(
1020 int *sense, int *p, double *N, double *diff_N, double *edgeN[12],
1021 double *diff_edgeN[12], int nb_integration_pts) {
1022
1023 FTensor::Index<'i', 3> i;
1024 FTensor::Index<'j', 3> j;
1025
1026 MoFEMFunctionBeginHot;
1027
1028 EntityType sub_type;
1029 int num_sub_ent_vertices;
1030 const short int *edges_conn[12];
1031 for (int e = 0; e != 12; ++e)
1032 edges_conn[e] =
1033 CN::SubEntityVertexIndices(MBHEX, 1, e, sub_type, num_sub_ent_vertices);
1034
1035 const short int *face_conn[6];
1036 for (int f = 0; f != 6; ++f)
1037 face_conn[f] =
1038 CN::SubEntityVertexIndices(MBHEX, 2, f, sub_type, num_sub_ent_vertices);
1039
1040 constexpr int quad_edge[12][2] = {{3, 1}, {0, 2}, {1, 3}, {2, 0},
1041 {4, 5}, {4, 5}, {4, 5}, {4, 5},
1042 {3, 1}, {0, 2}, {1, 3}, {2, 0}};
1043
1044 for (int qq = 0; qq != nb_integration_pts; qq++) {
1045
1046 double ksi[12];
1047 double diff_ksi[12 * 3];
1048
1049 const int shift = 8 * qq;
1050
1051 auto calc_ksi = [&]() {
1052 auto t_diff_ksi = getFTensor1FromPtr<3>(diff_ksi);
1053 for (int e = 0; e != 12; ++e) {
1054 if (p[e] > 0) {
1055 ksi[e] = 0;
1056 t_diff_ksi(i) = 0;
1057 for (int n = 0; n != 4; ++n) {
1058 const auto n1 = shift + face_conn[quad_edge[e][1]][n];
1059 const auto n0 = shift + face_conn[quad_edge[e][0]][n];
1060 ksi[e] += N[n1] - N[n0];
1061 const auto n03 = 3 * n0;
1062 const auto n13 = 3 * n1;
1063 for (int d = 0; d != 3; ++d)
1064 t_diff_ksi(d) += diff_N[n13 + d] - diff_N[n03 + d];
1065 }
1066
1067 ksi[e] *= sense[e];
1068 t_diff_ksi(i) *= sense[e];
1069
1070 ++t_diff_ksi;
1071 }
1072 }
1073 };
1074
1075 double mu[12];
1076 double diff_mu[12][3];
1077 auto calc_mu = [&]() {
1078 for (int e = 0; e != 12; ++e) {
1079 if (p[e] > 0) {
1080 const auto n0 = shift + edges_conn[e][0];
1081 const auto n1 = shift + edges_conn[e][1];
1082 mu[e] = N[n0] + N[n1];
1083 const auto n03 = 3 * n0;
1084 const auto n13 = 3 * n1;
1085 for (int d = 0; d != 3; ++d) {
1086 diff_mu[e][d] = diff_N[n03 + d] + diff_N[n13 + d];
1087 }
1088 }
1089 }
1090 };
1091
1092 auto calc_base = [&]() {
1093 auto t_diff_ksi = getFTensor1FromPtr<3>(diff_ksi);
1094 for (int ee = 0; ee != 12; ee++) {
1095 if (p[ee] > 0) {
1096
1097 double L[p[ee]];
1098 double diffL[3 * (p[ee])];
1099 CHKERR Legendre_polynomials(p[ee] - 1, ksi[ee], &(t_diff_ksi(0)), L,
1100 diffL, 3);
1101
1102 int qd_shift = p[ee] * qq;
1103 auto t_n = getFTensor1FromPtr<3>(&edgeN[ee][3 * qd_shift]);
1104 auto t_diff_n =
1105 getFTensor2HVecFromPtr<3, 3>(&diff_edgeN[ee][9 * qd_shift]);
1106
1107 for (int ii = 0; ii != p[ee]; ii++) {
1108
1109 const double a = mu[ee] * L[ii];
1110 auto t_diff_a = FTensor::Tensor1<const double, 3>{
1111
1112 diff_mu[ee][0] * L[ii] + diffL[0 * p[ee] + ii],
1113
1114 diff_mu[ee][1] * L[ii] + diffL[1 * p[ee] + ii],
1115
1116 diff_mu[ee][2] * L[ii] + diffL[2 * p[ee] + ii]
1117
1118 };
1119
1120 t_n(i) = (a / 2) * t_diff_ksi(i);
1121 t_diff_n(i, j) = (t_diff_a(j) / 2) * t_diff_ksi(i);
1122
1123 ++t_n;
1124 ++t_diff_n;
1125 }
1126
1127 ++t_diff_ksi;
1128 }
1129 }
1130 };
1131
1132 calc_ksi();
1133 calc_mu();
1134 calc_base();
1135 }
1136 MoFEMFunctionReturnHot(0);
1137}
1138
1139MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hcurl_FaceShapeFunctions_ONHEX(
1140 int *face_nodes, int *face_nodes_order, int *p, double *N, double *diffN,
1141 double *faceN[6][2], double *diff_faceN[6][2], int nb_integration_pts) {
1142 MoFEMFunctionBeginHot;
1143
1144 constexpr int opposite_face_node[6][4] = {
1145
1146 {3, 2, 6, 7},
1147
1148 {0, 3, 7, 4},
1149
1150 {1, 0, 4, 5},
1151
1152 {2, 1, 5, 6},
1153
1154 {4, 7, 6, 5},
1155
1156 {0, 1, 2, 3}
1157
1158 };
1159
1160 for (int face = 0; face != 6; face++) {
1161 if ((p[face] - 1) * p[face] > 0) {
1162
1163 const int n0 = face_nodes[4 * face + 0];
1164 const int n1 = face_nodes[4 * face + 1];
1165 const int n2 = face_nodes[4 * face + 2];
1166 const int n3 = face_nodes[4 * face + 3];
1167
1168 const int o0 = opposite_face_node[face][face_nodes_order[4 * face + 0]];
1169 const int o1 = opposite_face_node[face][face_nodes_order[4 * face + 1]];
1170 const int o2 = opposite_face_node[face][face_nodes_order[4 * face + 2]];
1171 const int o3 = opposite_face_node[face][face_nodes_order[4 * face + 3]];
1172
1173 int permute[(p[face] - 1) * p[face]][3];
1174 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0], p[face] - 1,
1175 p[face] - 2);
1176
1177 for (int q = 0; q != nb_integration_pts; ++q) {
1178
1179 const int shift = 8 * q;
1180
1181 const double shape0 = N[shift + n0];
1182 const double shape1 = N[shift + n1];
1183 const double shape2 = N[shift + n2];
1184 const double shape3 = N[shift + n3];
1185
1186 const double o_shape0 = N[shift + o0];
1187 const double o_shape1 = N[shift + o1];
1188 const double o_shape2 = N[shift + o2];
1189 const double o_shape3 = N[shift + o3];
1190
1191 const double ksi01 = (shape1 + shape2 - shape0 - shape3) +
1192 (o_shape1 + o_shape2 - o_shape0 - o_shape3);
1193 const double ksi12 = (shape2 + shape3 - shape1 - shape0) +
1194 (o_shape2 + o_shape3 - o_shape1 - o_shape0);
1195 const double mu = shape1 + shape2 + shape0 + shape3;
1196
1197 const int diff_shift = 3 * shift;
1198 double diff_ksi01[3], diff_ksi12[3], diff_mu[3];
1199 for (int d = 0; d != 3; ++d) {
1200 const double diff_shape0 = diffN[diff_shift + 3 * n0 + d];
1201 const double diff_shape1 = diffN[diff_shift + 3 * n1 + d];
1202 const double diff_shape2 = diffN[diff_shift + 3 * n2 + d];
1203 const double diff_shape3 = diffN[diff_shift + 3 * n3 + d];
1204 const double o_diff_shape0 = diffN[diff_shift + 3 * o0 + d];
1205 const double o_diff_shape1 = diffN[diff_shift + 3 * o1 + d];
1206 const double o_diff_shape2 = diffN[diff_shift + 3 * o2 + d];
1207 const double o_diff_shape3 = diffN[diff_shift + 3 * o3 + d];
1208
1209 diff_ksi01[d] =
1210 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) +
1211 (o_diff_shape1 + o_diff_shape2 - o_diff_shape0 - o_diff_shape3);
1212 diff_ksi12[d] =
1213 (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) +
1214 (o_diff_shape2 + o_diff_shape3 - o_diff_shape1 - o_diff_shape0);
1215 diff_mu[d] = (diff_shape1 + diff_shape2 + diff_shape0 + diff_shape3);
1216 }
1217
1218 int pq[2] = {p[face], p[face]};
1219 int qp[2] = {p[face], p[face]};
1220 double mu_ksi_eta[2] = {ksi01, ksi12};
1221 double mu_eta_ksi[2] = {ksi12, ksi01};
1222 double *diff_ksi_eta[2] = {diff_ksi01, diff_ksi12};
1223 double *diff_eta_ksi[2] = {diff_ksi12, diff_ksi01};
1224
1225 for (int family = 0; family != 2; ++family) {
1226
1227 const int pp = pq[family];
1228 const int qq = qp[family];
1229 const int nb_dofs = (pp - 1) * qq;
1230
1231 if (nb_dofs > 0) {
1232
1233 double zeta[pp + 2];
1234 double diff_zeta[3 * (pp + 2)];
1235 CHKERR Lobatto_polynomials(pp + 1, mu_ksi_eta[family],
1236 diff_ksi_eta[family], zeta, diff_zeta,
1237 3);
1238
1239 double eta[qq];
1240 double diff_eta[3 * qq];
1241 CHKERR Legendre_polynomials(qq - 1, mu_eta_ksi[family],
1242 diff_eta_ksi[family], eta, diff_eta, 3);
1243
1244 const int qd_shift = nb_dofs * q;
1245 double *t_n_ptr = &(faceN[face][family][3 * qd_shift]);
1246 double *t_diff_n_ptr = &(diff_faceN[face][family][9 * qd_shift]);
1247 auto t_n = getFTensor1FromPtr<3>(t_n_ptr);
1248 auto t_diff_n = getFTensor2HVecFromPtr<3, 3>(t_diff_n_ptr);
1249
1250 for (int n = 0; n != nb_dofs; n++) {
1251 int i = permute[n][0];
1252 int j = permute[n][1];
1253
1254 const double z = zeta[j + 2];
1255 const double e = eta[i];
1256 const double ze = z * e;
1257 const double a = ze * mu;
1258 double d_a[3];
1259
1260 for (int d = 0; d != 3; ++d)
1261 d_a[d] = (diff_zeta[d * (pp + 2) + j + 2] * e +
1262 z * diff_eta[d * qq + i]) *
1263 mu +
1264 ze * diff_mu[d];
1265
1266 for (int d = 0; d != 3; ++d) {
1267 t_n(d) = a * diff_eta_ksi[family][d];
1268 for (int m = 0; m != 3; ++m) {
1269 t_diff_n(d, m) = d_a[m] * diff_eta_ksi[family][d];
1270 }
1271 }
1272
1273 ++t_n;
1274 ++t_diff_n;
1275 }
1276 }
1277 }
1278 }
1279 }
1280 }
1281
1282 MoFEMFunctionReturnHot(0);
1283}
1284
1285MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hcurl_InteriorShapeFunctions_ONHEX(
1286 int *p, double *N, double *diffN, double *volN[], double *diff_volN[],
1287 int nb_integration_pts) {
1288 MoFEMFunctionBeginHot;
1289
1290 int pqr[3] = {p[0], p[1], p[2]};
1291 int qrp[3] = {p[1], p[2], p[0]};
1292 int rpq[3] = {p[2], p[0], p[1]};
1293
1294 const int nb_dofs_fm0 = (p[0] - 1) * p[1] * p[2];
1295 const int nb_dofs_fm1 = p[0] * (p[1] - 1) * p[2];
1296 const int nb_dofs_fm2 = p[0] * p[1] * (p[2] - 1);
1297 int permute_fm0[3 * nb_dofs_fm0];
1298 int permute_fm1[3 * nb_dofs_fm1];
1299 int permute_fm2[3 * nb_dofs_fm2];
1300
1301 std::array<int *, 3> permute = {&permute_fm0[0], &permute_fm1[0],
1302 &permute_fm2[0]};
1303
1304 for (int fam = 0; fam != 3; ++fam) {
1305 const int qqq = qrp[fam];
1306 const int rrr = rpq[fam];
1307 const int ppp = pqr[fam];
1308 CHKERR ::DemkowiczHexAndQuad::monom_ordering(permute[fam], ppp - 2, qqq - 1,
1309 rrr - 2);
1310 }
1311
1312 for (int qq = 0; qq != nb_integration_pts; ++qq) {
1313
1314 const int shift = 8 * qq;
1315 double ksi[3] = {0, 0, 0};
1316 double diff_ksi[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
1317 ::DemkowiczHexAndQuad::get_ksi_hex(shift, N, diffN, ksi, diff_ksi);
1318
1319 double ksi_eta_gma[3] = {ksi[0], ksi[1], ksi[2]};
1320 double eta_gma_ksi[3] = {ksi[1], ksi[2], ksi[0]};
1321 double gma_ksi_eta[3] = {ksi[2], ksi[0], ksi[1]};
1322
1323 double *diff_ksi_eta_gma[3] = {diff_ksi[0], diff_ksi[1], diff_ksi[2]};
1324 double *diff_eta_gma_ksi[3] = {diff_ksi[1], diff_ksi[2], diff_ksi[0]};
1325 double *diff_gma_ksi_eta[3] = {diff_ksi[2], diff_ksi[0], diff_ksi[1]};
1326
1327 int pqr[3] = {p[0], p[1], p[2]};
1328 int qrp[3] = {p[1], p[2], p[0]};
1329 int rpq[3] = {p[2], p[0], p[1]};
1330 for (int fam = 0; fam != 3; ++fam) {
1331
1332 const int qqq = qrp[fam];
1333 const int rrr = rpq[fam];
1334 const int ppp = pqr[fam];
1335
1336 const int nb_dofs = (ppp - 1) * qqq * (rrr - 1);
1337 if (nb_dofs > 0) {
1338
1339 double phi_j[ppp + 2];
1340 double diff_phi_j[3 * (ppp + 2)];
1341
1342 CHKERR Lobatto_polynomials(ppp + 1, ksi_eta_gma[fam],
1343 diff_ksi_eta_gma[fam], phi_j, diff_phi_j, 3);
1344
1345 double eta_i[qqq];
1346 double diff_eta_i[3 * qqq];
1347
1348 CHKERR Legendre_polynomials(qqq - 1, eta_gma_ksi[fam],
1349 diff_eta_gma_ksi[fam], eta_i, diff_eta_i,
1350 3);
1351
1352 double phi_k[rrr + 2];
1353 double diff_phi_k[3 * (rrr + 2)];
1354
1355 CHKERR Lobatto_polynomials(rrr + 1, gma_ksi_eta[fam],
1356 diff_gma_ksi_eta[fam], phi_k, diff_phi_k, 3);
1357
1358 int qd_shift = nb_dofs * qq;
1359 double *t_n_ptr = &volN[fam][3 * qd_shift];
1360 double *t_diff_n_ptr = &diff_volN[fam][3 * 3 * qd_shift];
1361 auto t_n = getFTensor1FromPtr<3>(t_n_ptr);
1362 auto t_diff_n = getFTensor2HVecFromPtr<3, 3>(t_diff_n_ptr);
1363
1364 int n = 0;
1365 for (; n != nb_dofs; n++) {
1366 int ii = permute[fam][3 * n + 0];
1367 int jj = permute[fam][3 * n + 1];
1368 int kk = permute[fam][3 * n + 2];
1369
1370 const double p_k = phi_k[kk + 2];
1371 const double p_j = phi_j[jj + 2];
1372 const double e = eta_i[ii];
1373 const double a = p_j * p_k * e;
1374
1375 const double d_a[] = {
1376 diff_phi_k[0 * (ppp + 2) + kk + 2] * p_j * e +
1377 p_k * diff_phi_j[0 * (rrr + 2) + jj + 2] * e +
1378 p_k * p_j * diff_eta_i[0 * qqq + ii],
1379
1380 diff_phi_k[1 * (ppp + 2) + kk + 2] * p_j * e +
1381 p_k * diff_phi_j[1 * (rrr + 2) + jj + 2] * e +
1382 p_k * p_j * diff_eta_i[1 * qqq + ii],
1383
1384 diff_phi_k[2 * (ppp + 2) + kk + 2] * p_j * e +
1385 p_k * diff_phi_j[2 * (rrr + 2) + jj + 2] * e +
1386 p_k * p_j * diff_eta_i[2 * qqq + ii]};
1387
1388 for (int d = 0; d != 3; ++d) {
1389 t_n(d) = a * diff_eta_gma_ksi[fam][d];
1390 for (int m = 0; m != 3; ++m) {
1391 t_diff_n(d, m) = d_a[m] * diff_eta_gma_ksi[fam][d];
1392 }
1393 }
1394 ++t_n;
1395 ++t_diff_n;
1396 }
1397 }
1398 }
1399 }
1400
1401 MoFEMFunctionReturnHot(0);
1402}
1403
1404MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hdiv_FaceShapeFunctions_ONHEX(
1405 int *face_nodes, int *face_nodes_order, int *p, double *N, double *diffN,
1406 double *faceN[6], double *div_faceN[6], int nb_integration_pts) {
1407 MoFEMFunctionBeginHot;
1408
1409 FTensor::Index<'i', 3> i;
1410 FTensor::Index<'j', 3> j;
1411 FTensor::Index<'k', 3> k;
1412
1413 constexpr int opposite_face_node[6][4] = {
1414
1415 {3, 2, 6, 7},
1416
1417 {0, 3, 7, 4},
1418
1419 {1, 0, 4, 5},
1420
1421 {2, 1, 5, 6},
1422
1423 {4, 7, 6, 5},
1424
1425 {0, 1, 2, 3}
1426
1427 };
1428
1429 for (int face = 0; face != 6; face++) {
1430
1431 const int nb_dofs =
1432 NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL(p[face], p[face]);
1433
1434 if (nb_dofs > 0) {
1435
1436 auto t_n = getFTensor1FromPtr<3>(faceN[face]);
1437 auto t_diff_n = getFTensor2HVecFromPtr<3, 3>(div_faceN[face]);
1438
1439 const int n0 = face_nodes[4 * face + 0];
1440 const int n1 = face_nodes[4 * face + 1];
1441 const int n2 = face_nodes[4 * face + 2];
1442 const int n3 = face_nodes[4 * face + 3];
1443
1444 const int o0 = opposite_face_node[face][face_nodes_order[4 * face + 0]];
1445 const int o1 = opposite_face_node[face][face_nodes_order[4 * face + 1]];
1446 const int o2 = opposite_face_node[face][face_nodes_order[4 * face + 2]];
1447 const int o3 = opposite_face_node[face][face_nodes_order[4 * face + 3]];
1448
1449 int permute[nb_dofs][3];
1450 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0], p[face] - 1,
1451 p[face] - 1);
1452
1453 for (int q = 0; q != nb_integration_pts; ++q) {
1454
1455 const int shift = 8 * q;
1456
1457 const double shape0 = N[shift + n0];
1458 const double shape1 = N[shift + n1];
1459 const double shape2 = N[shift + n2];
1460 const double shape3 = N[shift + n3];
1461
1462 const double o_shape0 = N[shift + o0];
1463 const double o_shape1 = N[shift + o1];
1464 const double o_shape2 = N[shift + o2];
1465 const double o_shape3 = N[shift + o3];
1466
1467 const double ksi01 = (shape1 + shape2 - shape0 - shape3) +
1468 (o_shape1 + o_shape2 - o_shape0 - o_shape3);
1469 const double ksi12 = (shape2 + shape3 - shape1 - shape0) +
1470 (o_shape2 + o_shape3 - o_shape1 - o_shape0);
1471 const double mu = shape1 + shape2 + shape0 + shape3;
1472
1473 const int diff_shift = 3 * shift;
1474 FTensor::Tensor1<double, 3> t_diff_ksi01;
1475 FTensor::Tensor1<double, 3> t_diff_ksi12;
1476 FTensor::Tensor1<double, 3> t_diff_mu;
1477
1478 for (int d = 0; d != 3; ++d) {
1479 const double diff_shape0 = diffN[diff_shift + 3 * n0 + d];
1480 const double diff_shape1 = diffN[diff_shift + 3 * n1 + d];
1481 const double diff_shape2 = diffN[diff_shift + 3 * n2 + d];
1482 const double diff_shape3 = diffN[diff_shift + 3 * n3 + d];
1483 const double o_diff_shape0 = diffN[diff_shift + 3 * o0 + d];
1484 const double o_diff_shape1 = diffN[diff_shift + 3 * o1 + d];
1485 const double o_diff_shape2 = diffN[diff_shift + 3 * o2 + d];
1486 const double o_diff_shape3 = diffN[diff_shift + 3 * o3 + d];
1487 t_diff_ksi01(d) =
1488 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) +
1489 (o_diff_shape1 + o_diff_shape2 - o_diff_shape0 - o_diff_shape3);
1490 t_diff_ksi12(d) =
1491 (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) +
1492 (o_diff_shape2 + o_diff_shape3 - o_diff_shape1 - o_diff_shape0);
1493 t_diff_mu(d) =
1494 (diff_shape1 + diff_shape2 + diff_shape0 + diff_shape3);
1495 }
1496
1497 FTensor::Tensor1<double, 3> t_cross;
1498 t_cross(i) =
1499 FTensor::levi_civita(i, j, k) * t_diff_ksi01(j) * t_diff_ksi12(k);
1500
1501 const auto p_zeta = p[face];
1502 const auto p_eta = p_zeta;
1503
1504 double zeta[p[face] + 1];
1505 double diff_zeta[3 * (p[face] + 1)];
1506 CHKERR Legendre_polynomials(p[face], ksi01, &t_diff_ksi01(0), zeta,
1507 diff_zeta, 3);
1508
1509 double eta[p_eta + 1];
1510 double diff_eta[3 * (p_eta + 1)];
1511 CHKERR Legendre_polynomials(p_eta, ksi12, &t_diff_ksi12(0), eta,
1512 diff_eta, 3);
1513
1514 for (int n = 0; n != nb_dofs; ++n) {
1515 int ii = permute[n][0];
1516 int jj = permute[n][1];
1517
1518 const double z = zeta[ii];
1519 const double e = eta[jj];
1520 const double ez = e * z;
1521
1522 auto t_diff_zeta = FTensor::Tensor1<double, 3>{
1523 diff_zeta[ii], diff_zeta[1 * (p_zeta + 1) + ii],
1524 diff_zeta[2 * (p_zeta + 1) + ii]};
1525 auto t_diff_eta = FTensor::Tensor1<double, 3>{
1526 diff_eta[jj], diff_eta[1 * (p_eta + 1) + jj],
1527 diff_eta[2 * (p_eta + 1) + jj]};
1528
1529 t_n(i) = t_cross(i) * ez * mu;
1530 t_diff_n(i, j) =
1531 t_cross(i) * ((t_diff_zeta(j) * e + z * t_diff_eta(j)) * mu +
1532 ez * t_diff_mu(j));
1533
1534 ++t_n;
1535 ++t_diff_n;
1536 }
1537 }
1538 }
1539 }
1540
1541 MoFEMFunctionReturnHot(0);
1542}
1543
1544MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hdiv_InteriorShapeFunctions_ONHEX(
1545 int *p, double *N, double *diffN, double *volN[3], double *diff_volN[3],
1546 int nb_integration_pts) {
1547 MoFEMFunctionBeginHot;
1548
1549 int pqr[3] = {p[0], p[1], p[2]};
1550 int qrp[3] = {p[1], p[2], p[0]};
1551 int rpq[3] = {p[2], p[0], p[1]};
1552
1553 int perm_fam0[3 * (pqr[0] - 1) * qrp[0] * rpq[0]];
1554 int perm_fam1[3 * (pqr[1] - 1) * qrp[1] * rpq[1]];
1555 int perm_fam2[3 * (pqr[2] - 1) * qrp[2] * rpq[2]];
1556
1557 std::array<int *, 3> permute = {perm_fam0, perm_fam1, perm_fam2};
1558 for (int fam = 0; fam != 3; ++fam) {
1559 const int ppp = pqr[fam];
1560 const int qqq = qrp[fam];
1561 const int rrr = rpq[fam];
1562 CHKERR ::DemkowiczHexAndQuad::monom_ordering(permute[fam], ppp - 2, qqq - 1,
1563 rrr - 1);
1564 }
1565
1566 FTensor::Index<'i', 3> i;
1567 FTensor::Index<'j', 3> j;
1568 FTensor::Index<'k', 3> k;
1569
1570 FTensor::Tensor1<double, 3> t_cross;
1571
1572 // = {
1573 // {1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}};
1574
1575 for (int qq = 0; qq != nb_integration_pts; ++qq) {
1576
1577 const int shift = 8 * qq;
1578 double ksi[3] = {0, 0, 0};
1579 double diff_ksi[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
1580 ::DemkowiczHexAndQuad::get_ksi_hex(shift, N, diffN, ksi, diff_ksi);
1581
1582 double ksi_eta_gma[3] = {ksi[0], ksi[1], ksi[2]};
1583 double eta_gma_ksi[3] = {ksi[1], ksi[2], ksi[0]};
1584 double gma_ksi_eta[3] = {ksi[2], ksi[0], ksi[1]};
1585
1586 double *diff_ksi_eta_gma[3] = {diff_ksi[0], diff_ksi[1], diff_ksi[2]};
1587 double *diff_eta_gma_ksi[3] = {diff_ksi[1], diff_ksi[2], diff_ksi[0]};
1588 double *diff_gma_ksi_eta[3] = {diff_ksi[2], diff_ksi[0], diff_ksi[1]};
1589
1590 for (int fam = 0; fam != 3; ++fam) {
1591
1592 const int ppp = pqr[fam];
1593 const int qqq = qrp[fam];
1594 const int rrr = rpq[fam];
1595
1596 const int nb_dofs = (ppp - 1) * qqq * rrr;
1597 if (nb_dofs > 0) {
1598
1599 FTensor::Tensor1<double, 3> t_e1{diff_eta_gma_ksi[fam][0],
1600 diff_eta_gma_ksi[fam][1],
1601 diff_eta_gma_ksi[fam][2]};
1602 FTensor::Tensor1<double, 3> t_e2{diff_gma_ksi_eta[fam][0],
1603 diff_gma_ksi_eta[fam][1],
1604 diff_gma_ksi_eta[fam][2]};
1605
1606 t_cross(i) = FTensor::levi_civita(i, j, k) * t_e1(j) * t_e2(k);
1607
1608 double eta_i[ppp + 2];
1609 double diff_eta_i[3 * (ppp + 2)];
1610
1611 CHKERR Lobatto_polynomials(ppp + 1, ksi_eta_gma[fam],
1612 diff_ksi_eta_gma[fam], eta_i, diff_eta_i, 3);
1613
1614 double phi_j[qqq];
1615 double diff_phi_j[3 * qqq];
1616
1617 CHKERR Legendre_polynomials(qqq - 1, eta_gma_ksi[fam],
1618 diff_eta_gma_ksi[fam], phi_j, diff_phi_j,
1619 3);
1620
1621 double phi_k[rrr];
1622 double diff_phi_k[3 * rrr];
1623
1624 CHKERR Legendre_polynomials(rrr - 1, gma_ksi_eta[fam],
1625 diff_gma_ksi_eta[fam], phi_k, diff_phi_k,
1626 3);
1627
1628 int qd_shift = nb_dofs * qq;
1629 double *t_n_ptr = &volN[fam][3 * qd_shift];
1630 double *t_diff_n_ptr = &diff_volN[fam][9 * qd_shift];
1631 auto t_n = getFTensor1FromPtr<3>(t_n_ptr);
1632 auto t_diff_n = getFTensor2HVecFromPtr<3, 3>(t_diff_n_ptr);
1633
1634 for (int n = 0; n != nb_dofs; n++) {
1635 int ii = permute[fam][3 * n + 0];
1636 int jj = permute[fam][3 * n + 1];
1637 int kk = permute[fam][3 * n + 2];
1638
1639 const double e_i = eta_i[ii + 2];
1640 const double p_j = phi_j[jj];
1641 const double p_k = phi_k[kk];
1642
1643 const double p_jk = p_j * p_k;
1644 const double ep_ij = e_i * p_j;
1645 const double ep_ik = e_i * p_k;
1646
1647 const double a = e_i * p_jk;
1648
1649 FTensor::Tensor1<double, 3> t_d_a;
1650 for (int d = 0; d != 3; ++d)
1651 t_d_a(d) = diff_eta_i[d * (ppp + 2) + ii + 2] * p_jk +
1652 diff_phi_j[d * qqq + jj] * ep_ik +
1653 diff_phi_k[d * rrr + kk] * ep_ij;
1654
1655 t_n(i) = a * t_cross(i);
1656 t_diff_n(i, j) = t_cross(i) * t_d_a(j);
1657
1658 ++t_n;
1659 ++t_diff_n;
1660 }
1661 }
1662 }
1663 }
1664
1665 MoFEMFunctionReturnHot(0);
1666}
J
static Index< 'J', 3 > J
Definition: BasicFeTools.hpp:93
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65
eta
double eta
Definition: free_surface.cpp:170
Legendre_polynomials
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
Definition: base_functions.c:15
Lobatto_polynomials
PetscErrorCode Lobatto_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Lobatto base functions .
Definition: base_functions.c:199
NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL
#define NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL(P, Q)
Definition: h1_hdiv_hcurl_l2.h:143
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
DemkowiczHexAndQuad::monom_ordering
MoFEMErrorCode monom_ordering(int *perm, int p, int q, int r=0)
Definition: EdgeQuadHexPolynomials.cpp:11
DemkowiczHexAndQuad::get_ksi_hex
static void get_ksi_hex(int shift, double *N, double *N_diff, double ksi[3], double diff_ksi[3][3])
Definition: EdgeQuadHexPolynomials.cpp:103
FTensor::d
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
FTensor::permute
U permute(const Tensor2_Expr< A, T, Dim0_0, Dim0_1, i0, j0 > &, const Tensor2_Expr< B, U, Dim1_0, Dim1_1, i1, j1 > &rhs, const int N0, const int N1)
Definition: permute.hpp:11
MoFEM::DemkowiczHexAndQuad::Hcurl_EdgeShapeFunctions_ONHEX
MoFEMErrorCode Hcurl_EdgeShapeFunctions_ONHEX(int *sense, int *p, double *N, double *N_diff, double *edgeN[12], double *diff_edgeN[12], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:1019
MoFEM::DemkowiczHexAndQuad::L2_FaceShapeFunctions_ONQUAD
MoFEMErrorCode L2_FaceShapeFunctions_ONQUAD(int *p, double *N, double *diffN, double *face_buble, double *diff_face_bubble, int nb_integration_pts)
L2 Face base functions on Quad.
Definition: EdgeQuadHexPolynomials.cpp:335
MoFEM::DemkowiczHexAndQuad::H1_EdgeShapeFunctions_ONHEX
MoFEMErrorCode H1_EdgeShapeFunctions_ONHEX(int *sense, int *p, double *N, double *N_diff, double *edgeN[12], double *diff_edgeN[12], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:707
MoFEM::DemkowiczHexAndQuad::H1_FaceShapeFunctions_ONQUAD
MoFEMErrorCode H1_FaceShapeFunctions_ONQUAD(int *face_nodes, int *p, double *N, double *diffN, double *faceN, double *diff_faceN, int nb_integration_pts)
H1 Face bubble functions on Quad.
Definition: EdgeQuadHexPolynomials.cpp:275
MoFEM::DemkowiczHexAndQuad::L2_InteriorShapeFunctions_ONHEX
MoFEMErrorCode L2_InteriorShapeFunctions_ONHEX(const int *p, double *N, double *N_diff, double *volN, double *diff_volN, int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:966
MoFEM::DemkowiczHexAndQuad::Hcurl_EdgeShapeFunctions_ONQUAD
MoFEMErrorCode Hcurl_EdgeShapeFunctions_ONQUAD(int *sense, int *p, double *N, double *diffN, double *edgeN[4], double *curl_edgeN[4], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:393
MoFEM::DemkowiczHexAndQuad::H1_FaceShapeFunctions_ONHEX
MoFEMErrorCode H1_FaceShapeFunctions_ONHEX(int *face_nodes, int *face_nodes_order, int *p, double *N, double *N_diff, double *faceN[6], double *diff_faceN[6], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:795
MoFEM::DemkowiczHexAndQuad::Hcurl_FaceShapeFunctions_ONQUAD
MoFEMErrorCode Hcurl_FaceShapeFunctions_ONQUAD(int *face_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:484
MoFEM::DemkowiczHexAndQuad::Hdiv_FaceShapeFunctions_ONHEX
MoFEMErrorCode Hdiv_FaceShapeFunctions_ONHEX(int *face_nodes, int *face_nodes_order, int *p, double *N, double *diffN, double *faceN[6], double *div_faceN[6], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:1404
MoFEM::DemkowiczHexAndQuad::H1_EdgeShapeFunctions_ONQUAD
MoFEMErrorCode H1_EdgeShapeFunctions_ONQUAD(int *sense, int *p, double *N, double *diffN, double *edgeN[4], double *edgeDiffN[4], int nb_integration_pts)
H1 Edge base functions on Quad.
Definition: EdgeQuadHexPolynomials.cpp:202
MoFEM::DemkowiczHexAndQuad::Hcurl_FaceShapeFunctions_ONHEX
MoFEMErrorCode Hcurl_FaceShapeFunctions_ONHEX(int *face_nodes, int *face_nodes_order, int *p, double *N, double *N_diff, double *faceN[6][2], double *diff_faceN[6][2], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:1139
MoFEM::DemkowiczHexAndQuad::L2_ShapeFunctions_ONSEGMENT
MoFEMErrorCode L2_ShapeFunctions_ONSEGMENT(int p, double *N, double *diffN, double *funN, double *funDiffN, int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:166
MoFEM::DemkowiczHexAndQuad::H1_BubbleShapeFunctions_ONSEGMENT
MoFEMErrorCode H1_BubbleShapeFunctions_ONSEGMENT(int p, double *N, double *diffN, double *bubbleN, double *diff_bubbleN, int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:141
MoFEM::DemkowiczHexAndQuad::H1_InteriorShapeFunctions_ONHEX
MoFEMErrorCode H1_InteriorShapeFunctions_ONHEX(const int *p, double *N, double *N_diff, double *faceN, double *diff_faceN, int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:905
MoFEM::DemkowiczHexAndQuad::Hdiv_InteriorShapeFunctions_ONHEX
MoFEMErrorCode Hdiv_InteriorShapeFunctions_ONHEX(int *p, double *N, double *N_diff, double *bubleN[3], double *div_bubleN[3], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:1544
MoFEM::DemkowiczHexAndQuad::Hdiv_FaceShapeFunctions_ONQUAD
MoFEMErrorCode Hdiv_FaceShapeFunctions_ONQUAD(int *face_nodes, int *p, double *N, double *diffN, double *faceN, double *diff_faceN, int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:590
MoFEM::DemkowiczHexAndQuad::Hcurl_InteriorShapeFunctions_ONHEX
MoFEMErrorCode Hcurl_InteriorShapeFunctions_ONHEX(int *p, double *N, double *N_diff, double *volN[3], double *diff_volN[3], int nb_integration_pts)
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::getFTensor2HVecFromPtr< 3, 3 >
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:1047
MoFEM::getFTensor2HVecFromPtr< 3, 2 >
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2HVecFromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:1036
MoFEM::edges_conn
static constexpr int edges_conn[]
Definition: EntityRefine.cpp:18
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:177
N
const int N
Definition: speed_test.cpp:3
Generated by Doxygen 1.9.5 and hosted at University of Glasgow