v0.10.0
Macros | Functions
h1_hdiv_hcurl_l2.h File Reference

Functions to approximate hierarchical spaces. More...

Go to the source code of this file.

Macros

#define DEPRECATED
 
#define NBVOLUMETET_L2(P)   ((P + 1) * (P + 2) * (P + 3) / 6)
 Number of base functions on tetrahedron for L2 space. More...
 
#define NBFACETRI_L2(P)   ((P + 1) * (P + 2) / 2)
 Number of base functions on triangle for L2 space. More...
 
#define NBEDGE_L2(P)   (P + 1)
 Number of base functions on edge fro L2 space. More...
 
#define NBEDGE_H1(P)   (((P) > 1) ? (P - 1) : 0)
 Numer of base function on edge for H1 space. More...
 
#define NBFACETRI_H1(P)   (((P) > 2) ? ((P - 2) * (P - 1) / 2) : 0)
 Number of base function on triangle for H1 space. More...
 
#define NBFACEQUAD_H1(P)   (((P) > 1) ? ((P - 1) * (P - 1)) : 0)
 Number of base functions on quad for H1 space. More...
 
#define NBFACEQUAD_L2(P)   (((P) >= 0) ? (P + 1) * (P + 1) : 0)
 Number of base functions on quad for L2 space. More...
 
#define NBVOLUMETET_H1(P)   (((P) > 3) ? ((P - 3) * (P - 2) * (P - 1) / 6) : 0)
 Number of base functions on tetrahedron for H1 space. More...
 
#define NBVOLUMEPRISM_H1(P)   ((P > 3) ? ((P - 2) * (P - 2) * (P - 2)) : 0)
 Number of base functions on prism for H1 space. More...
 
#define NBEDGE_AINSWORTH_HCURL(P)   (((P) > 0) ? (P + 1) : 0)
 
#define NBFACETRI_AINSWORTH_EDGE_HCURL(P)   (((P) > 1) ? P - 1 : 0)
 
#define NBFACETRI_AINSWORTH_FACE_HCURL(P)   (((P) > 2) ? (P - 1) * (P - 2) : 0)
 
#define NBFACETRI_AINSWORTH_HCURL(P)   ((P > 1) ? ((P)-1) * (P + 1) : 0)
 
#define NBVOLUMETET_AINSWORTH_FACE_HCURL(P)    (((P) > 2) ? (2 * (P - 1) * (P - 2)) : 0)
 
#define NBVOLUMETET_AINSWORTH_TET_HCURL(P)    (((P) > 3) ? ((P - 3) * (P - 2) * (P - 1) / 2) : 0)
 
#define NBVOLUMETET_AINSWORTH_HCURL(P)    (((P) > 2) ? (P - 2) * (P - 1) * (P + 1) / 2 : 0)
 
#define NBEDGE_DEMKOWICZ_HCURL(P)   (((P) > 0) ? (P) : 0)
 
#define NBFACETRI_DEMKOWICZ_HCURL(P)   (((P) > 1) ? (P) * ((P)-1) : 0)
 
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)    (((P) > 2) ? ((P) * ((P)-1) * ((P)-2) / 2) : 0)
 
#define NBFACEQUAD_DEMKOWICZ_FAMILY_QUAD_HCURL(P, Q)    (((P) > 0 && (Q) > 1) ? P * (Q - 1) : 0)
 Number of base functions on quad for Hcurl space. More...
 
#define NBFACEQUAD_DEMKOWICZ_HCURL(P)    (2 * NBFACEQUAD_DEMKOWICZ_FAMILY_QUAD_HCURL(P, P))
 
#define NBEDGE_HDIV(P)   (0)
 
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)   (((P) > 0) ? (P) : 0)
 
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)   (((P) > 2) ? (P - 1) * (P - 2) / 2 : 0)
 
#define NBFACETRI_AINSWORTH_HDIV(P)   (((P) > 0) ? (P + 1) * (P + 2) / 2 : 0)
 
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)   (((P) > 1) ? (P - 1) : 0)
 
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)   (((P) > 2) ? (P - 1) * (P - 2) : 0)
 
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)    (((P) > 3) ? (P - 3) * (P - 2) * (P - 1) / 2 : 0)
 
#define NBVOLUMETET_AINSWORTH_HDIV(P)    (((P) > 1) ? (P - 1) * (P + 1) * (P + 2) / 2 : 0)
 
#define NBFACETRI_DEMKOWICZ_HDIV(P)   ((P > 0) ? (P) * (P + 1) / 2 : 0)
 
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)    (((P) > 1) ? (P) * (P - 1) * (P + 1) / 2 : 0)
 

Functions

PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTRI (int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 Get base functions on triangle for L2 space. More...
 
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET (int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 Get base functions on tetrahedron for L2 space. More...
 
DEPRECATED PetscErrorCode L2_ShapeFunctions_MBTRI (int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
DEPRECATED PetscErrorCode L2_ShapeFunctions_MBTET (int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode L2_VolumeShapeDiffMBTETinvJ (int base_p, int p, double *volume_diffN, double *invJac, double *volume_diffNinvJac, int GDIM)
 
PetscErrorCode H1_EdgeShapeFunctions_MBTRI (int *sense, int *p, double *N, double *diffN, double *edgeN[3], double *diff_edgeN[3], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 H1_EdgeShapeFunctions_MBTRI. More...
 
PetscErrorCode H1_FaceShapeFunctions_MBTRI (const int *face_nodes, int p, double *N, double *diffN, double *faceN, double *diff_faceN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_EdgeShapeFunctions_MBTET (int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_FaceShapeFunctions_MBTET (int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_VolumeShapeFunctions_MBTET (int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_EdgeShapeDiffMBTETinvJ (int *base_p, int *p, double *edge_diffN[], double *invJac, double *edge_diffNinvJac[], int GDIM)
 
PetscErrorCode H1_FaceShapeDiffMBTETinvJ (int *base_p, int *p, double *face_diffN[], double *invJac, double *face_diffNinvJac[], int GDIM)
 
PetscErrorCode H1_VolumeShapeDiffMBTETinvJ (int base_p, int p, double *volume_diffN, double *invJac, double *volume_diffNinvJac, int GDIM)
 
PetscErrorCode H1_EdgeGradientOfDeformation_hierarchical (int p, double *diffN, double *dofs, double *F)
 
PetscErrorCode H1_FaceGradientOfDeformation_hierarchical (int p, double *diffN, double *dofs, double *F)
 
PetscErrorCode H1_VolumeGradientOfDeformation_hierarchical (int p, double *diffN, double *dofs, double *F)
 
DEPRECATED PetscErrorCode H1_EdgeGradientOfDeformation_hierachical (int p, double *diffN, double *dofs, double *F)
 
DEPRECATED PetscErrorCode H1_FaceGradientOfDeformation_hierachical (int p, double *diffN, double *dofs, double *F)
 
DEPRECATED PetscErrorCode H1_VolumeGradientOfDeformation_hierachical (int p, double *diffN, double *dofs, double *F)
 
PetscErrorCode H1_QuadShapeFunctions_MBPRISM (int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_VolumeShapeFunctions_MBPRISM (int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_QuadShapeFunctions_MBQUAD (int *faces_nodes, int p, double *N, double *diffN, double *faceN, double *diff_faceN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_EdgeShapeFunctions_MBQUAD (int *sense, int *p, double *N, double *diffN, double *edgeN[4], double *diff_edgeN[4], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 

Detailed Description

Functions to approximate hierarchical spaces.

\FIXME: Name Shape Functions is used, in that context is more appropriate to use base functions. Need to be changed.

Definition in file h1_hdiv_hcurl_l2.h.

Macro Definition Documentation

◆ DEPRECATED

#define DEPRECATED

Definition at line 23 of file h1_hdiv_hcurl_l2.h.

◆ NBEDGE_AINSWORTH_HCURL

#define NBEDGE_AINSWORTH_HCURL (   P)    (((P) > 0) ? (P + 1) : 0)

Definition at line 84 of file h1_hdiv_hcurl_l2.h.

◆ NBEDGE_DEMKOWICZ_HCURL

#define NBEDGE_DEMKOWICZ_HCURL (   P)    (((P) > 0) ? (P) : 0)
Examples
bases_on_reference_rectangle.cpp.

Definition at line 95 of file h1_hdiv_hcurl_l2.h.

◆ NBEDGE_H1

#define NBEDGE_H1 (   P)    (((P) > 1) ? (P - 1) : 0)

Numer of base function on edge for H1 space.

Examples
bernstein_bezier_generate_base.cpp, and edge_and_bubble_shape_functions_on_quad.cpp.

Definition at line 55 of file h1_hdiv_hcurl_l2.h.

◆ NBEDGE_HDIV

#define NBEDGE_HDIV (   P)    (0)

Definition at line 110 of file h1_hdiv_hcurl_l2.h.

◆ NBEDGE_L2

#define NBEDGE_L2 (   P)    (P + 1)

Number of base functions on edge fro L2 space.

Definition at line 48 of file h1_hdiv_hcurl_l2.h.

◆ NBFACEQUAD_DEMKOWICZ_FAMILY_QUAD_HCURL

#define NBFACEQUAD_DEMKOWICZ_FAMILY_QUAD_HCURL (   P,
 
)     (((P) > 0 && (Q) > 1) ? P * (Q - 1) : 0)

Number of base functions on quad for Hcurl space.

Examples
bases_on_reference_rectangle.cpp.

Definition at line 103 of file h1_hdiv_hcurl_l2.h.

◆ NBFACEQUAD_DEMKOWICZ_HCURL

#define NBFACEQUAD_DEMKOWICZ_HCURL (   P)     (2 * NBFACEQUAD_DEMKOWICZ_FAMILY_QUAD_HCURL(P, P))
Examples
bases_on_reference_rectangle.cpp.

Definition at line 105 of file h1_hdiv_hcurl_l2.h.

◆ NBFACEQUAD_H1

#define NBFACEQUAD_H1 (   P)    (((P) > 1) ? ((P - 1) * (P - 1)) : 0)

Number of base functions on quad for H1 space.

Examples
edge_and_bubble_shape_functions_on_quad.cpp.

Definition at line 65 of file h1_hdiv_hcurl_l2.h.

◆ NBFACEQUAD_L2

#define NBFACEQUAD_L2 (   P)    (((P) >= 0) ? (P + 1) * (P + 1) : 0)

Number of base functions on quad for L2 space.

Definition at line 70 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_AINSWORTH_EDGE_HCURL

#define NBFACETRI_AINSWORTH_EDGE_HCURL (   P)    (((P) > 1) ? P - 1 : 0)

Definition at line 85 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_AINSWORTH_EDGE_HDIV

#define NBFACETRI_AINSWORTH_EDGE_HDIV (   P)    (((P) > 0) ? (P) : 0)

Definition at line 111 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_AINSWORTH_FACE_HCURL

#define NBFACETRI_AINSWORTH_FACE_HCURL (   P)    (((P) > 2) ? (P - 1) * (P - 2) : 0)

Definition at line 86 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_AINSWORTH_FACE_HDIV

#define NBFACETRI_AINSWORTH_FACE_HDIV (   P)    (((P) > 2) ? (P - 1) * (P - 2) / 2 : 0)

Definition at line 112 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_AINSWORTH_HCURL

#define NBFACETRI_AINSWORTH_HCURL (   P)    ((P > 1) ? ((P)-1) * (P + 1) : 0)

Definition at line 87 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_AINSWORTH_HDIV

#define NBFACETRI_AINSWORTH_HDIV (   P)    (((P) > 0) ? (P + 1) * (P + 2) / 2 : 0)

Definition at line 113 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_DEMKOWICZ_HCURL

#define NBFACETRI_DEMKOWICZ_HCURL (   P)    (((P) > 1) ? (P) * ((P)-1) : 0)

Definition at line 96 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_DEMKOWICZ_HDIV

#define NBFACETRI_DEMKOWICZ_HDIV (   P)    ((P > 0) ? (P) * (P + 1) / 2 : 0)

◆ NBFACETRI_H1

#define NBFACETRI_H1 (   P)    (((P) > 2) ? ((P - 2) * (P - 1) / 2) : 0)

Number of base function on triangle for H1 space.

Examples
bernstein_bezier_generate_base.cpp.

Definition at line 60 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_L2

#define NBFACETRI_L2 (   P)    ((P + 1) * (P + 2) / 2)

Number of base functions on triangle for L2 space.

Definition at line 42 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMEPRISM_H1

#define NBVOLUMEPRISM_H1 (   P)    ((P > 3) ? ((P - 2) * (P - 2) * (P - 2)) : 0)

Number of base functions on prism for H1 space.

Definition at line 80 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMETET_AINSWORTH_EDGE_HDIV

#define NBVOLUMETET_AINSWORTH_EDGE_HDIV (   P)    (((P) > 1) ? (P - 1) : 0)

Definition at line 114 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMETET_AINSWORTH_FACE_HCURL

#define NBVOLUMETET_AINSWORTH_FACE_HCURL (   P)     (((P) > 2) ? (2 * (P - 1) * (P - 2)) : 0)

Definition at line 88 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMETET_AINSWORTH_FACE_HDIV

#define NBVOLUMETET_AINSWORTH_FACE_HDIV (   P)    (((P) > 2) ? (P - 1) * (P - 2) : 0)

Definition at line 115 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMETET_AINSWORTH_HCURL

#define NBVOLUMETET_AINSWORTH_HCURL (   P)     (((P) > 2) ? (P - 2) * (P - 1) * (P + 1) / 2 : 0)

Definition at line 92 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMETET_AINSWORTH_HDIV

#define NBVOLUMETET_AINSWORTH_HDIV (   P)     (((P) > 1) ? (P - 1) * (P + 1) * (P + 2) / 2 : 0)

Definition at line 118 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMETET_AINSWORTH_TET_HCURL

#define NBVOLUMETET_AINSWORTH_TET_HCURL (   P)     (((P) > 3) ? ((P - 3) * (P - 2) * (P - 1) / 2) : 0)

Definition at line 90 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMETET_AINSWORTH_VOLUME_HDIV

#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV (   P)     (((P) > 3) ? (P - 3) * (P - 2) * (P - 1) / 2 : 0)

Definition at line 116 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMETET_DEMKOWICZ_HCURL

#define NBVOLUMETET_DEMKOWICZ_HCURL (   P)     (((P) > 2) ? ((P) * ((P)-1) * ((P)-2) / 2) : 0)

Definition at line 97 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMETET_DEMKOWICZ_HDIV

#define NBVOLUMETET_DEMKOWICZ_HDIV (   P)     (((P) > 1) ? (P) * (P - 1) * (P + 1) / 2 : 0)

◆ NBVOLUMETET_H1

#define NBVOLUMETET_H1 (   P)    (((P) > 3) ? ((P - 3) * (P - 2) * (P - 1) / 6) : 0)

Number of base functions on tetrahedron for H1 space.

Examples
bernstein_bezier_generate_base.cpp.

Definition at line 75 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMETET_L2

#define NBVOLUMETET_L2 (   P)    ((P + 1) * (P + 2) * (P + 3) / 6)

Number of base functions on tetrahedron for L2 space.

Definition at line 37 of file h1_hdiv_hcurl_l2.h.

Function Documentation

◆ H1_EdgeGradientOfDeformation_hierachical()

DEPRECATED PetscErrorCode H1_EdgeGradientOfDeformation_hierachical ( int  p,
double diffN,
double dofs,
double F 
)
Deprecated:
use H1_EdgeGradientOfDeformation_hierarchical

Definition at line 1220 of file h1.c.

1222  {
1223  return H1_EdgeGradientOfDeformation_hierarchical(p, diffN, dofs, F);
1224 }

◆ H1_EdgeGradientOfDeformation_hierarchical()

PetscErrorCode H1_EdgeGradientOfDeformation_hierarchical ( int  p,
double diffN,
double dofs,
double F 
)

Definition at line 553 of file h1.c.

555  {
557  int col, row = 0;
558  for (; row < 3; row++)
559  for (col = 0; col < 3; col++)
560  F[3 * row + col] =
561  cblas_ddot(NBEDGE_H1(p), &diffN[col], 3, &dofs[row], 3);
563 }

◆ H1_EdgeShapeDiffMBTETinvJ()

PetscErrorCode H1_EdgeShapeDiffMBTETinvJ ( int *  base_p,
int *  p,
double edge_diffN[],
double invJac,
double edge_diffNinvJac[],
int  GDIM 
)

Definition at line 500 of file h1.c.

502  {
504  int ee = 0, ii, gg;
505  for (; ee < 6; ee++) {
506  for (ii = 0; ii < NBEDGE_H1(p[ee]); ii++) {
507  for (gg = 0; gg < GDIM; gg++) {
508  int shift1 = NBEDGE_H1(base_p[ee]) * gg;
509  int shift2 = NBEDGE_H1(p[ee]) * gg;
511  &(edge_diffN[ee])[3 * shift1 + 3 * ii], 1, 0.,
512  &(edge_diffNinvJac[ee])[3 * shift2 + 3 * ii], 1);
513  }
514  }
515  }
517 }

◆ H1_EdgeShapeFunctions_MBQUAD()

PetscErrorCode H1_EdgeShapeFunctions_MBQUAD ( int *  sense,
int *  p,
double N,
double diffN,
double edgeN[4],
double diff_edgeN[4],
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)
Examples
edge_and_bubble_shape_functions_on_quad.cpp.

Definition at line 1035 of file h1.c.

1040  {
1042 
1043  double *edgeN01 = NULL, *edgeN12 = NULL, *edgeN23 = NULL, *edgeN30 = NULL;
1044  if (edgeN != NULL) {
1045  edgeN01 = edgeN[0];
1046  edgeN12 = edgeN[1];
1047  edgeN23 = edgeN[2];
1048  edgeN30 = edgeN[3];
1049  }
1050  double *diff_edgeN01 = NULL, *diff_edgeN12 = NULL, *diff_edgeN23 = NULL,
1051  *diff_edgeN30 = NULL;
1052  if (diff_edgeN != NULL) {
1053  diff_edgeN01 = diff_edgeN[0];
1054  diff_edgeN12 = diff_edgeN[1];
1055  diff_edgeN23 = diff_edgeN[2];
1056  diff_edgeN30 = diff_edgeN[3];
1057  }
1058  int P[4];
1059  int ee = 0;
1060  for (; ee < 4; ee++)
1061  P[ee] = NBEDGE_H1(p[ee]);
1062 
1063  int n0 = 0;
1064  int n1 = 1;
1065  int n2 = 2;
1066  int n3 = 3;
1067 
1068  int ii = 0;
1069  for (; ii < GDIM; ii++) {
1070  int node_shift = ii * 4;
1071  int node_diff_shift = 2 * node_shift;
1072 
1073  double shape0 = N[node_shift + n0];
1074  double shape1 = N[node_shift + n1];
1075  double shape2 = N[node_shift + n2];
1076  double shape3 = N[node_shift + n3];
1077 
1078  double ksi01 = (shape1 + shape2 - shape0 - shape3) * sense[n0];
1079  double ksi12 = (shape2 + shape3 - shape1 - shape0) * sense[n1];
1080  double ksi23 = (shape3 + shape0 - shape2 - shape1) * sense[n2];
1081  double ksi30 = (shape0 + shape1 - shape3 - shape2) * sense[n3];
1082 
1083  double extrude_zeta01 = shape0 + shape1;
1084  double extrude_ksi12 = shape1 + shape2;
1085  double extrude_zeta23 = shape2 + shape3;
1086  double extrude_ksi30 = shape0 + shape3;
1087 
1088  double bubble_ksi = extrude_ksi12 * extrude_ksi30;
1089  double bubble_zeta = extrude_zeta01 * extrude_zeta23;
1090 
1091  double diff_ksi01[2], diff_ksi12[2], diff_ksi23[2], diff_ksi30[2];
1092  double diff_extrude_zeta01[2];
1093  double diff_extrude_ksi12[2];
1094  double diff_extrude_zeta23[2];
1095  double diff_extrude_ksi30[2];
1096  double diff_bubble_ksi[2];
1097  double diff_bubble_zeta[2];
1098 
1099  int d = 0;
1100  for (; d < 2; d++) {
1101  double diff_shape0 = diffN[node_diff_shift + 2 * n0 + d];
1102  double diff_shape1 = diffN[node_diff_shift + 2 * n1 + d];
1103  double diff_shape2 = diffN[node_diff_shift + 2 * n2 + d];
1104  double diff_shape3 = diffN[node_diff_shift + 2 * n3 + d];
1105  diff_ksi01[d] =
1106  (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) * sense[n0];
1107  diff_ksi12[d] =
1108  (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) * sense[n1];
1109  diff_ksi23[d] =
1110  (diff_shape3 + diff_shape0 - diff_shape2 - diff_shape1) * sense[n2];
1111  diff_ksi30[d] =
1112  (diff_shape0 + diff_shape1 - diff_shape3 - diff_shape2) * sense[n3];
1113  diff_extrude_zeta01[d] = diff_shape0 + diff_shape1;
1114  diff_extrude_ksi12[d] = diff_shape1 + diff_shape2;
1115  diff_extrude_zeta23[d] = diff_shape2 + diff_shape3;
1116  diff_extrude_ksi30[d] = diff_shape0 + diff_shape3;
1117  diff_bubble_ksi[d] = diff_extrude_ksi12[d] * extrude_ksi30 +
1118  extrude_ksi12 * diff_extrude_ksi30[d];
1119  diff_bubble_zeta[d] = diff_extrude_zeta01[d] * extrude_zeta23 +
1120  extrude_zeta01 * diff_extrude_zeta23[d];
1121  }
1122 
1123  double L01[p[0] + 1], L12[p[1] + 1], L23[p[2] + 1], L30[p[3] + 1];
1124  double diffL01[2 * (p[0] + 1)], diffL12[2 * (p[1] + 1)],
1125  diffL23[2 * (p[2] + 1)], diffL30[2 * (p[3] + 1)];
1126  ierr = base_polynomials(p[0], ksi01, diff_ksi01, L01, diffL01, 2);
1127  CHKERRQ(ierr);
1128  ierr = base_polynomials(p[1], ksi12, diff_ksi12, L12, diffL12, 2);
1129  CHKERRQ(ierr);
1130  ierr = base_polynomials(p[2], ksi23, diff_ksi23, L23, diffL23, 2);
1131  CHKERRQ(ierr);
1132  ierr = base_polynomials(p[3], ksi30, diff_ksi30, L30, diffL30, 2);
1133  CHKERRQ(ierr);
1134 
1135  int shift;
1136  if (edgeN != NULL) {
1137  // edge01
1138  shift = ii * (P[0]);
1139  cblas_dcopy(P[0], L01, 1, &edgeN01[shift], 1);
1140  cblas_dscal(P[0], bubble_ksi * extrude_zeta01, &edgeN01[shift], 1);
1141  // edge12
1142  shift = ii * (P[1]);
1143  cblas_dcopy(P[1], L12, 1, &edgeN12[shift], 1);
1144  cblas_dscal(P[1], bubble_zeta * extrude_ksi12, &edgeN12[shift], 1);
1145  // edge23
1146  shift = ii * (P[2]);
1147  cblas_dcopy(P[2], L23, 1, &edgeN23[shift], 1);
1148  cblas_dscal(P[2], bubble_ksi * extrude_zeta23, &edgeN23[shift], 1);
1149  // edge30
1150  shift = ii * (P[3]);
1151  cblas_dcopy(P[3], L30, 1, &edgeN30[shift], 1);
1152  cblas_dscal(P[3], bubble_zeta * extrude_ksi30, &edgeN30[shift], 1);
1153  }
1154  if (diff_edgeN != NULL) {
1155  if (P[0] > 0) {
1156  // edge01
1157  shift = ii * (P[0]);
1158  bzero(&diff_edgeN01[2 * shift], sizeof(double) * 2 * (P[0]));
1159  int d = 0;
1160  for (; d != 2; ++d) {
1161  cblas_daxpy(P[0], bubble_ksi * extrude_zeta01,
1162  &diffL01[d * (p[0] + 1)], 1, &diff_edgeN01[2 * shift + d],
1163  2);
1164  cblas_daxpy(P[0],
1165  diff_bubble_ksi[d] * extrude_zeta01 +
1166  bubble_ksi * diff_extrude_zeta01[d],
1167  L01, 1, &diff_edgeN01[2 * shift + d], 2);
1168  }
1169  }
1170  if (P[1] > 0) {
1171  // edge12
1172  shift = ii * (P[1]);
1173  bzero(&diff_edgeN12[2 * shift], sizeof(double) * 2 * (P[1]));
1174  int d = 0;
1175  for (; d != 2; ++d) {
1176  cblas_daxpy(P[1], bubble_zeta * extrude_ksi12,
1177  &diffL12[d * (p[1] + 1)], 1, &diff_edgeN12[2 * shift + d],
1178  2);
1179  cblas_daxpy(P[1],
1180  diff_bubble_zeta[d] * extrude_ksi12 +
1181  bubble_zeta * diff_extrude_ksi12[d],
1182  L12, 1, &diff_edgeN12[2 * shift + d], 2);
1183  }
1184  }
1185  if (P[2] > 0) {
1186  // edge23
1187  shift = ii * (P[2]);
1188  bzero(&diff_edgeN23[2 * shift], sizeof(double) * 2 * (P[2]));
1189  int d = 0;
1190  for (; d != 2; ++d) {
1191  cblas_daxpy(P[2], bubble_ksi * extrude_zeta23,
1192  &diffL23[d * (p[2] + 1)], 1, &diff_edgeN23[2 * shift + d],
1193  2);
1194  cblas_daxpy(P[2],
1195  diff_bubble_ksi[d] * extrude_zeta23 +
1196  bubble_ksi * diff_extrude_zeta23[d],
1197  L23, 1, &diff_edgeN23[2 * shift + d], 2);
1198  }
1199  }
1200  if (P[3] > 0) {
1201  // edge30
1202  shift = ii * (P[3]);
1203  bzero(&diff_edgeN30[2 * shift], sizeof(double) * 2 * (P[3]));
1204  int d = 0;
1205  for (; d != 2; ++d) {
1206  cblas_daxpy(P[3], bubble_zeta * extrude_ksi30,
1207  &diffL30[d * (p[3] + 1)], 1, &diff_edgeN30[2 * shift + d],
1208  2);
1209  cblas_daxpy(P[3],
1210  diff_bubble_zeta[d] * extrude_ksi30 +
1211  bubble_zeta * diff_extrude_ksi30[d],
1212  L30, 1, &diff_edgeN30[2 * shift + d], 2);
1213  }
1214  }
1215  }
1216  }
1218 }

◆ H1_EdgeShapeFunctions_MBTET()

PetscErrorCode H1_EdgeShapeFunctions_MBTET ( int *  sense,
int *  p,
double N,
double diffN,
double edgeN[],
double diff_edgeN[],
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)

Definition at line 245 of file h1.c.

250  {
252 
253  int P[6];
254  int ee = 0;
255  for (; ee < 6; ee++)
256  P[ee] = NBEDGE_H1(p[ee]);
257  int edges_nodes[2 * 6] = {0, 1, 1, 2, 2, 0, 0, 3, 1, 3, 2, 3};
258  double diff_ksi01[3], diff_ksi12[3], diff_ksi20[3];
259  double diff_ksi03[3], diff_ksi13[3], diff_ksi23[3];
260  double *edges_diff_ksi[6] = {diff_ksi01, diff_ksi12, diff_ksi20,
261  diff_ksi03, diff_ksi13, diff_ksi23};
262  for (ee = 0; ee < 6; ee++) {
263  int dd = 0;
264  for (; dd < 3; dd++) {
265  edges_diff_ksi[ee][dd] = (diffN[edges_nodes[2 * ee + 1] * 3 + dd] -
266  diffN[edges_nodes[2 * ee + 0] * 3 + dd]) *
267  sense[ee];
268  }
269  }
270  int ii = 0;
271  for (; ii < GDIM; ii++) {
272  int node_shift = ii * 4;
273  double edges_ksi[6];
274  for (ee = 0; ee < 6; ee++) {
275  edges_ksi[ee] = (N[node_shift + edges_nodes[2 * ee + 1]] -
276  N[node_shift + edges_nodes[2 * ee + 0]]) *
277  sense[ee];
278  }
279  for (ee = 0; ee < 6; ee++) {
280  if (P[ee] == 0)
281  continue;
282  double L[p[ee] + 1], diffL[3 * (p[ee] + 1)];
283  ierr = base_polynomials(p[ee], edges_ksi[ee], edges_diff_ksi[ee], L,
284  diffL, 3);
285  CHKERRQ(ierr);
286  double v = N[node_shift + edges_nodes[2 * ee + 0]] *
287  N[node_shift + edges_nodes[2 * ee + 1]];
288  if (edgeN != NULL)
289  if (edgeN[ee] != NULL) {
290  int shift = ii * P[ee];
291  cblas_dcopy(P[ee], L, 1, &edgeN[ee][shift], 1);
292  cblas_dscal(P[ee],
293  N[node_shift + edges_nodes[2 * ee + 0]] *
294  N[node_shift + edges_nodes[2 * ee + 1]],
295  &edgeN[ee][shift], 1);
296  }
297  if (diff_edgeN != NULL)
298  if (diff_edgeN[ee] != NULL) {
299  int shift = ii * P[ee];
300  bzero(&diff_edgeN[ee][3 * shift], sizeof(double) * 3 * P[ee]);
301  int dd = 0;
302  for (; dd < 3; dd++) {
303  cblas_daxpy(P[ee], v, &diffL[dd * (p[ee] + 1)], 1,
304  &diff_edgeN[ee][3 * shift + dd], 3);
305  cblas_daxpy(P[ee],
306  diffN[3 * edges_nodes[2 * ee + 0] + dd] *
307  N[node_shift + edges_nodes[2 * ee + 1]] +
308  N[node_shift + edges_nodes[2 * ee + 0]] *
309  diffN[3 * edges_nodes[2 * ee + 1] + dd],
310  L, 1, &diff_edgeN[ee][3 * shift + dd], 3);
311  }
312  }
313  }
314  }
316 }

◆ H1_EdgeShapeFunctions_MBTRI()

PetscErrorCode H1_EdgeShapeFunctions_MBTRI ( int *  sense,
int *  p,
double N,
double diffN,
double edgeN[3],
double diff_edgeN[3],
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)

H1_EdgeShapeFunctions_MBTRI.

Parameters
senseof edges, it is array of integers dim 3 (3-edges of triangle)
pof edges

Definition at line 27 of file h1.c.

32  {
34 
35  double *edgeN01 = NULL, *edgeN12 = NULL, *edgeN20 = NULL;
36  if (edgeN != NULL) {
37  edgeN01 = edgeN[0];
38  edgeN12 = edgeN[1];
39  edgeN20 = edgeN[2];
40  }
41  double *diff_edgeN01 = NULL, *diff_edgeN12 = NULL, *diff_edgeN20 = NULL;
42  if (diff_edgeN != NULL) {
43  diff_edgeN01 = diff_edgeN[0];
44  diff_edgeN12 = diff_edgeN[1];
45  diff_edgeN20 = diff_edgeN[2];
46  }
47  int P[3];
48  int ee = 0;
49  for (; ee < 3; ee++) {
50  P[ee] = NBEDGE_H1(p[ee]);
51  }
52  int dd = 0;
53  double diff_ksi01[2], diff_ksi12[2], diff_ksi20[2];
54  if (diff_edgeN != NULL) {
55  for (; dd < 2; dd++) {
56  diff_ksi01[dd] = (diffN[1 * 2 + dd] - diffN[0 * 2 + dd]) * sense[0];
57  diff_ksi12[dd] = (diffN[2 * 2 + dd] - diffN[1 * 2 + dd]) * sense[1];
58  diff_ksi20[dd] = (diffN[0 * 2 + dd] - diffN[2 * 2 + dd]) * sense[2];
59  }
60  }
61  int ii = 0;
62  for (; ii < GDIM; ii++) {
63  int node_shift = ii * 3;
64  double ksi01 = (N[node_shift + 1] - N[node_shift + 0]) * sense[0];
65  double ksi12 = (N[node_shift + 2] - N[node_shift + 1]) * sense[1];
66  double ksi20 = (N[node_shift + 0] - N[node_shift + 2]) * sense[2];
67  double L01[p[0] + 1], L12[p[1] + 1], L20[p[2] + 1];
68  double diffL01[2 * (p[0] + 1)], diffL12[2 * (p[1] + 1)],
69  diffL20[2 * (p[2] + 1)];
70  ierr = base_polynomials(p[0], ksi01, diff_ksi01, L01, diffL01, 2);
71  CHKERRQ(ierr);
72  ierr = base_polynomials(p[1], ksi12, diff_ksi12, L12, diffL12, 2);
73  CHKERRQ(ierr);
74  ierr = base_polynomials(p[2], ksi20, diff_ksi20, L20, diffL20, 2);
75  CHKERRQ(ierr);
76  int shift;
77  if (edgeN != NULL) {
78  // edge01
79  shift = ii * (P[0]);
80  cblas_dcopy(P[0], L01, 1, &edgeN01[shift], 1);
81  cblas_dscal(P[0], N[node_shift + 0] * N[node_shift + 1], &edgeN01[shift],
82  1);
83  // edge12
84  shift = ii * (P[1]);
85  cblas_dcopy(P[1], L12, 1, &edgeN12[shift], 1);
86  cblas_dscal(P[1], N[node_shift + 1] * N[node_shift + 2], &edgeN12[shift],
87  1);
88  // edge20
89  shift = ii * (P[2]);
90  cblas_dcopy(P[2], L20, 1, &edgeN20[shift], 1);
91  cblas_dscal(P[2], N[node_shift + 0] * N[node_shift + 2], &edgeN20[shift],
92  1);
93  }
94  if (diff_edgeN != NULL) {
95  if (P[0] > 0) {
96  // edge01
97  shift = ii * (P[0]);
98  bzero(&diff_edgeN01[2 * shift], sizeof(double) * 2 * (P[0]));
99  // diffX
100  cblas_daxpy(P[0], N[node_shift + 0] * N[node_shift + 1],
101  &diffL01[0 * (p[0] + 1)], 1, &diff_edgeN01[2 * shift + 0],
102  2);
103  cblas_daxpy(P[0],
104  diffN[2 * 0 + 0] * N[node_shift + 1] +
105  N[node_shift + 0] * diffN[2 * 1 + 0],
106  L01, 1, &diff_edgeN01[2 * shift + 0], 2);
107  // diff Y
108  cblas_daxpy(P[0], N[node_shift + 0] * N[node_shift + 1],
109  &diffL01[1 * (p[0] + 1)], 1, &diff_edgeN01[2 * shift + 1],
110  2);
111  cblas_daxpy(P[0],
112  diffN[2 * 0 + 1] * N[node_shift + 1] +
113  N[node_shift + 0] * diffN[2 * 1 + 1],
114  L01, 1, &diff_edgeN01[2 * shift + 1], 2);
115  }
116  if (P[1] > 0) {
117  // edge12
118  shift = ii * (P[1]);
119  bzero(&diff_edgeN12[2 * shift], sizeof(double) * 2 * (P[1]));
120  // diffX
121  cblas_daxpy(P[1], N[node_shift + 1] * N[node_shift + 2],
122  &diffL12[0 * (p[1] + 1)], 1, &diff_edgeN12[2 * shift + 0],
123  2);
124  cblas_daxpy(P[1],
125  diffN[2 * 1 + 0] * N[node_shift + 2] +
126  N[node_shift + 1] * diffN[2 * 2 + 0],
127  L12, 1, &diff_edgeN12[2 * shift + 0], 2);
128  // diffY
129  cblas_daxpy(P[1], N[node_shift + 1] * N[node_shift + 2],
130  &diffL12[1 * (p[1] + 1)], 1, &diff_edgeN12[2 * shift + 1],
131  2);
132  cblas_daxpy(P[1],
133  diffN[2 * 1 + 1] * N[node_shift + 2] +
134  N[node_shift + 1] * diffN[2 * 2 + 1],
135  L12, 1, &diff_edgeN12[2 * shift + 1], 2);
136  }
137  if (P[2] > 0) {
138  // edge20
139  shift = ii * (P[2]);
140  bzero(&diff_edgeN20[2 * shift], sizeof(double) * 2 * (P[2]));
141  // diffX
142  cblas_daxpy(P[2], N[node_shift + 0] * N[node_shift + 2],
143  &diffL20[0 * (p[2] + 1)], 1, &diff_edgeN20[2 * shift + 0],
144  2);
145  cblas_daxpy(P[2],
146  diffN[2 * 0 + 0] * N[node_shift + 2] +
147  N[node_shift + 0] * diffN[2 * 2 + 0],
148  L20, 1, &diff_edgeN20[2 * shift + 0], 2);
149  // diffY
150  cblas_daxpy(P[2], N[node_shift + 0] * N[node_shift + 2],
151  &diffL20[1 * (p[2] + 1)], 1, &diff_edgeN20[2 * shift + 1],
152  2);
153  cblas_daxpy(P[2],
154  diffN[2 * 0 + 1] * N[node_shift + 2] +
155  N[node_shift + 0] * diffN[2 * 2 + 1],
156  L20, 1, &diff_edgeN20[2 * shift + 1], 2);
157  }
158  }
159  }
161 }

◆ H1_FaceGradientOfDeformation_hierachical()

DEPRECATED PetscErrorCode H1_FaceGradientOfDeformation_hierachical ( int  p,
double diffN,
double dofs,
double F 
)
Deprecated:
use H1_FaceGradientOfDeformation_hierarchical

Definition at line 1225 of file h1.c.

1227  {
1228  return H1_FaceGradientOfDeformation_hierarchical(p, diffN, dofs, F);
1229 }

◆ H1_FaceGradientOfDeformation_hierarchical()

PetscErrorCode H1_FaceGradientOfDeformation_hierarchical ( int  p,
double diffN,
double dofs,
double F 
)

Definition at line 564 of file h1.c.

566  {
568  int col, row = 0;
569  for (; row < 3; row++)
570  for (col = 0; col < 3; col++)
571  F[3 * row + col] =
572  cblas_ddot(NBFACETRI_H1(p), &diffN[col], 3, &dofs[row], 3);
574 }

◆ H1_FaceShapeDiffMBTETinvJ()

PetscErrorCode H1_FaceShapeDiffMBTETinvJ ( int *  base_p,
int *  p,
double face_diffN[],
double invJac,
double face_diffNinvJac[],
int  GDIM 
)

Definition at line 518 of file h1.c.

520  {
522  int ff = 0, ii, gg;
523  for (; ff < 4; ff++) {
524  for (ii = 0; ii < NBFACETRI_H1(p[ff]); ii++) {
525  for (gg = 0; gg < GDIM; gg++) {
526  int shift1 = NBFACETRI_H1(base_p[ff]) * gg;
527  int shift2 = NBFACETRI_H1(p[ff]) * gg;
529  &(face_diffN[ff])[3 * shift1 + 3 * ii], 1, 0.,
530  &(face_diffNinvJac[ff])[3 * shift2 + 3 * ii], 1);
531  }
532  }
533  }
535 }

◆ H1_FaceShapeFunctions_MBTET()

PetscErrorCode H1_FaceShapeFunctions_MBTET ( int *  faces_nodes,
int *  p,
double N,
double diffN,
double faceN[],
double diff_faceN[],
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)

Definition at line 317 of file h1.c.

322  {
324 
325  int P[4];
326  int ff = 0;
327  for (; ff < 4; ff++)
328  P[ff] = NBFACETRI_H1(p[ff]);
329  double diff_ksiL0F0[3], diff_ksiL1F0[3];
330  double diff_ksiL0F1[3], diff_ksiL1F1[3];
331  double diff_ksiL0F2[3], diff_ksiL1F2[3];
332  double diff_ksiL0F3[3], diff_ksiL1F3[3];
333  double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL1F0, diff_ksiL0F1,
334  diff_ksiL1F1, diff_ksiL0F2, diff_ksiL1F2,
335  diff_ksiL0F3, diff_ksiL1F3};
336  int dd = 0;
337  for (; dd < 3; dd++) {
338  for (ff = 0; ff < 4; ff++) {
339  diff_ksi_faces[ff * 2 + 0][dd] =
340  (diffN[faces_nodes[3 * ff + 1] * 3 + dd] -
341  diffN[faces_nodes[3 * ff + 0] * 3 + dd]);
342  diff_ksi_faces[ff * 2 + 1][dd] =
343  (diffN[faces_nodes[3 * ff + 2] * 3 + dd] -
344  diffN[faces_nodes[3 * ff + 0] * 3 + dd]);
345  }
346  }
347  int ii = 0;
348  for (; ii < GDIM; ii++) {
349  int node_shift = ii * 4;
350  double ksi_faces[8];
351  for (ff = 0; ff < 4; ff++) {
352  ksi_faces[2 * ff + 0] = N[node_shift + faces_nodes[3 * ff + 1]] -
353  N[node_shift + faces_nodes[3 * ff + 0]];
354  ksi_faces[2 * ff + 1] = N[node_shift + faces_nodes[3 * ff + 2]] -
355  N[node_shift + faces_nodes[3 * ff + 0]];
356  }
357  int shift;
358  for (ff = 0; ff < 4; ff++) {
359  if (P[ff] == 0)
360  continue;
361  double L0[p[ff] + 1], L1[p[ff] + 1];
362  double diffL0[3 * (p[ff] + 1)], diffL1[3 * (p[ff] + 1)];
363  ierr = base_polynomials(p[ff], ksi_faces[ff * 2 + 0],
364  diff_ksi_faces[ff * 2 + 0], L0, diffL0, 3);
365  CHKERRQ(ierr);
366  ierr = base_polynomials(p[ff], ksi_faces[ff * 2 + 1],
367  diff_ksi_faces[ff * 2 + 1], L1, diffL1, 3);
368  CHKERRQ(ierr);
369  double v = N[node_shift + faces_nodes[3 * ff + 0]] *
370  N[node_shift + faces_nodes[3 * ff + 1]] *
371  N[node_shift + faces_nodes[3 * ff + 2]];
372  double v2[3] = {0, 0, 0};
373  dd = 0;
374  double n1n2 = N[node_shift + faces_nodes[3 * ff + 1]] *
375  N[node_shift + faces_nodes[3 * ff + 2]];
376  double n0n2 = N[node_shift + faces_nodes[3 * ff + 0]] *
377  N[node_shift + faces_nodes[3 * ff + 2]];
378  double n0n1 = N[node_shift + faces_nodes[3 * ff + 0]] *
379  N[node_shift + faces_nodes[3 * ff + 1]];
380  for (; dd < 3; dd++) {
381  v2[dd] = diffN[3 * faces_nodes[3 * ff + 0] + dd] * n1n2 +
382  diffN[3 * faces_nodes[3 * ff + 1] + dd] * n0n2 +
383  diffN[3 * faces_nodes[3 * ff + 2] + dd] * n0n1;
384  }
385  shift = ii * P[ff];
386  int jj = 0;
387  int oo = 0;
388  for (; oo <= (p[ff] - 3); oo++) {
389  int pp0 = 0;
390  for (; pp0 <= oo; pp0++) {
391  int pp1 = oo - pp0;
392  if (pp1 >= 0) {
393  if (faceN != NULL)
394  if (faceN[ff] != NULL) {
395  faceN[ff][shift + jj] = v * L0[pp0] * L1[pp1];
396  }
397  if (diff_faceN != NULL)
398  if (diff_faceN[ff] != NULL) {
399  dd = 0;
400  double L0L1 = L0[pp0] * L1[pp1];
401  for (; dd < 3; dd++) {
402  diff_faceN[ff][3 * shift + 3 * jj + dd] =
403  (L0[pp0] * diffL1[dd * (p[ff] + 1) + pp1] +
404  diffL0[dd * (p[ff] + 1) + pp0] * L1[pp1]) *
405  v;
406  diff_faceN[ff][3 * shift + 3 * jj + dd] += L0L1 * v2[dd];
407  }
408  }
409  jj++;
410  }
411  }
412  }
413  if (jj != P[ff])
414  SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P[ff]);
415  }
416  }
418 }

◆ H1_FaceShapeFunctions_MBTRI()

PetscErrorCode H1_FaceShapeFunctions_MBTRI ( const int *  face_nodes,
int  p,
double N,
double diffN,
double faceN,
double diff_faceN,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)

Definition at line 162 of file h1.c.

167  {
169 
170  int P = NBFACETRI_H1(p);
171  if (P == 0)
173  double diff_ksiL0[2], diff_ksiL1[2];
174  double *diff_ksi_faces[] = {diff_ksiL0, diff_ksiL1};
175  int dd = 0;
176  if (diff_faceN != NULL) {
177  for (; dd < 2; dd++) {
178  diff_ksi_faces[0][dd] =
179  diffN[face_nodes[1] * 2 + dd] - diffN[face_nodes[0] * 2 + dd];
180  diff_ksi_faces[1][dd] =
181  diffN[face_nodes[2] * 2 + dd] - diffN[face_nodes[0] * 2 + dd];
182  }
183  }
184  int ii = 0;
185  for (; ii < GDIM; ii++) {
186  int node_shift = ii * 3;
187  double ksi_faces[2];
188  ksi_faces[0] =
189  N[node_shift + face_nodes[1]] - N[node_shift + face_nodes[0]];
190  ksi_faces[1] =
191  N[node_shift + face_nodes[2]] - N[node_shift + face_nodes[0]];
192  double L0[p + 1], L1[p + 1];
193  double diffL0[2 * (p + 1)], diffL1[2 * (p + 1)];
194  ierr = base_polynomials(p, ksi_faces[0], diff_ksi_faces[0], L0, diffL0, 2);
195  CHKERRQ(ierr);
196  ierr = base_polynomials(p, ksi_faces[1], diff_ksi_faces[1], L1, diffL1, 2);
197  CHKERRQ(ierr);
198  double v = N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[1]] *
199  N[node_shift + face_nodes[2]];
200  double v2[2] = {0, 0};
201  if (diff_faceN != NULL) {
202  double n1n2 =
203  N[node_shift + face_nodes[1]] * N[node_shift + face_nodes[2]];
204  double n0n2 =
205  N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[2]];
206  double n0n1 =
207  N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[1]];
208  dd = 0;
209  for (; dd < 2; dd++) {
210  v2[dd] = diffN[face_nodes[0] * 2 + dd] * n1n2 +
211  diffN[face_nodes[1] * 2 + dd] * n0n2 +
212  diffN[face_nodes[2] * 2 + dd] * n0n1;
213  }
214  }
215  int shift = ii * P;
216  int jj = 0;
217  int oo = 0;
218  for (; oo <= (p - 3); oo++) {
219  int pp0 = 0;
220  for (; pp0 <= oo; pp0++) {
221  int pp1 = oo - pp0;
222  if (pp1 >= 0) {
223  if (faceN != NULL) {
224  faceN[shift + jj] = v * L0[pp0] * L1[pp1];
225  }
226  if (diff_faceN != NULL) {
227  dd = 0;
228  for (; dd < 2; dd++) {
229  double *val = &diff_faceN[2 * shift + 2 * jj + dd];
230  *val = (L0[pp0] * diffL1[dd * (p + 1) + pp1] +
231  diffL0[dd * (p + 1) + pp0] * L1[pp1]) *
232  v;
233  *val += L0[pp0] * L1[pp1] * v2[dd];
234  }
235  }
236  jj++;
237  }
238  }
239  }
240  if (jj != P)
241  SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P);
242  }
244 }

◆ H1_QuadShapeFunctions_MBPRISM()

PetscErrorCode H1_QuadShapeFunctions_MBPRISM ( int *  faces_nodes,
int *  p,
double N,
double diffN,
double faceN[],
double diff_faceN[],
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)

Definition at line 586 of file h1.c.

591  {
593  // TODO: return separately components of the tensor product between two edges
594 
595  int P[3];
596  int ff = 0;
597  for (; ff < 3; ff++)
598  P[ff] = NBFACEQUAD_H1(p[ff]);
599 
600  double ksi_faces[6];
601  double diff_ksiL0F0[3], diff_ksiL3F0[3];
602  double diff_ksiL0F1[3], diff_ksiL3F1[3];
603  double diff_ksiL0F2[3], diff_ksiL3F2[3];
604  double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL3F0, diff_ksiL0F1,
605  diff_ksiL3F1, diff_ksiL0F2, diff_ksiL3F2};
606  int ii = 0;
607  for (; ii < GDIM; ii++) {
608  int node_shift = ii * 6;
609  int node_diff_shift = 3 * node_shift;
610  int ff = 0;
611  for (; ff < 3; ff++) {
612  if (P[ff] == 0)
613  continue;
614  int n0 = faces_nodes[4 * ff + 0];
615  int n1 = faces_nodes[4 * ff + 1];
616  int n2 = faces_nodes[4 * ff + 2];
617  int n3 = faces_nodes[4 * ff + 3];
618  int e0 = 2 * ff + 0;
619  int e1 = 2 * ff + 1;
620  double ksi0 = N[node_shift + n0] + N[node_shift + n3];
621  double ksi1 = N[node_shift + n1] + N[node_shift + n2];
622  double eta0 = N[node_shift + n0] + N[node_shift + n1];
623  double eta1 = N[node_shift + n2] + N[node_shift + n3];
624  ksi_faces[e0] = ksi1 - ksi0;
625  ksi_faces[e1] = eta1 - eta0;
626 
627  int dd = 0;
628  for (; dd < 3; dd++) {
629  double diff_ksi0 = diffN[node_diff_shift + 3 * n0 + dd] +
630  diffN[node_diff_shift + 3 * n3 + dd];
631  double diff_ksi1 = diffN[node_diff_shift + 3 * n1 + dd] +
632  diffN[node_diff_shift + 3 * n2 + dd];
633  double diff_eta0 = diffN[node_diff_shift + 3 * n0 + dd] +
634  diffN[node_diff_shift + 3 * n1 + dd];
635  double diff_eta1 = diffN[node_diff_shift + 3 * n2 + dd] +
636  diffN[node_diff_shift + 3 * n3 + dd];
637  diff_ksi_faces[e0][dd] = diff_ksi1 - diff_ksi0;
638  diff_ksi_faces[e1][dd] = diff_eta1 - diff_eta0;
639  }
640  double L0[p[ff] + 1], L1[p[ff] + 1];
641  double diffL0[3 * (p[ff] + 1)], diffL1[3 * (p[ff] + 1)];
642  ierr = base_polynomials(p[ff], ksi_faces[e0], diff_ksi_faces[e0], L0,
643  diffL0, 3);
644  CHKERRQ(ierr);
645  ierr = base_polynomials(p[ff], ksi_faces[e1], diff_ksi_faces[e1], L1,
646  diffL1, 3);
647  CHKERRQ(ierr);
648 
649  double v = N[node_shift + n0] * N[node_shift + n2] +
650  N[node_shift + n1] * N[node_shift + n3];
651 
652  double diff_v[3];
653  dd = 0;
654  for (; dd < 3; ++dd)
655  diff_v[dd] =
656 
657  diffN[node_diff_shift + 3 * n0 + dd] * N[node_shift + n2] +
658  N[node_shift + n0] * diffN[node_diff_shift + 3 * n2 + dd] +
659 
660  diffN[node_diff_shift + 3 * n1 + dd] * N[node_shift + n3] +
661  N[node_shift + n1] * diffN[node_diff_shift + 3 * n3 + dd];
662 
663  int shift;
664  shift = ii * P[ff];
665 
666  if (faceN != NULL) {
667  if (faceN[ff] != NULL) {
668  int jj = 0;
669  int oo = 0;
670  for (; oo <= p[ff] - 2; ++oo) {
671  int pp0 = 0;
672  for (; pp0 < oo; ++pp0) {
673  int pp1 = oo;
674  faceN[ff][shift + jj] = L0[pp0] * L1[pp1] * v;
675  ++jj;
676  faceN[ff][shift + jj] = L0[pp1] * L1[pp0] * v;
677  ++jj;
678  }
679  faceN[ff][shift + jj] = L0[oo] * L1[oo] * v;
680  ++jj;
681  }
682  if (jj != P[ff])
683  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
684  "Inconsistent implementation (bug in the code) %d != %d",
685  jj, P);
686  }
687  }
688 
689  if (diff_faceN != NULL) {
690  if (diff_faceN[ff] != NULL) {
691  int jj = 0;
692  int oo = 0;
693  for (; oo <= p[ff] - 2; ++oo) {
694  int pp0 = 0;
695  for (; pp0 < oo; ++pp0) {
696  int pp1 = oo;
697  int dd;
698  for (dd = 0; dd < 3; dd++) {
699  diff_faceN[ff][3 * shift + 3 * jj + dd] =
700  (L0[pp0] * diffL1[dd * (p[ff] + 1) + pp1] +
701  diffL0[dd * (p[ff] + 1) + pp0] * L1[pp1]) *
702  v +
703  L0[pp0] * L1[pp1] * diff_v[dd];
704  }
705  ++jj;
706  for (dd = 0; dd < 3; dd++) {
707  diff_faceN[ff][3 * shift + 3 * jj + dd] =
708  (L0[pp1] * diffL1[dd * (p[ff] + 1) + pp0] +
709  diffL0[dd * (p[ff] + 1) + pp1] * L1[pp0]) *
710  v +
711  L0[pp1] * L1[pp0] * diff_v[dd];
712  }
713  ++jj;
714  }
715  for (dd = 0; dd < 3; dd++) {
716  diff_faceN[ff][3 * shift + 3 * jj + dd] =
717  (L0[oo] * diffL1[dd * (p[ff] + 1) + oo] +
718  diffL0[dd * (p[ff] + 1) + oo] * L1[oo]) *
719  v +
720  L0[oo] * L1[oo] * diff_v[dd];
721  }
722  ++jj;
723  }
724  if (jj != P[ff])
725  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
726  "Inconsistent implementation (bug in the code) %d != %d",
727  jj, P);
728  }
729  }
730  }
731  }
733 }

◆ H1_QuadShapeFunctions_MBQUAD()

PetscErrorCode H1_QuadShapeFunctions_MBQUAD ( int *  faces_nodes,
int  p,
double N,
double diffN,
double faceN,
double diff_faceN,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)
Examples
edge_and_bubble_shape_functions_on_quad.cpp.

Definition at line 903 of file h1.c.

908  {
910  // TODO: return separately components of the tensor product between two edges
911 
912  int P = NBFACEQUAD_H1(p);
913  if (P == 0)
915  double diff_ksiL0F0[2], diff_ksiL1F0[2];
916  double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL1F0};
917  int ii = 0;
918  for (; ii < GDIM; ii++) {
919  int node_shift = ii * 4;
920  int node_diff_shift = 2 * node_shift;
921 
922  int n0 = faces_nodes[0];
923  int n1 = faces_nodes[1];
924  int n2 = faces_nodes[2];
925  int n3 = faces_nodes[3];
926  double ksi0 = N[node_shift + n0] + N[node_shift + n3];
927  double ksi1 = N[node_shift + n1] + N[node_shift + n2];
928  double eta0 = N[node_shift + n0] + N[node_shift + n1];
929  double eta1 = N[node_shift + n2] + N[node_shift + n3];
930  double ksi_faces = ksi1 - ksi0;
931  double eta_faces = eta1 - eta0;
932  double diff_ksi_faces[2];
933  double diff_eta_faces[2];
934  int dd = 0;
935  for (; dd < 2; dd++) {
936  double diff_ksi0 = diffN[node_diff_shift + 2 * n0 + dd] +
937  diffN[node_diff_shift + 2 * n3 + dd];
938  double diff_ksi1 = diffN[node_diff_shift + 2 * n1 + dd] +
939  diffN[node_diff_shift + 2 * n2 + dd];
940  double diff_eta0 = diffN[node_diff_shift + 2 * n0 + dd] +
941  diffN[node_diff_shift + 2 * n1 + dd];
942  double diff_eta1 = diffN[node_diff_shift + 2 * n2 + dd] +
943  diffN[node_diff_shift + 2 * n3 + dd];
944 
945  diff_ksi_faces[dd] = diff_ksi1 - diff_ksi0;
946  diff_eta_faces[dd] = diff_eta1 - diff_eta0;
947  }
948  double L0[p + 1], L1[p + 1];
949  double diffL0[2 * (p + 1)], diffL1[2 * (p + 1)];
950  ierr = base_polynomials(p, ksi_faces, diff_ksi_faces, L0, diffL0, 2);
951  CHKERRQ(ierr);
952  ierr = base_polynomials(p, eta_faces, diff_eta_faces, L1, diffL1, 2);
953  CHKERRQ(ierr);
954 
955  double v = N[node_shift + n0] * N[node_shift + n2] +
956  N[node_shift + n1] * N[node_shift + n3];
957 
958  double diff_v[2];
959  dd = 0;
960  for (; dd < 2; ++dd)
961  diff_v[dd] =
962 
963  diffN[node_diff_shift + 2 * n0 + dd] * N[node_shift + n2] +
964  N[node_shift + n0] * diffN[node_diff_shift + 2 * n2 + dd] +
965 
966  diffN[node_diff_shift + 2 * n1 + dd] * N[node_shift + n3] +
967  N[node_shift + n1] * diffN[node_diff_shift + 2 * n3 + dd];
968 
969  const int shift = ii * P;
970 
971  if (faceN != NULL) {
972  int jj = 0;
973  int oo = 0;
974  for (; oo <= p - 2; ++oo) {
975  int pp0 = 0;
976  for (; pp0 < oo; ++pp0) {
977  int pp1 = oo;
978  faceN[shift + jj] = L0[pp0] * L1[pp1] * v;
979  ++jj;
980  faceN[shift + jj] = L0[pp1] * L1[pp0] * v;
981  ++jj;
982  }
983  faceN[shift + jj] = L0[oo] * L1[oo] * v;
984  ++jj;
985  }
986  if (jj != P)
987  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
988  "Inconsistent implementation (bug in the code) %d != %d", jj,
989  P);
990  }
991 
992  if (diff_faceN != NULL) {
993  int jj = 0;
994  int oo = 0;
995  for (; oo <= p - 2; ++oo) {
996  int pp0 = 0;
997  for (; pp0 < oo; ++pp0) {
998  int pp1 = oo;
999  int dd;
1000  for (dd = 0; dd < 2; dd++) {
1001  diff_faceN[2 * shift + 2 * jj + dd] =
1002  (L0[pp0] * diffL1[dd * (p + 1) + pp1] +
1003  diffL0[dd * (p + 1) + pp0] * L1[pp1]) *
1004  v +
1005  L0[pp0] * L1[pp1] * diff_v[dd];
1006  }
1007  ++jj;
1008  for (dd = 0; dd < 2; dd++) {
1009  diff_faceN[2 * shift + 2 * jj + dd] =
1010  (L0[pp1] * diffL1[dd * (p + 1) + pp0] +
1011  diffL0[dd * (p + 1) + pp1] * L1[pp0]) *
1012  v +
1013  L0[pp1] * L1[pp0] * diff_v[dd];
1014  }
1015  ++jj;
1016  }
1017  for (dd = 0; dd < 2; dd++) {
1018  diff_faceN[2 * shift + 2 * jj + dd] =
1019  (L0[oo] * diffL1[dd * (p + 1) + oo] +
1020  diffL0[dd * (p + 1) + oo] * L1[oo]) *
1021  v +
1022  L0[oo] * L1[oo] * diff_v[dd];
1023  }
1024  ++jj;
1025  }
1026  if (jj != P)
1027  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1028  "Inconsistent implementation (bug in the code) %d != %d", jj,
1029  P);
1030  }
1031  }
1033 }

