v0.14.0
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 
7 using namespace MoFEM;
8 
9 namespace DemkowiczHexAndQuad {
10 
11 MoFEMErrorCode 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 
103 static 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 */
141 MoFEMErrorCode 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 
166 MoFEMErrorCode 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 
202 MoFEMErrorCode 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 
275 MoFEMErrorCode 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 
335 MoFEMErrorCode 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 
393 MoFEMErrorCode 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 
484 MoFEMErrorCode 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 
590 MoFEMErrorCode 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 
707 MoFEMErrorCode 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 
795 MoFEMErrorCode 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 
905 MoFEMErrorCode 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 
966 MoFEMErrorCode 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 
1019 MoFEMErrorCode 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 
1139 MoFEMErrorCode 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 
1285 MoFEMErrorCode 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 
1404 MoFEMErrorCode 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_GENERAL(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 
1544 MoFEMErrorCode 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 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FTensor::Tensor1< double, 3 >
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
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
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
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
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:126
sdf.r
int r
Definition: sdf.py:8
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_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::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
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:182
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_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
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
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
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::getFTensor2HVecFromPtr< 3, 3 >
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:968
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::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
eta
double eta
Definition: free_surface.cpp:170
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::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
DemkowiczHexAndQuad
Definition: EdgeQuadHexPolynomials.cpp:9
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
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::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
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
FTensor::Index< 'i', 3 >
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
convert.n
n
Definition: convert.py:82
N
const int N
Definition: speed_test.cpp:3
DemkowiczHexAndQuad::monom_ordering
MoFEMErrorCode monom_ordering(int *perm, int p, int q, int r=0)
Definition: EdgeQuadHexPolynomials.cpp:11
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
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::getFTensor2HVecFromPtr< 3, 2 >
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2HVecFromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:957
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::edges_conn
static constexpr int edges_conn[]
Definition: EntityRefine.cpp:18
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
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GENERAL
#define NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GENERAL(P, Q)
Definition: h1_hdiv_hcurl_l2.h:143
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
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::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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow