v0.14.0
Functions
MoFEM::DemkowiczHexAndQuad Namespace Reference

Functions

MoFEMErrorCode Legendre_polynomials01 (int p, double s, double *L)
 
MoFEMErrorCode Integrated_Legendre01 (int p, double s, double *L, double *diffL)
 
MoFEMErrorCode Face_orientMat (int *face_nodes, double orientMat[2][2])
 
MoFEMErrorCode H1_BubbleShapeFunctions_ONSEGMENT (int p, double *N, double *diffN, double *bubbleN, double *diff_bubbleN, int nb_integration_pts)
 
MoFEMErrorCode L2_ShapeFunctions_ONSEGMENT (int p, double *N, double *diffN, double *funN, double *funDiffN, 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. More...
 
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. More...
 
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. More...
 
MoFEMErrorCode Hcurl_EdgeShapeFunctions_ONQUAD (int *sense, int *p, double *N, double *diffN, double *edgeN[4], double *curl_edgeN[4], 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_ONQUAD (int *face_nodes, int *p, double *N, double *diffN, double *faceN, double *diff_faceN, int nb_integration_pts)
 
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_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 H1_InteriorShapeFunctions_ONHEX (const int *p, double *N, double *N_diff, double *faceN, double *diff_faceN, int nb_integration_pts)
 
MoFEMErrorCode L2_InteriorShapeFunctions_ONHEX (const int *p, double *N, double *N_diff, double *volN, double *diff_volN, int nb_integration_pts)
 
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 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 Hcurl_InteriorShapeFunctions_ONHEX (int *p, double *N, double *N_diff, double *volN[3], double *diff_volN[3], 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 Hdiv_InteriorShapeFunctions_ONHEX (int *p, double *N, double *N_diff, double *bubleN[3], double *div_bubleN[3], int nb_integration_pts)
 

Function Documentation

◆ Face_orientMat()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Face_orientMat ( int *  face_nodes,
double  orientMat[2][2] 
)

◆ H1_BubbleShapeFunctions_ONSEGMENT()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::H1_BubbleShapeFunctions_ONSEGMENT ( int  p,
double N,
double diffN,
double bubbleN,
double diff_bubbleN,
int  nb_integration_pts 
)

Definition at line 141 of file EdgeQuadHexPolynomials.cpp.

143  {
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  }
164 }

◆ H1_EdgeShapeFunctions_ONHEX()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::H1_EdgeShapeFunctions_ONHEX ( int *  sense,
int *  p,
double N,
double N_diff,
double edgeN[12],
double diff_edgeN[12],
int  nb_integration_pts 
)

Definition at line 707 of file EdgeQuadHexPolynomials.cpp.

709  {
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  }
793 }

◆ H1_EdgeShapeFunctions_ONQUAD()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::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.

Function generates hierarchical base of H1 comforting functions on a 2D Quad.

Parameters
sensearray of orientation of edges (take 1 or -1)
parray of orders (in each direction) of base functions
Narray vertex shape functions evaluated at each integration point
diffNderivatives of vertex shape functions
Returns
edgeN base functions on edges at qd pts
diff_edgeN derivatives of edge shape functions at qd pts
Parameters
nb_integration_ptsnumber of integration points
base_polynomialspolynomial base function (f.e. Legendre of Integrated Legendre)
Returns
error code

Definition at line 202 of file EdgeQuadHexPolynomials.cpp.

204  {
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  }
273 }

◆ H1_FaceShapeFunctions_ONHEX()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::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 at line 795 of file EdgeQuadHexPolynomials.cpp.

797  {
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];
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  }
903 }

◆ H1_FaceShapeFunctions_ONQUAD()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::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.

Function generates hierarchical base of H1 comforting functions on a 2D quad.

Parameters
face_nodesface nodes order
parray of orders (in each direction) of base functions
Narray vertex shape functions evaluated at each integration point
diffNderivatives of vertex shape functions
Returns
edgeN base functions on edges at qd pts
diff_edgeN derivatives of edge shape functions at qd pts
Parameters
nb_integration_ptsnumber of integration points
base_polynomialspolynomial base function (f.e. Legendre of Integrated Legendre)
Returns
error code

Definition at line 275 of file EdgeQuadHexPolynomials.cpp.

277  {
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 
288  int permute[(p[0] - 1) * (p[1] - 1)][3];
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  }
333 }

◆ H1_InteriorShapeFunctions_ONHEX()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::H1_InteriorShapeFunctions_ONHEX ( const int *  p,
double N,
double N_diff,
double faceN,
double diff_faceN,
int  nb_integration_pts 
)

Definition at line 905 of file EdgeQuadHexPolynomials.cpp.

907  {
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];
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  }
964 }