◆ H1_VolumeGradientOfDeformation_hierachical()

DEPRECATED PetscErrorCode H1_VolumeGradientOfDeformation_hierachical ( int  p,
double diffN,
double dofs,
double F 
)
Deprecated:
use H1_VolumeGradientOfDeformation_hierarchical

Definition at line 1230 of file h1.c.

1232  {
1233  return H1_VolumeGradientOfDeformation_hierarchical(p, diffN, dofs, F);
1234 }

◆ H1_VolumeGradientOfDeformation_hierarchical()

PetscErrorCode H1_VolumeGradientOfDeformation_hierarchical ( int  p,
double diffN,
double dofs,
double F 
)

Definition at line 575 of file h1.c.

577  {
579  int col, row = 0;
580  for (; row < 3; row++)
581  for (col = 0; col < 3; col++)
582  F[3 * row + col] =
583  cblas_ddot(NBVOLUMETET_H1(p), &diffN[col], 3, &dofs[row], 3);
585 }

◆ H1_VolumeShapeDiffMBTETinvJ()

PetscErrorCode H1_VolumeShapeDiffMBTETinvJ ( int  base_p,
int  p,
double volume_diffN,
double invJac,
double volume_diffNinvJac,
int  GDIM 
)

Definition at line 536 of file h1.c.

