16 for (
int m = 0;
m != std::max(std::max(
p, q), r) + 1; ++
m) {
18 const int i = std::min(
m,
p);
19 const int j = std::min(
m, q);
20 const int k = std::min(
m, r);
23 for (
int jj = 0; jj !=
j; ++jj) {
24 for (
int kk = 0; kk !=
k; ++kk) {
32 for (
int ii = 0; ii !=
i; ++ii) {
33 for (
int kk = 0; kk !=
k; ++kk) {
41 for (
int ii = 0; ii !=
i; ++ii) {
42 for (
int jj = 0; jj !=
j; ++jj) {
50 for (
int ii = 0; ii !=
i; ++ii) {
57 for (
int jj = 0; jj !=
j; ++jj) {
64 for (
int kk = 0; kk !=
k; ++kk) {
76 for (
int m = 0;
m != std::max(
p, q) + 1; ++
m) {
78 const int i = std::min(
m,
p);
79 const int j = std::min(
m, q);
82 for (
int ii = 0; ii !=
i; ++ii) {
89 for (
int jj = 0; jj !=
j; ++jj) {
104 double ksi[3],
double diff_ksi[3][3]) {
106 constexpr std::array<size_t, 4> ksi_nodes[2][3] = {
108 {{1, 2, 6, 5}, {3, 2, 6, 7}, {4, 5, 6, 7}},
110 {{0, 3, 7, 4}, {0, 1, 5, 4}, {0, 1, 2, 3}}
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];
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]) {
125 diff_ksi[
i][d] += N_diff[3 * (shift +
n) + d];
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];
142 int p,
double *
N,
double *diffN,
double *bubbleN,
double *diff_bubbleN,
143 int nb_integration_pts) {
145 const int nb_dofs = (
p - 1);
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++) {
152 double mu =
N[shift + n1] -
N[shift + n0];
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];
167 int p,
double *
N,
double *diffN,
double *funN,
double *funDiffN,
168 int nb_integration_pts) {
170 const int nb_dofs =
p + 1;
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++) {
177 double mu =
N[shift + n1] -
N[shift + n0];
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];
203 int *sense,
int *
p,
double *
N,
double *diffN,
double *edgeN[4],
204 double *diff_edgeN[4],
int nb_integration_pts) {
207 constexpr int n0 = 0;
208 constexpr int n1 = 1;
209 constexpr int n2 = 2;
210 constexpr int n3 = 3;
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];
223 const double mu[] = {ksi01, ksi12, ksi23, ksi30};
224 const double mu_const[] = {shape0 + shape1, shape1 + shape2,
225 shape2 + shape3, shape3 + shape0};
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;
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;
254 for (
int e = 0; e != 4; e++) {
255 const int nb_dofs =
p[e] - 1;
258 double diffL[2 * (
p[e] + 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];
276 int *faces_nodes,
int *
p,
double *
N,
double *diffN,
double *faceN,
277 double *diff_faceN,
int nb_integration_pts) {
279 if (
p[0] > 1 &&
p[1] > 1) {
280 const int nb_dofs = (
p[0] - 1) * (
p[1] - 1);
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];
288 int permute[(
p[0] - 1) * (
p[1] - 1)][3];
289 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0],
p[0] - 2,
291 for (
int q = 0; q != nb_integration_pts; q++) {
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);
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);
312 double L01[
p[0] + 2];
313 double diffL01[2 * (
p[0] + 2)];
315 double L12[
p[1] + 2];
316 double diffL12[2 * (
p[1] + 2)];
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];
336 int *
p,
double *
N,
double *diffN,
double *faceN,
double *diff_faceN,
337 int nb_integration_pts) {
339 const int nb_dofs = (
p[0] + 1) * (
p[1] + 1);
342 int permute[nb_dofs][3];
343 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0],
p[0],
p[1]);
345 constexpr int n0 = 0;
346 constexpr int n1 = 1;
347 constexpr int n2 = 2;
348 constexpr int n3 = 3;
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);
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);
370 double L01[
p[0] + 2];
371 double diffL01[2 * (
p[0] + 2)];
373 double L12[
p[1] + 2];
374 double diffL12[2 * (
p[1] + 2)];
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];
394 int *sense,
int *
p,
double *
N,
double *diffN,
double *edgeN[4],
395 double *diff_edgeN[4],
int nb_integration_pts) {
398 constexpr int n0 = 0;
399 constexpr int n1 = 1;
400 constexpr int n2 = 2;
401 constexpr int n3 = 3;
403 for (
int q = 0; q != nb_integration_pts; q++) {
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];
415 const double mu[] = {ksi01, ksi12, ksi23, ksi30};
416 const double mu_const[] = {shape0 + shape1, shape1 + shape2,
417 shape2 + shape3, shape3 + shape0};
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;
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;
446 for (
int e = 0; e != 4; e++) {
450 double diffL[2 *
p[e]];
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);
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]};
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];
485 int *face_nodes,
int *
p,
double *
N,
double *diffN,
double *faceN[],
486 double *diff_faceN[],
int nb_integration_pts) {
489 const int pq[2] = {
p[0],
p[1]};
490 const int qp[2] = {
p[1],
p[0]};
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];
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);
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];
509 for (
int q = 0; q != nb_integration_pts; q++) {
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);
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);
530 double mu_ksi_eta[2] = {ksi01, ksi12};
531 double mu_eta_ksi[2] = {ksi12, ksi01};
533 double *diff_ksi_eta[2] = {diff_ksi01, diff_ksi12};
534 double *diff_eta_ksi[2] = {diff_ksi12, diff_ksi01};
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;
544 double diff_zeta[2 * (pp + 2)];
546 diff_ksi_eta[family],
zeta, diff_zeta, 2);
549 double diff_eta[2 * qq];
551 diff_eta_ksi[family],
eta, diff_eta, 2);
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);
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];
563 const double z =
zeta[
j + 2];
564 const double e =
eta[
i];
565 const double a = z * e;
566 const double d_a[] = {
568 diff_zeta[0 * (pp + 2) +
j + 2] * e + z * diff_eta[0 * qq +
i],
570 diff_zeta[1 * (pp + 2) +
j + 2] * e + z * diff_eta[1 * qq +
i]};
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];
591 int *face_nodes,
int *
p,
double *
N,
double *diffN,
double *faceN,
592 double *diff_faceN,
int nb_integration_pts) {
595 const int nb_dofs = (
p[0] *
p[1]);
599 int permute[nb_dofs][3];
600 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0],
p[0] - 1,
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];
612 auto t_n = getFTensor1FromPtr<3>(faceN);
615 for (
int q = 0; q != nb_integration_pts; q++) {
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);
625 const int diff_shift = 2 * shift;
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);
645 double zeta[
p[0] + 1];
646 double diff_zeta[2 * (
p[0] + 1)];
650 double eta[
p[1] + 1];
651 double diff_eta[2 * (
p[1] + 1)];
655 for (
int n = 0;
n != nb_dofs; ++
n) {
656 int ii = permute[
n][0];
657 int jj = permute[
n][1];
659 const double z =
zeta[ii];
660 const double e =
eta[jj];
661 const double ez = e * z;
664 diff_zeta[ii], diff_zeta[(
p[0] + 1) + ii]};
666 diff_eta[jj], diff_eta[(
p[1] + 1) + jj]};
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);
708 int *sense,
int *
p,
double *
N,
double *diff_N,
double *edgeN[12],
709 double *diff_edgeN[12],
int nb_integration_pts) {
713 int num_sub_ent_vertices;
715 for (
int e = 0; e != 12; ++e)
717 CN::SubEntityVertexIndices(MBHEX, 1, e, sub_type, num_sub_ent_vertices);
719 const short int *face_conn[6];
720 for (
int f = 0; f != 6; ++f)
722 CN::SubEntityVertexIndices(MBHEX, 2, f, sub_type, num_sub_ent_vertices);
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}};
728 for (
int qq = 0; qq != nb_integration_pts; ++qq) {
730 const int shift = 8 * qq;
733 double diff_ksi[12][3];
734 for (
int e = 0; e != 12; ++e) {
737 for (
int d = 0; d != 3; ++d)
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];
750 for (
int d = 0; d != 3; ++d)
751 diff_ksi[e][d] *= sense[e];
755 double diff_mu[12][3];
756 for (
int e = 0; e != 12; ++e) {
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];
767 for (
int e = 0; e != 12; e++) {
769 const int nb_dofs = (
p[e] - 1);
773 double diffL[3 * (
p[e] + 2)];
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] =
782 diff_mu[e][d] *
L[
n + 2]
786 mu[e] * diffL[d * (
p[e] + 2) +
n + 2];
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) {
800 constexpr int opposite_face_node[6][4] = {
816 for (
int face = 0; face != 6; face++) {
819 const int nb_dofs = (
p[face] - 1) * (
p[face] - 1);
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];
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]];
831 int permute[nb_dofs][3];
832 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0],
p[face] - 2,
835 for (
int qq = 0; qq != nb_integration_pts; qq++) {
837 const int shift = 8 * qq;
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];
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];
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;
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];
868 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) +
869 (o_diff_shape1 + o_diff_shape2 - o_diff_shape0 - o_diff_shape3);
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);
876 double L01[
p[face] + 2];
877 double diffL01[3 * (
p[face] + 2)];
880 double L12[
p[face] + 2];
881 double diffL12[3 * (
p[face] + 2)];
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]) *
906 const int *
p,
double *
N,
double *N_diff,
double *faceN,
double *diff_faceN,
907 int nb_integration_pts) {
910 const int nb_bases = (
p[0] - 1) * (
p[1] - 1) * (
p[2] - 1);
914 int permute[nb_bases][3];
915 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0],
p[0] - 2,
925 for (
int qq = 0; qq != nb_integration_pts; ++qq) {
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}};
933 double diffL0[3 * (
p[0] + 2)];
935 double diffL1[3 * (
p[1] + 2)];
937 double diffL2[3 * (
p[2] + 2)];
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];
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];
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);
967 const int *
p,
double *
N,
double *N_diff,
double *volN,
double *diff_volN,
968 int nb_integration_pts) {
971 const int nb_bases = (
p[0] + 1) * (
p[1] + 1) * (
p[2] + 1);
974 int permute[nb_bases][3];
975 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0],
p[0],
p[1],
979 double diffL0[3 * (
p[0] + 2)];
981 double diffL1[3 * (
p[1] + 2)];
983 double diffL2[3 * (
p[2] + 2)];
985 for (
int qq = 0; qq != nb_integration_pts; qq++) {
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}};
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];
1002 const double p1p2 = P1[jj] * P2[kk];
1003 const double p0p1 = P0[ii] * P1[jj];
1004 const double p0p2 = P0[ii] * P2[kk];
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];
1020 int *sense,
int *
p,
double *
N,
double *diff_N,
double *edgeN[12],
1021 double *diff_edgeN[12],
int nb_integration_pts) {
1029 int num_sub_ent_vertices;
1031 for (
int e = 0; e != 12; ++e)
1033 CN::SubEntityVertexIndices(MBHEX, 1, e, sub_type, num_sub_ent_vertices);
1035 const short int *face_conn[6];
1036 for (
int f = 0; f != 6; ++f)
1038 CN::SubEntityVertexIndices(MBHEX, 2, f, sub_type, num_sub_ent_vertices);
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}};
1044 for (
int qq = 0; qq != nb_integration_pts; qq++) {
1047 double diff_ksi[12 * 3];
1049 const int shift = 8 * qq;
1051 auto calc_ksi = [&]() {
1052 auto t_diff_ksi = getFTensor1FromPtr<3>(diff_ksi);
1053 for (
int e = 0; e != 12; ++e) {
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];
1068 t_diff_ksi(
i) *= sense[e];
1076 double diff_mu[12][3];
1077 auto calc_mu = [&]() {
1078 for (
int e = 0; e != 12; ++e) {
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];
1092 auto calc_base = [&]() {
1093 auto t_diff_ksi = getFTensor1FromPtr<3>(diff_ksi);
1094 for (
int ee = 0; ee != 12; ee++) {
1098 double diffL[3 * (
p[ee])];
1102 int qd_shift =
p[ee] * qq;
1103 auto t_n = getFTensor1FromPtr<3>(&edgeN[ee][3 * qd_shift]);
1107 for (
int ii = 0; ii !=
p[ee]; ii++) {
1109 const double a =
mu[ee] *
L[ii];
1112 diff_mu[ee][0] *
L[ii] + diffL[0 *
p[ee] + ii],
1114 diff_mu[ee][1] *
L[ii] + diffL[1 *
p[ee] + ii],
1116 diff_mu[ee][2] *
L[ii] + diffL[2 *
p[ee] + ii]
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);
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) {
1144 constexpr int opposite_face_node[6][4] = {
1160 for (
int face = 0; face != 6; face++) {
1161 if ((
p[face] - 1) *
p[face] > 0) {
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];
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]];
1173 int permute[(
p[face] - 1) *
p[face]][3];
1174 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0],
p[face] - 1,
1177 for (
int q = 0; q != nb_integration_pts; ++q) {
1179 const int shift = 8 * q;
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];
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];
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;
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];
1210 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) +
1211 (o_diff_shape1 + o_diff_shape2 - o_diff_shape0 - o_diff_shape3);
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);
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};
1225 for (
int family = 0; family != 2; ++family) {
1227 const int pp = pq[family];
1228 const int qq = qp[family];
1229 const int nb_dofs = (pp - 1) * qq;
1233 double zeta[pp + 2];
1234 double diff_zeta[3 * (pp + 2)];
1236 diff_ksi_eta[family],
zeta, diff_zeta,
1240 double diff_eta[3 * qq];
1242 diff_eta_ksi[family],
eta, diff_eta, 3);
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);
1250 for (
int n = 0;
n != nb_dofs;
n++) {
1251 int i = permute[
n][0];
1252 int j = permute[
n][1];
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;
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]) *
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];
1286 int *
p,
double *
N,
double *diffN,
double *volN[],
double *diff_volN[],
1287 int nb_integration_pts) {
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]};
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];
1301 std::array<int *, 3>
permute = {&permute_fm0[0], &permute_fm1[0],
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,
1312 for (
int qq = 0; qq != nb_integration_pts; ++qq) {
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}};
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]};
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]};
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) {
1332 const int qqq = qrp[fam];
1333 const int rrr = rpq[fam];
1334 const int ppp = pqr[fam];
1336 const int nb_dofs = (ppp - 1) * qqq * (rrr - 1);
1339 double phi_j[ppp + 2];
1340 double diff_phi_j[3 * (ppp + 2)];
1343 diff_ksi_eta_gma[fam], phi_j, diff_phi_j, 3);
1346 double diff_eta_i[3 * qqq];
1349 diff_eta_gma_ksi[fam], eta_i, diff_eta_i,
1352 double phi_k[rrr + 2];
1353 double diff_phi_k[3 * (rrr + 2)];
1356 diff_gma_ksi_eta[fam], phi_k, diff_phi_k, 3);
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);
1365 for (;
n != nb_dofs;
n++) {
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;
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],
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],
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]};
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];
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) {
1413 constexpr int opposite_face_node[6][4] = {
1429 for (
int face = 0; face != 6; face++) {
1431 const int nb_dofs = (
p[face] *
p[face]);
1435 auto t_n = getFTensor1FromPtr<3>(faceN[face]);
1438 const int n0 = face_nodes[4 * face + 0];
1439 const int n1 = face_nodes[4 * face + 1];
1440 const int n2 = face_nodes[4 * face + 2];
1441 const int n3 = face_nodes[4 * face + 3];
1443 const int o0 = opposite_face_node[face][face_nodes_order[4 * face + 0]];
1444 const int o1 = opposite_face_node[face][face_nodes_order[4 * face + 1]];
1445 const int o2 = opposite_face_node[face][face_nodes_order[4 * face + 2]];
1446 const int o3 = opposite_face_node[face][face_nodes_order[4 * face + 3]];
1448 int permute[nb_dofs][3];
1449 CHKERR ::DemkowiczHexAndQuad::monom_ordering(&permute[0][0],
p[face] - 1,
1452 for (
int q = 0; q != nb_integration_pts; ++q) {
1454 const int shift = 8 * q;
1456 const double shape0 =
N[shift + n0];
1457 const double shape1 =
N[shift + n1];
1458 const double shape2 =
N[shift + n2];
1459 const double shape3 =
N[shift + n3];
1461 const double o_shape0 =
N[shift + o0];
1462 const double o_shape1 =
N[shift + o1];
1463 const double o_shape2 =
N[shift + o2];
1464 const double o_shape3 =
N[shift + o3];
1466 const double ksi01 = (shape1 + shape2 - shape0 - shape3) +
1467 (o_shape1 + o_shape2 - o_shape0 - o_shape3);
1468 const double ksi12 = (shape2 + shape3 - shape1 - shape0) +
1469 (o_shape2 + o_shape3 - o_shape1 - o_shape0);
1470 const double mu = shape1 + shape2 + shape0 + shape3;
1472 const int diff_shift = 3 * shift;
1477 for (
int d = 0; d != 3; ++d) {
1478 const double diff_shape0 = diffN[diff_shift + 3 * n0 + d];
1479 const double diff_shape1 = diffN[diff_shift + 3 * n1 + d];
1480 const double diff_shape2 = diffN[diff_shift + 3 * n2 + d];
1481 const double diff_shape3 = diffN[diff_shift + 3 * n3 + d];
1482 const double o_diff_shape0 = diffN[diff_shift + 3 * o0 + d];
1483 const double o_diff_shape1 = diffN[diff_shift + 3 * o1 + d];
1484 const double o_diff_shape2 = diffN[diff_shift + 3 * o2 + d];
1485 const double o_diff_shape3 = diffN[diff_shift + 3 * o3 + d];
1487 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) +
1488 (o_diff_shape1 + o_diff_shape2 - o_diff_shape0 - o_diff_shape3);
1490 (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) +
1491 (o_diff_shape2 + o_diff_shape3 - o_diff_shape1 - o_diff_shape0);
1493 (diff_shape1 + diff_shape2 + diff_shape0 + diff_shape3);
1500 double zeta[
p[0] + 1];
1501 double diff_zeta[3 * (
p[0] + 1)];
1505 double eta[
p[1] + 1];
1506 double diff_eta[3 * (
p[1] + 1)];
1510 for (
int n = 0;
n != nb_dofs; ++
n) {
1511 int ii = permute[
n][0];
1512 int jj = permute[
n][1];
1514 const double z =
zeta[ii];
1515 const double e =
eta[jj];
1516 const double ez = e * z;
1519 diff_zeta[ii], diff_zeta[1 * (
p[0] + 1) + ii],
1520 diff_zeta[2 * (
p[0] + 1) + ii]};
1522 diff_eta[jj], diff_eta[1 * (
p[1] + 1) + jj],
1523 diff_eta[2 * (
p[1] + 1) + jj]};
1525 t_n(
i) = t_cross(
i) * ez *
mu;
1527 t_cross(
i) * ((t_diff_zeta(
j) * e + z * t_diff_eta(
j)) *
mu +
1541 int *
p,
double *
N,
double *diffN,
double *volN[3],
double *diff_volN[3],
1542 int nb_integration_pts) {
1545 int pqr[3] = {
p[0],
p[1],
p[2]};
1546 int qrp[3] = {
p[1],
p[2],
p[0]};
1547 int rpq[3] = {
p[2],
p[0],
p[1]};
1549 int perm_fam0[3 * (pqr[0] - 1) * qrp[0] * rpq[0]];
1550 int perm_fam1[3 * (pqr[1] - 1) * qrp[1] * rpq[1]];
1551 int perm_fam2[3 * (pqr[2] - 1) * qrp[2] * rpq[2]];
1553 std::array<int *, 3> permute = {perm_fam0, perm_fam1, perm_fam2};
1554 for (
int fam = 0; fam != 3; ++fam) {
1555 const int ppp = pqr[fam];
1556 const int qqq = qrp[fam];
1557 const int rrr = rpq[fam];
1558 CHKERR ::DemkowiczHexAndQuad::monom_ordering(permute[fam], ppp - 2, qqq - 1,
1571 for (
int qq = 0; qq != nb_integration_pts; ++qq) {
1573 const int shift = 8 * qq;
1574 double ksi[3] = {0, 0, 0};
1575 double diff_ksi[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
1578 double ksi_eta_gma[3] = {ksi[0], ksi[1], ksi[2]};
1579 double eta_gma_ksi[3] = {ksi[1], ksi[2], ksi[0]};
1580 double gma_ksi_eta[3] = {ksi[2], ksi[0], ksi[1]};
1582 double *diff_ksi_eta_gma[3] = {diff_ksi[0], diff_ksi[1], diff_ksi[2]};
1583 double *diff_eta_gma_ksi[3] = {diff_ksi[1], diff_ksi[2], diff_ksi[0]};
1584 double *diff_gma_ksi_eta[3] = {diff_ksi[2], diff_ksi[0], diff_ksi[1]};
1586 for (
int fam = 0; fam != 3; ++fam) {
1588 const int ppp = pqr[fam];
1589 const int qqq = qrp[fam];
1590 const int rrr = rpq[fam];
1592 const int nb_dofs = (ppp - 1) * qqq * rrr;
1596 diff_eta_gma_ksi[fam][1],
1597 diff_eta_gma_ksi[fam][2]};
1599 diff_gma_ksi_eta[fam][1],
1600 diff_gma_ksi_eta[fam][2]};
1604 double eta_i[ppp + 2];
1605 double diff_eta_i[3 * (ppp + 2)];
1608 diff_ksi_eta_gma[fam], eta_i, diff_eta_i, 3);
1611 double diff_phi_j[3 * qqq];
1614 diff_eta_gma_ksi[fam], phi_j, diff_phi_j,
1618 double diff_phi_k[3 * rrr];
1621 diff_gma_ksi_eta[fam], phi_k, diff_phi_k,
1624 int qd_shift = nb_dofs * qq;
1625 double *t_n_ptr = &volN[fam][3 * qd_shift];
1626 double *t_diff_n_ptr = &diff_volN[fam][9 * qd_shift];
1627 auto t_n = getFTensor1FromPtr<3>(t_n_ptr);
1630 for (
int n = 0;
n != nb_dofs;
n++) {
1631 int ii = permute[fam][3 *
n + 0];
1632 int jj = permute[fam][3 *
n + 1];
1633 int kk = permute[fam][3 *
n + 2];
1635 const double e_i = eta_i[ii + 2];
1636 const double p_j = phi_j[jj];
1637 const double p_k = phi_k[kk];
1639 const double p_jk = p_j * p_k;
1640 const double ep_ij = e_i * p_j;
1641 const double ep_ik = e_i * p_k;
1643 const double a = e_i * p_jk;
1646 for (
int d = 0; d != 3; ++d)
1647 t_d_a(d) = diff_eta_i[d * (ppp + 2) + ii + 2] * p_jk +
1648 diff_phi_j[d * qqq + jj] * ep_ik +
1649 diff_phi_k[d * rrr + kk] * ep_ij;
1651 t_n(
i) =
a * t_cross(
i);
1652 t_diff_n(
i,
j) = t_cross(
i) * t_d_a(
j);
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
PetscErrorCode Lobatto_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Lobatto base functions .
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
MoFEMErrorCode monom_ordering(int *perm, int p, int q, int r=0)
static void get_ksi_hex(int shift, double *N, double *N_diff, double ksi[3], double diff_ksi[3][3])
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)
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
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)
MoFEMErrorCode Hcurl_EdgeShapeFunctions_ONHEX(int *sense, int *p, double *N, double *N_diff, double *edgeN[12], double *diff_edgeN[12], int nb_integration_pts)
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.
MoFEMErrorCode H1_EdgeShapeFunctions_ONHEX(int *sense, int *p, double *N, double *N_diff, double *edgeN[12], double *diff_edgeN[12], int nb_integration_pts)
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.
MoFEMErrorCode L2_InteriorShapeFunctions_ONHEX(const int *p, double *N, double *N_diff, double *volN, double *diff_volN, int nb_integration_pts)
MoFEMErrorCode Hcurl_EdgeShapeFunctions_ONQUAD(int *sense, int *p, double *N, double *diffN, double *edgeN[4], double *curl_edgeN[4], int nb_integration_pts)
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)
MoFEMErrorCode Hcurl_FaceShapeFunctions_ONQUAD(int *face_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int nb_integration_pts)
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)
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.
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)
MoFEMErrorCode L2_ShapeFunctions_ONSEGMENT(int p, double *N, double *diffN, double *funN, double *funDiffN, int nb_integration_pts)
MoFEMErrorCode H1_BubbleShapeFunctions_ONSEGMENT(int p, double *N, double *diffN, double *bubbleN, double *diff_bubbleN, int nb_integration_pts)
MoFEMErrorCode H1_InteriorShapeFunctions_ONHEX(const int *p, double *N, double *N_diff, double *faceN, double *diff_faceN, int nb_integration_pts)
MoFEMErrorCode Hdiv_InteriorShapeFunctions_ONHEX(int *p, double *N, double *N_diff, double *bubleN[3], double *div_bubleN[3], int nb_integration_pts)
MoFEMErrorCode Hdiv_FaceShapeFunctions_ONQUAD(int *face_nodes, int *p, double *N, double *diffN, double *faceN, double *diff_faceN, int nb_integration_pts)
MoFEMErrorCode Hcurl_InteriorShapeFunctions_ONHEX(int *p, double *N, double *N_diff, double *volN[3], double *diff_volN[3], int nb_integration_pts)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2FromPtr< 3, 3 >(double *ptr)
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
static constexpr int edges_conn[]