◆ Hcurl_EdgeShapeFunctions_ONHEX()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hcurl_EdgeShapeFunctions_ONHEX ( int *  sense,
int *  p,
double N,
double N_diff,
double edgeN[12],
double diff_edgeN[12],
int  nb_integration_pts 
)

Definition at line 1019 of file EdgeQuadHexPolynomials.cpp.

1021  {
1022 
1025 
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  }
1137 }

◆ Hcurl_EdgeShapeFunctions_ONQUAD()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hcurl_EdgeShapeFunctions_ONQUAD ( int *  sense,
int *  p,
double N,
double diffN,
double edgeN[4],
double curl_edgeN[4],
int  nb_integration_pts 
)

Definition at line 393 of file EdgeQuadHexPolynomials.cpp.

395  {
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  }
482 }

◆ Hcurl_FaceShapeFunctions_ONHEX()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::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 at line 1139 of file EdgeQuadHexPolynomials.cpp.

1141  {
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];
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 
1283 }

◆ Hcurl_FaceShapeFunctions_ONQUAD()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hcurl_FaceShapeFunctions_ONQUAD ( int *  face_nodes,
int *  p,
double N,
double diffN,
double faceN[],
double diff_faceN[],
int  nb_integration_pts 
)

Definition at line 484 of file EdgeQuadHexPolynomials.cpp.

486  {
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];
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  }
588 }

◆ Hcurl_InteriorShapeFunctions_ONHEX()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hcurl_InteriorShapeFunctions_ONHEX ( int *  p,
double N,
double N_diff,
double volN[3],
double diff_volN[3],
int  nb_integration_pts 
)

◆ Hdiv_FaceShapeFunctions_ONHEX()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::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 at line 1404 of file EdgeQuadHexPolynomials.cpp.

1406  {
1408 
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];
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 
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 
1542 }

◆ Hdiv_FaceShapeFunctions_ONQUAD()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hdiv_FaceShapeFunctions_ONQUAD ( int *  face_nodes,
int *  p,
double N,
double diffN,
double faceN,
double diff_faceN,
int  nb_integration_pts 
)

Definition at line 590 of file EdgeQuadHexPolynomials.cpp.

592  {
594 
595  const int nb_dofs = (p[0] * p[1]);
596 
597  if (nb_dofs > 0) {
598 
599  int permute[nb_dofs][3];
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 
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 
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 
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  }
678 }

◆ Hdiv_InteriorShapeFunctions_ONHEX()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Hdiv_InteriorShapeFunctions_ONHEX ( int *  p,
double N,
double N_diff,
double bubleN[3],
double div_bubleN[3],
int  nb_integration_pts 
)

Definition at line 1544 of file EdgeQuadHexPolynomials.cpp.

1546  {
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];
1563  rrr - 1);
1564  }
1565 
1569 
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 
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 
1666 }

◆ Integrated_Legendre01()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Integrated_Legendre01 ( int  p,
double  s,
double L,
double diffL 
)

◆ L2_FaceShapeFunctions_ONQUAD()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::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.

Function generates hierarchical base of H1 comforting functions on a 2D quad.

Parameters
parray of orders (in each direction) of base functions
Narray vertex shape functions evaluated at each integration point
diffNderivatives of vertex shape functions
Returns
edgeN base functions on edges at qd pts
diff_edgeN derivatives of edge shape functions at qd pts
Parameters
nb_integration_ptsnumber of integration points
base_polynomialspolynomial base function (f.e. Legendre of Integrated Legendre)
Returns
error code

Definition at line 335 of file EdgeQuadHexPolynomials.cpp.

337  {
339  const int nb_dofs = (p[0] + 1) * (p[1] + 1);
340  if (nb_dofs > 0) {
341 
342  int permute[nb_dofs][3];
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  }
391 }

◆ L2_InteriorShapeFunctions_ONHEX()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::L2_InteriorShapeFunctions_ONHEX ( const int *  p,
double N,
double N_diff,
double volN,
double diff_volN,
int  nb_integration_pts 
)

Definition at line 966 of file EdgeQuadHexPolynomials.cpp.

968  {
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];
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  }
1017 }

◆ L2_ShapeFunctions_ONSEGMENT()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::L2_ShapeFunctions_ONSEGMENT ( int  p,
double N,
double diffN,
double funN,
double funDiffN,
int  nb_integration_pts 
)

Definition at line 166 of file EdgeQuadHexPolynomials.cpp.

168  {
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  }
189 }

◆ Legendre_polynomials01()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Legendre_polynomials01 ( int  p,
double  s,
double L 
)
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
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
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
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
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
eta
double eta
Definition: free_surface.cpp:170
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
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
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
FTensor::Index< 'i', 3 >
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::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
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