539  {
541  int ii, gg;
542  for (ii = 0; ii < NBVOLUMETET_H1(p); ii++) {
543  for (gg = 0; gg < GDIM; gg++) {
544  int shift1 = NBVOLUMETET_H1(base_p) * gg;
545  int shift2 = NBVOLUMETET_H1(p) * gg;
547  &(volume_diffN)[3 * shift1 + 3 * ii], 1, 0.,
548  &(volume_diffNinvJac)[3 * shift2 + 3 * ii], 1);
549  }
550  }
552 }

◆ H1_VolumeShapeFunctions_MBPRISM()

PetscErrorCode H1_VolumeShapeFunctions_MBPRISM ( int  p,
double N,
double diffN,
double volumeN,
double diff_volumeN,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)

Definition at line 734 of file h1.c.

739  {
741 
742  int P = NBVOLUMEPRISM_H1(p);
743  if (P == 0)
745  double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
746  int ii = 0;
747  for (; ii < GDIM; ii++) {
748  int node_shift = ii * 6;
749  int node_diff_shift = ii * 18;
750 
751  double n0 = N[node_shift + 0];
752  double n1 = N[node_shift + 1];
753  double n2 = N[node_shift + 2];
754  double n3 = N[node_shift + 3];
755  double n4 = N[node_shift + 4];
756  double n5 = N[node_shift + 5];
757 
758  double ksiL0 = n1 + n4 - n0 - n3;
759  double ksiL1 = n2 + n5 - n0 - n3;
760  double ksiL2 = (n3 + n4 + n5) - (n0 + n1 + n2);
761 
762  int dd = 0;
763  for (; dd < 3; dd++) {
764  double diff_n0 = diffN[node_diff_shift + 3 * 0 + dd];
765  double diff_n1 = diffN[node_diff_shift + 3 * 1 + dd];
766  double diff_n2 = diffN[node_diff_shift + 3 * 2 + dd];
767  double diff_n3 = diffN[node_diff_shift + 3 * 3 + dd];
768  double diff_n4 = diffN[node_diff_shift + 3 * 4 + dd];
769  double diff_n5 = diffN[node_diff_shift + 3 * 5 + dd];
770  diff_ksiL0[dd] = (diff_n1 + diff_n4) - (diff_n0 + diff_n3);
771  diff_ksiL1[dd] = (diff_n2 + diff_n5) - (diff_n0 + diff_n3);
772  diff_ksiL2[dd] =
773  (diff_n3 + diff_n4 + diff_n5) - (diff_n0 + diff_n1 + diff_n2);
774  }
775 
776  double L0[p + 1], L1[p + 1], L2[p + 1];
777  double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
778  ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
779  CHKERRQ(ierr);
780  ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
781  CHKERRQ(ierr);
782  ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
783  CHKERRQ(ierr);
784 
785  double v_tri = (n0 + n3) * (n1 + n4) * (n2 + n5);
786  double v_edge = (n0 + n1 + n2) * (n3 + n4 + n5);
787  double v = v_tri * v_edge;
788 
789  double diff_v_tri[3];
790  double diff_v_edge[3];
791  double diff_v[3];
792  dd = 0;
793  for (; dd < 3; dd++) {
794  double diff_n0 = diffN[node_diff_shift + 3 * 0 + dd];
795  double diff_n1 = diffN[node_diff_shift + 3 * 1 + dd];
796  double diff_n2 = diffN[node_diff_shift + 3 * 2 + dd];
797  double diff_n3 = diffN[node_diff_shift + 3 * 3 + dd];
798  double diff_n4 = diffN[node_diff_shift + 3 * 4 + dd];
799  double diff_n5 = diffN[node_diff_shift + 3 * 5 + dd];
800  diff_v_tri[dd] = (diff_n0 + diff_n3) * (n1 + n4) * (n2 + n5) +
801  (diff_n1 + diff_n4) * (n0 + n3) * (n2 + n5) +
802  (diff_n2 + diff_n5) * (n0 + n3) * (n1 + n4);
803  diff_v_edge[dd] = (diff_n0 + diff_n1 + diff_n2) * (n3 + n4 + n5) +
804  (diff_n3 + diff_n4 + diff_n5) * (n0 + n1 + n2);
805  diff_v[dd] = diff_v_tri[dd] * v_edge + v_tri * diff_v_edge[dd];
806  }
807 
808  int shift = ii * P;
809 
810  if (volumeN != NULL) {
811  int jj = 0;
812  int oo, pp0, pp1, zz;
813  for (oo = 0; oo <= p - 3; ++oo) {
814  for (pp0 = 0; pp0 < oo; ++pp0) {
815  for (pp1 = 0; pp1 < oo; ++pp1) {
816  const int perm[3][3] = {
817  {oo, pp0, pp1}, {pp0, oo, pp1}, {pp0, pp1, oo}};
818  for (zz = 0; zz != 3; ++zz) {
819  volumeN[shift + jj] =
820  L0[perm[zz][0]] * L1[perm[zz][1]] * L2[perm[zz][2]] * v;
821  ++jj;
822  }
823  }
824  const int perm[3][3] = {{pp0, oo, oo}, {oo, pp0, oo}, {oo, oo, pp0}};
825  for (zz = 0; zz != 3; ++zz) {
826  volumeN[shift + jj] =
827  L0[perm[zz][0]] * L1[perm[zz][1]] * L2[perm[zz][2]] * v;
828  ++jj;
829  }
830  }
831  volumeN[shift + jj] = L0[oo] * L1[oo] * L2[oo] * v;
832  ++jj;
833  }
834  if (jj != P)
835  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
836  "Inconsistent implementation (bug in the code) %d != %d", jj,
837  P);
838  }
839 
840  if (diff_volumeN != NULL) {
841  int jj = 0;
842  int oo, pp0, pp1, zz, dd;
843  for (oo = 0; oo <= p - 3; ++oo) {
844  for (pp0 = 0; pp0 < oo; ++pp0) {
845  for (pp1 = 0; pp1 < oo; ++pp1) {
846  const int perm[3][3] = {
847  {oo, pp0, pp1}, {pp0, oo, pp1}, {pp0, pp1, oo}};
848  for (zz = 0; zz != 3; ++zz) {
849  const int i = perm[zz][0];
850  const int j = perm[zz][1];
851  const int k = perm[zz][2];
852  for (dd = 0; dd != 3; ++dd) {
853  diff_volumeN[3 * shift + 3 * jj + dd] =
854  (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
855  L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
856  L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
857  v +
858  L0[i] * L1[j] * L2[k] * diff_v[dd];
859  }
860  ++jj;
861  }
862  }
863  const int perm[3][3] = {{pp0, oo, oo}, {oo, pp0, oo}, {oo, oo, pp0}};
864  for (zz = 0; zz != 3; ++zz) {
865  const int i = perm[zz][0];
866  const int j = perm[zz][1];
867  const int k = perm[zz][2];
868  for (dd = 0; dd != 3; ++dd) {
869  diff_volumeN[3 * shift + 3 * jj + dd] =
870  (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
871  L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
872  L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
873  v +
874  L0[i] * L1[j] * L2[k] * diff_v[dd];
875  }
876  ++jj;
877  }
878  }
879 
880  const int i = oo;
881  const int j = oo;
882  const int k = oo;
883  for (dd = 0; dd != 3; ++dd) {
884  diff_volumeN[3 * shift + 3 * jj + dd] =
885  (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
886  L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
887  L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
888  v +
889  L0[i] * L1[j] * L2[k] * diff_v[dd];
890  }
891 
892  ++jj;
893  }
894  if (jj != P)
895  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
896  "Inconsistent implementation (bug in the code) %d != %d", jj,
897  P);
898  }
899  }
901 }

◆ H1_VolumeShapeFunctions_MBTET()

PetscErrorCode H1_VolumeShapeFunctions_MBTET ( int  p,
double N,
double diffN,
double volumeN,
double diff_volumeN,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)

Definition at line 419 of file h1.c.

424  {
426 
427  int P = NBVOLUMETET_H1(p);
428  if (P == 0)
430  double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
431  int dd = 0;
432  for (; dd < 3; dd++) {
433  diff_ksiL0[dd] = (diffN[1 * 3 + dd] - diffN[0 * 3 + dd]);
434  diff_ksiL1[dd] = (diffN[2 * 3 + dd] - diffN[0 * 3 + dd]);
435  diff_ksiL2[dd] = (diffN[3 * 3 + dd] - diffN[0 * 3 + dd]);
436  }
437  int ii = 0;
438  for (; ii < GDIM; ii++) {
439  int node_shift = ii * 4;
440  double ksiL0 = N[node_shift + 1] - N[node_shift + 0];
441  double ksiL1 = N[node_shift + 2] - N[node_shift + 0];
442  double ksiL2 = N[node_shift + 3] - N[node_shift + 0];
443  double L0[p + 1], L1[p + 1], L2[p + 1];
444  double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
445  ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
446  CHKERRQ(ierr);
447  ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
448  CHKERRQ(ierr);
449  ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
450  CHKERRQ(ierr);
451  double v = N[node_shift + 0] * N[node_shift + 1] * N[node_shift + 2] *
452  N[node_shift + 3];
453  double v2[3] = {0, 0, 0};
454  dd = 0;
455  for (; dd < 3; dd++) {
456  v2[dd] = diffN[3 * 0 + dd] * N[node_shift + 1] * N[node_shift + 2] *
457  N[node_shift + 3] +
458  N[node_shift + 0] * diffN[3 * 1 + dd] * N[node_shift + 2] *
459  N[node_shift + 3] +
460  N[node_shift + 0] * N[node_shift + 1] * diffN[3 * 2 + dd] *
461  N[node_shift + 3] +
462  N[node_shift + 0] * N[node_shift + 1] * N[node_shift + 2] *
463  diffN[3 * 3 + dd];
464  }
465  int shift = ii * P;
466  int jj = 0;
467  int oo = 0;
468  for (; oo <= (p - 4); oo++) {
469  int pp0 = 0;
470  for (; pp0 <= oo; pp0++) {
471  int pp1 = 0;
472  for (; (pp0 + pp1) <= oo; pp1++) {
473  int pp2 = oo - pp0 - pp1;
474  if (pp2 >= 0) {
475  if (volumeN != NULL) {
476  volumeN[shift + jj] = L0[pp0] * L1[pp1] * L2[pp2] * v;
477  }
478  if (diff_volumeN != NULL) {
479  dd = 0;
480  for (; dd < 3; dd++) {
481  diff_volumeN[3 * shift + 3 * jj + dd] =
482  (diffL0[dd * (p + 1) + pp0] * L1[pp1] * L2[pp2] +
483  L0[pp0] * diffL1[dd * (p + 1) + pp1] * L2[pp2] +
484  L0[pp0] * L1[pp1] * diffL2[dd * (p + 1) + pp2]) *
485  v;
486  diff_volumeN[3 * shift + 3 * jj + dd] +=
487  L0[pp0] * L1[pp1] * L2[pp2] * v2[dd];
488  }
489  }
490  jj++;
491  }
492  }
493  }
494  }
495  if (jj != P)
496  SETERRQ1(PETSC_COMM_SELF, 1, "wrong order %d", jj);
497  }
499 }

◆ L2_Ainsworth_ShapeFunctions_MBTET()

PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET ( int  p,
double N,
double diffN,
double L2N,
double diff_L2N,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)

Get base functions on tetrahedron for L2 space.

Parameters
ppolynomial order
Nbarycentric coordinates (shape functions) at integration points
diffNdirevatives of barycentric coordinates, i.e. direvatives of shape functions
L2Nvalues of L2 base at integration points
diff_L2Ndirvatives of base functions at integration points
GDIMnumber of integration points
base_polynomialspolynomial base used to construct L2 base on element
Returns
PetscErrorCode

Definition at line 84 of file l2.c.

88  {
90 
91  int P = NBVOLUMETET_L2(p);
92  if (P == 0)
94  double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
95  int dd = 0;
96  if (diffN != NULL) {
97  for (; dd < 3; dd++) {
98  diff_ksiL0[dd] = (diffN[1 * 3 + dd] - diffN[0 * 3 + dd]);
99  diff_ksiL1[dd] = (diffN[2 * 3 + dd] - diffN[0 * 3 + dd]);
100  diff_ksiL2[dd] = (diffN[3 * 3 + dd] - diffN[0 * 3 + dd]);
101  }
102  }
103  int ii = 0;
104  for (; ii != GDIM; ++ii) {
105  int node_shift = ii * 4;
106  double ksiL0 = N[node_shift + 1] - N[node_shift + 0];
107  double ksiL1 = N[node_shift + 2] - N[node_shift + 0];
108  double ksiL2 = N[node_shift + 3] - N[node_shift + 0];
109  double L0[p + 1], L1[p + 1], L2[p + 1];
110  double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
111  if (diffN != NULL) {
112  ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
113  CHKERRQ(ierr);
114  ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
115  CHKERRQ(ierr);
116  ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
117  CHKERRQ(ierr);
118  } else {
119  ierr = base_polynomials(p, ksiL0, NULL, L0, NULL, 3);
120  CHKERRQ(ierr);
121  ierr = base_polynomials(p, ksiL1, NULL, L1, NULL, 3);
122  CHKERRQ(ierr);
123  ierr = base_polynomials(p, ksiL2, NULL, L2, NULL, 3);
124  CHKERRQ(ierr);
125  }
126  int shift = ii * P;
127  int jj = 0;
128  int oo = 0;
129  for (; oo <= p; oo++) {
130  int pp0 = 0;
131  for (; pp0 <= oo; pp0++) {
132  int pp1 = 0;
133  for (; (pp0 + pp1) <= oo; pp1++) {
134  int pp2 = oo - pp0 - pp1;
135  if (pp2 >= 0) {
136  if (L2N != NULL) {
137  L2N[shift + jj] = L0[pp0] * L1[pp1] * L2[pp2];
138  }
139  if (diff_L2N != NULL) {
140  int dd = 0;
141  for (; dd < 3; dd++) {
142  diff_L2N[3 * shift + 3 * jj + dd] =
143  diffL0[dd * (p + 1) + pp0] * L1[pp1] * L2[pp2] +
144  L0[pp0] * diffL1[dd * (p + 1) + pp1] * L2[pp2] +
145  L0[pp0] * L1[pp1] * diffL2[dd * (p + 1) + pp2];
146  }
147  }
148  jj++;
149  }
150  }
151  }
152  }
153  if (jj != P)
154  SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P);
155  }
157 }

