v0.14.0
Loading...
Searching...
No Matches
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 {
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}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
Lobatto_polynomials
PetscErrorCode Lobatto_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Lobatto base functions .
Definition: base_functions.c:199
N
const int N
Definition: speed_test.cpp:3

◆ 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 {
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}
FTensor::d
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
MoFEM::edges_conn
static constexpr int edges_conn[]
Definition: EntityRefine.cpp:18

◆ 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 {
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}

◆ 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 {
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}
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

◆ 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
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}

◆ 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 {
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}
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
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

◆ 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
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 getFTensor2FromPtr<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}
a
constexpr double a
Definition: approx_sphere.cpp:30
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
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::getFTensor2FromPtr< 3, 3 >
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2FromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:959

◆ 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 {
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 = getFTensor2FromPtr<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}
MoFEM::getFTensor2FromPtr< 3, 2 >
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:948

◆ 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 {
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 = getFTensor2FromPtr<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}
m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65
eta
double eta
Definition: free_surface.cpp:170
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:173

◆ 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 {
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 = getFTensor2FromPtr<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}

◆ 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 {
1407 MoFEMFunctionBeginHot;
1408
1409 FTensor::Index<'i', 3> i;
1410 FTensor::Index<'j', 3> j;
1411 FTensor::Index<'k', 3> k;
1412
1413 constexpr int opposite_face_node[6][4] = {
1414
1415 {3, 2, 6, 7},
1416
1417 {0, 3, 7, 4},
1418
1419 {1, 0, 4, 5},
1420
1421 {2, 1, 5, 6},
1422
1423 {4, 7, 6, 5},
1424
1425 {0, 1, 2, 3}
1426
1427 };
1428
1429 for (int face = 0; face != 6; face++) {
1430
1431 const int nb_dofs =
1432 NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL(p[face], p[face]);
1433
1434 if (nb_dofs > 0) {
1435
1436 auto t_n = getFTensor1FromPtr<3>(faceN[face]);
1437 auto t_diff_n = getFTensor2FromPtr<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}
NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL
#define NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL(P, Q)
Definition: h1_hdiv_hcurl_l2.h:143
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
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

◆ 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 {
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 = getFTensor2FromPtr<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}
J
static Index< 'J', 3 > J
Definition: BasicFeTools.hpp:93

◆ 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 {
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 = getFTensor2FromPtr<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}

◆ 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 {
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}

◆ 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 {
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}

◆ 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 {
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}

◆ Legendre_polynomials01()

MoFEMErrorCode MoFEM::DemkowiczHexAndQuad::Legendre_polynomials01 ( int  p,
double  s,
double *  L 
)
Generated by Doxygen 1.9.5 and hosted at University of Glasgow