◆ L2_Ainsworth_ShapeFunctions_MBTRI()

PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTRI ( int  p,
double N,
double diffN,
double L2N,
double diff_L2N,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)

Get base functions on triangle for L2 space.

Parameters
ppolynomial order
Nbarycentric coordinates (shape functions) at integration points
diffNdirevatives of barycentric coordinates, i.e. direvatives of shape functions
L2Nvalues of L2 base at integration points
diff_L2Ndirvatives of base functions at integration points
GDIMnumber of integration points
base_polynomialspolynomial base used to construct L2 base on element
Returns
PetscErrorCode

Definition at line 29 of file l2.c.

33  {
35 
36  int P = NBFACETRI_L2(p);
37  if (P == 0)
39  double diff_ksiL01[2], diff_ksiL20[2];
40  int dd = 0;
41  for (; dd < 2; dd++) {
42  diff_ksiL01[dd] = (diffN[1 * 2 + dd] - diffN[0 * 2 + dd]);
43  diff_ksiL20[dd] = (diffN[2 * 2 + dd] - diffN[0 * 2 + dd]);
44  }
45  int ii = 0;
46  for (; ii != GDIM; ++ii) {
47  int node_shift = ii * 3;
48  double ksiL01 = N[node_shift + 1] - N[node_shift + 0];
49  double ksiL20 = N[node_shift + 2] - N[node_shift + 0];
50  double L01[p + 1], L20[p + 1];
51  double diffL01[2 * (p + 1)], diffL20[2 * (p + 1)];
52  ierr = base_polynomials(p, ksiL01, diff_ksiL01, L01, diffL01, 2);
53  CHKERRQ(ierr);
54  ierr = base_polynomials(p, ksiL20, diff_ksiL20, L20, diffL20, 2);
55  CHKERRQ(ierr);
56  int shift = ii * P;
57  int jj = 0;
58  int oo = 0;
59  for (; oo <= p; oo++) {
60  int pp0 = 0;
61  for (; pp0 <= oo; pp0++) {
62  int pp1 = oo - pp0;
63  if (pp1 >= 0) {
64  if (L2N != NULL) {
65  L2N[shift + jj] = L01[pp0] * L20[pp1];
66  }
67  if (diff_L2N != NULL) {
68  int dd = 0;
69  for (; dd < 2; dd++) {
70  diff_L2N[2 * shift + 2 * jj + dd] =
71  diffL01[dd * (p + 1) + pp0] * L20[pp1] +
72  L01[pp0] * diffL20[dd * (p + 1) + pp1];
73  }
74  }
75  jj++;
76  }
77  }
78  }
79  if (jj != P)
80  SETERRQ1(PETSC_COMM_SELF, 1, "wrong order %d", jj);
81  }
83 }

◆ L2_ShapeFunctions_MBTET()

DEPRECATED PetscErrorCode L2_ShapeFunctions_MBTET ( int  p,
double N,
double diffN,
double L2N,
double diff_L2N,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)
Deprecated:
Use L2_Ainsworth_ShapeFunctions_MBTET

Definition at line 185 of file l2.c.

189  {
190  return L2_Ainsworth_ShapeFunctions_MBTET(p, N, diffN, L2N, diff_L2N, GDIM,
191  base_polynomials);
192 }

◆ L2_ShapeFunctions_MBTRI()

DEPRECATED PetscErrorCode L2_ShapeFunctions_MBTRI ( int  p,
double N,
double diffN,
double L2N,
double diff_L2N,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim base_polynomials 
)
Deprecated:
Use L2_Ainsworth_ShapeFunctions_MBTRI

Definition at line 176 of file l2.c.

180  {
181  return L2_Ainsworth_ShapeFunctions_MBTRI(p, N, diffN, L2N, diff_L2N, GDIM,
182  base_polynomials);
183 }

◆ L2_VolumeShapeDiffMBTETinvJ()

PetscErrorCode L2_VolumeShapeDiffMBTETinvJ ( int  base_p,
int  p,
double volume_diffN,
double invJac,
double volume_diffNinvJac,
int  GDIM 
)

Definition at line 158 of file l2.c.

161  {
163  int ii, gg;
164  for (ii = 0; ii != NBVOLUMETET_L2(p); ++ii) {
165  for (gg = 0; gg < GDIM; gg++) {
166  int shift1 = NBVOLUMETET_L2(base_p) * gg;
167  int shift2 = NBVOLUMETET_L2(p) * gg;
169  &(volume_diffN)[3 * shift1 + 3 * ii], 1, 0.,
170  &(volume_diffNinvJac)[3 * shift2 + 3 * ii], 1);
171  }
172  }
174 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
cblas_daxpy
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
ierr
static PetscErrorCode ierr
Definition: l2.c:27
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
L
static Index< 'L', 3 > L
Definition: BasicFeTools.hpp:85
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
L2
@ L2
field with C-1 continuity
Definition: definitions.h:180
cblas_dgemv
void cblas_dgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const double alpha, const double *A, const int lda, const double *X, const int incX, const double beta, double *Y, const int incY)
Definition: cblas_dgemv.c:11
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
H1_EdgeGradientOfDeformation_hierarchical
PetscErrorCode H1_EdgeGradientOfDeformation_hierarchical(int p, double *diffN, double *dofs, double *F)
Definition: h1.c:553
NBVOLUMETET_L2
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:37
L2_Ainsworth_ShapeFunctions_MBTRI
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTRI(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on triangle for L2 space.
Definition: l2.c:29
CblasTrans
@ CblasTrans
Definition: cblas.h:11
CblasRowMajor
@ CblasRowMajor
Definition: cblas.h:10
H1_FaceGradientOfDeformation_hierarchical
PetscErrorCode H1_FaceGradientOfDeformation_hierarchical(int p, double *diffN, double *dofs, double *F)
Definition: h1.c:564
invJac
MatrixDouble invJac
Definition: HookeElement.hpp:698
NBFACEQUAD_H1
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
Definition: h1_hdiv_hcurl_l2.h:65
N
const int N
Definition: speed_test.cpp:3
ierr
static PetscErrorCode ierr
Definition: h1.c:25
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:80
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
NBFACETRI_L2
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
Definition: h1_hdiv_hcurl_l2.h:42
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
H1_VolumeGradientOfDeformation_hierarchical
PetscErrorCode H1_VolumeGradientOfDeformation_hierarchical(int p, double *diffN, double *dofs, double *F)
Definition: h1.c:575
cblas_ddot
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
NBVOLUMEPRISM_H1
#define NBVOLUMEPRISM_H1(P)
Number of base functions on prism for H1 space.
Definition: h1_hdiv_hcurl_l2.h:80
L2_Ainsworth_ShapeFunctions_MBTET
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on tetrahedron for L2 space.
Definition: l2.c:84
i
FTensor::Index< 'i', 3 > i
Definition: matrix_function.cpp:18
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
cblas_dscal
void cblas_dscal(const int N, const double alpha, double *X, const int incX)
Definition: cblas_dscal.c:11
cblas_dcopy
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33