v0.13.1
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 NBVOLUMEHEX_L2_GENERAL(P, Q, R)   ((P + 1) * (Q + 1) * (R + 1))
 Number of base functions on hexahedron for L2 space. More...
 
#define NBVOLUMEHEX_L2(P)   (NBVOLUMEHEX_L2_GENERAL(P, P, P))
 Number of base functions on hexahedron 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 NBVOLUMEHEX_H1_GENERAL(P, Q, R)    ((((P) > 1) && ((Q) > 1) && ((R) > 1)) ? (((P)-1) * ((Q)-1) * ((R)-1)) : 0)
 Number of base functions on hex for H1 space. More...
 
#define NBVOLUMEHEX_H1(P)   (NBVOLUMEHEX_H1_GENERAL(P, P, P))
 Number of base functions on hex 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_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_HCURL(P, P))
 
#define NBVOLUMEHEX_DEMKOWICZ_FAMILY_HCURL(P, Q, R)    ((P > 0) && (Q > 1) && (R > 1) ? ((P) * (Q - 1) * (R - 1)) : 0)
 
#define NBVOLUMEHEX_DEMKOWICZ_HCURL(P)    (3 * NBVOLUMEHEX_DEMKOWICZ_FAMILY_HCURL(P, 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)
 
#define NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL(P, Q)    (((P) > 0 && (Q) > 0) ? ((P) * (Q)) : 0)
 
#define NBFACEQUAD_DEMKOWICZ_HDIV(P)    (NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL(P, P))
 
#define NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV(P, Q, R)    ((((P) > 0) && ((Q) > 0) && ((R) > 0)) ? ((P - 1) * Q * R) : 0)
 
#define NBVOLUMEHEX_DEMKOWICZ_HDIV(P)    (3 * NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV(P, P, P))
 

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...
 
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)
 
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 107 of file h1_hdiv_hcurl_l2.h.

◆ NBEDGE_DEMKOWICZ_HCURL

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

Definition at line 118 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 65 of file h1_hdiv_hcurl_l2.h.

◆ NBEDGE_HDIV

#define NBEDGE_HDIV (   P)    (0)

Definition at line 139 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 58 of file h1_hdiv_hcurl_l2.h.

◆ NBFACEQUAD_DEMKOWICZ_FAMILY_HCURL

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

Number of base functions on quad for Hcurl space.

Definition at line 126 of file h1_hdiv_hcurl_l2.h.

◆ NBFACEQUAD_DEMKOWICZ_HCURL

#define NBFACEQUAD_DEMKOWICZ_HCURL (   P)     (2 * NBFACEQUAD_DEMKOWICZ_FAMILY_HCURL(P, P))

Definition at line 128 of file h1_hdiv_hcurl_l2.h.

◆ NBFACEQUAD_DEMKOWICZ_HDIV

#define NBFACEQUAD_DEMKOWICZ_HDIV (   P)     (NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL(P, P))

Definition at line 155 of file h1_hdiv_hcurl_l2.h.

◆ NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL

#define NBFACEQUAD_DEMKOWICZ_QUAD_HDIV_GEMERAL (   P,
 
)     (((P) > 0 && (Q) > 0) ? ((P) * (Q)) : 0)

Definition at line 153 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 75 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 80 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 108 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_AINSWORTH_EDGE_HDIV

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

Definition at line 140 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 109 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 141 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 110 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 142 of file h1_hdiv_hcurl_l2.h.

◆ NBFACETRI_DEMKOWICZ_HCURL

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

Definition at line 119 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 70 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 52 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMEHEX_DEMKOWICZ_FAMILY_HCURL

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

Definition at line 131 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV

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

Definition at line 157 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMEHEX_DEMKOWICZ_HCURL

#define NBVOLUMEHEX_DEMKOWICZ_HCURL (   P)     (3 * NBVOLUMEHEX_DEMKOWICZ_FAMILY_HCURL(P, P, P))

Definition at line 134 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMEHEX_DEMKOWICZ_HDIV

#define NBVOLUMEHEX_DEMKOWICZ_HDIV (   P)     (3 * NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV(P, P, P))

Definition at line 159 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMEHEX_H1

#define NBVOLUMEHEX_H1 (   P)    (NBVOLUMEHEX_H1_GENERAL(P, P, P))

Number of base functions on hex for H1 space.

Definition at line 103 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMEHEX_H1_GENERAL

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

Number of base functions on hex for H1 space.

Definition at line 96 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMEHEX_L2

#define NBVOLUMEHEX_L2 (   P)    (NBVOLUMEHEX_L2_GENERAL(P, P, P))

Number of base functions on hexahedron for L2 space.

Definition at line 47 of file h1_hdiv_hcurl_l2.h.

◆ NBVOLUMEHEX_L2_GENERAL

#define NBVOLUMEHEX_L2_GENERAL (   P,
  Q,
 
)    ((P + 1) * (Q + 1) * (R + 1))

Number of base functions on hexahedron 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 90 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 143 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 111 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 144 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 115 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 147 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 113 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 145 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 120 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 85 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_hierarchical()

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

Definition at line 619 of file h1.c.

621 {
623 int col, row = 0;
624 for (; row < 3; row++)
625 for (col = 0; col < 3; col++)
626 F[3 * row + col] =
627 cblas_ddot(NBEDGE_H1(p), &diffN[col], 3, &dofs[row], 3);
629}
static Index< 'p', 3 > p
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.

◆ H1_EdgeShapeDiffMBTETinvJ()

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

Definition at line 566 of file h1.c.

568 {
570 int ee = 0, ii, gg;
571 for (; ee < 6; ee++) {
572 for (ii = 0; ii < NBEDGE_H1(p[ee]); ii++) {
573 for (gg = 0; gg < GDIM; gg++) {
574 int shift1 = NBEDGE_H1(base_p[ee]) * gg;
575 int shift2 = NBEDGE_H1(p[ee]) * gg;
576 cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
577 &(edge_diffN[ee])[3 * shift1 + 3 * ii], 1, 0.,
578 &(edge_diffNinvJac[ee])[3 * shift2 + 3 * ii], 1);
579 }
580 }
581 }
583}
MatrixDouble invJac

◆ 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 1101 of file h1.c.

1106 {
1108
1109 double *edgeN01 = NULL, *edgeN12 = NULL, *edgeN23 = NULL, *edgeN30 = NULL;
1110 if (edgeN != NULL) {
1111 edgeN01 = edgeN[0];
1112 edgeN12 = edgeN[1];
1113 edgeN23 = edgeN[2];
1114 edgeN30 = edgeN[3];
1115 }
1116 double *diff_edgeN01 = NULL, *diff_edgeN12 = NULL, *diff_edgeN23 = NULL,
1117 *diff_edgeN30 = NULL;
1118 if (diff_edgeN != NULL) {
1119 diff_edgeN01 = diff_edgeN[0];
1120 diff_edgeN12 = diff_edgeN[1];
1121 diff_edgeN23 = diff_edgeN[2];
1122 diff_edgeN30 = diff_edgeN[3];
1123 }
1124 int P[4];
1125 int ee = 0;
1126 for (; ee < 4; ee++)
1127 P[ee] = NBEDGE_H1(p[ee]);
1128
1129 int n0 = 0;
1130 int n1 = 1;
1131 int n2 = 2;
1132 int n3 = 3;
1133
1134 int ii = 0;
1135 for (; ii < GDIM; ii++) {
1136 int node_shift = ii * 4;
1137 int node_diff_shift = 2 * node_shift;
1138
1139 double shape0 = N[node_shift + n0];
1140 double shape1 = N[node_shift + n1];
1141 double shape2 = N[node_shift + n2];
1142 double shape3 = N[node_shift + n3];
1143
1144 double ksi01 = (shape1 + shape2 - shape0 - shape3) * sense[n0];
1145 double ksi12 = (shape2 + shape3 - shape1 - shape0) * sense[n1];
1146 double ksi23 = (shape3 + shape0 - shape2 - shape1) * sense[n2];
1147 double ksi30 = (shape0 + shape1 - shape3 - shape2) * sense[n3];
1148
1149 double extrude_zeta01 = shape0 + shape1;
1150 double extrude_ksi12 = shape1 + shape2;
1151 double extrude_zeta23 = shape2 + shape3;
1152 double extrude_ksi30 = shape0 + shape3;
1153
1154 double bubble_ksi = extrude_ksi12 * extrude_ksi30;
1155 double bubble_zeta = extrude_zeta01 * extrude_zeta23;
1156
1157 double diff_ksi01[2], diff_ksi12[2], diff_ksi23[2], diff_ksi30[2];
1158 double diff_extrude_zeta01[2];
1159 double diff_extrude_ksi12[2];
1160 double diff_extrude_zeta23[2];
1161 double diff_extrude_ksi30[2];
1162 double diff_bubble_ksi[2];
1163 double diff_bubble_zeta[2];
1164
1165 int d = 0;
1166 for (; d < 2; d++) {
1167 double diff_shape0 = diffN[node_diff_shift + 2 * n0 + d];
1168 double diff_shape1 = diffN[node_diff_shift + 2 * n1 + d];
1169 double diff_shape2 = diffN[node_diff_shift + 2 * n2 + d];
1170 double diff_shape3 = diffN[node_diff_shift + 2 * n3 + d];
1171 diff_ksi01[d] =
1172 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) * sense[n0];
1173 diff_ksi12[d] =
1174 (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) * sense[n1];
1175 diff_ksi23[d] =
1176 (diff_shape3 + diff_shape0 - diff_shape2 - diff_shape1) * sense[n2];
1177 diff_ksi30[d] =
1178 (diff_shape0 + diff_shape1 - diff_shape3 - diff_shape2) * sense[n3];
1179 diff_extrude_zeta01[d] = diff_shape0 + diff_shape1;
1180 diff_extrude_ksi12[d] = diff_shape1 + diff_shape2;
1181 diff_extrude_zeta23[d] = diff_shape2 + diff_shape3;
1182 diff_extrude_ksi30[d] = diff_shape0 + diff_shape3;
1183 diff_bubble_ksi[d] = diff_extrude_ksi12[d] * extrude_ksi30 +
1184 extrude_ksi12 * diff_extrude_ksi30[d];
1185 diff_bubble_zeta[d] = diff_extrude_zeta01[d] * extrude_zeta23 +
1186 extrude_zeta01 * diff_extrude_zeta23[d];
1187 }
1188
1189 double L01[p[0] + 1], L12[p[1] + 1], L23[p[2] + 1], L30[p[3] + 1];
1190 double diffL01[2 * (p[0] + 1)], diffL12[2 * (p[1] + 1)],
1191 diffL23[2 * (p[2] + 1)], diffL30[2 * (p[3] + 1)];
1192 ierr = base_polynomials(p[0], ksi01, diff_ksi01, L01, diffL01, 2);
1193 CHKERRQ(ierr);
1194 ierr = base_polynomials(p[1], ksi12, diff_ksi12, L12, diffL12, 2);
1195 CHKERRQ(ierr);
1196 ierr = base_polynomials(p[2], ksi23, diff_ksi23, L23, diffL23, 2);
1197 CHKERRQ(ierr);
1198 ierr = base_polynomials(p[3], ksi30, diff_ksi30, L30, diffL30, 2);
1199 CHKERRQ(ierr);
1200
1201 int shift;
1202 if (edgeN != NULL) {
1203 // edge01
1204 shift = ii * (P[0]);
1205 cblas_dcopy(P[0], L01, 1, &edgeN01[shift], 1);
1206 cblas_dscal(P[0], bubble_ksi * extrude_zeta01, &edgeN01[shift], 1);
1207 // edge12
1208 shift = ii * (P[1]);
1209 cblas_dcopy(P[1], L12, 1, &edgeN12[shift], 1);
1210 cblas_dscal(P[1], bubble_zeta * extrude_ksi12, &edgeN12[shift], 1);
1211 // edge23
1212 shift = ii * (P[2]);
1213 cblas_dcopy(P[2], L23, 1, &edgeN23[shift], 1);
1214 cblas_dscal(P[2], bubble_ksi * extrude_zeta23, &edgeN23[shift], 1);
1215 // edge30
1216 shift = ii * (P[3]);
1217 cblas_dcopy(P[3], L30, 1, &edgeN30[shift], 1);
1218 cblas_dscal(P[3], bubble_zeta * extrude_ksi30, &edgeN30[shift], 1);
1219 }
1220 if (diff_edgeN != NULL) {
1221 if (P[0] > 0) {
1222 // edge01
1223 shift = ii * (P[0]);
1224 bzero(&diff_edgeN01[2 * shift], sizeof(double) * 2 * (P[0]));
1225 int d = 0;
1226 for (; d != 2; ++d) {
1227 cblas_daxpy(P[0], bubble_ksi * extrude_zeta01,
1228 &diffL01[d * (p[0] + 1)], 1, &diff_edgeN01[2 * shift + d],
1229 2);
1230 cblas_daxpy(P[0],
1231 diff_bubble_ksi[d] * extrude_zeta01 +
1232 bubble_ksi * diff_extrude_zeta01[d],
1233 L01, 1, &diff_edgeN01[2 * shift + d], 2);
1234 }
1235 }
1236 if (P[1] > 0) {
1237 // edge12
1238 shift = ii * (P[1]);
1239 bzero(&diff_edgeN12[2 * shift], sizeof(double) * 2 * (P[1]));
1240 int d = 0;
1241 for (; d != 2; ++d) {
1242 cblas_daxpy(P[1], bubble_zeta * extrude_ksi12,
1243 &diffL12[d * (p[1] + 1)], 1, &diff_edgeN12[2 * shift + d],
1244 2);
1245 cblas_daxpy(P[1],
1246 diff_bubble_zeta[d] * extrude_ksi12 +
1247 bubble_zeta * diff_extrude_ksi12[d],
1248 L12, 1, &diff_edgeN12[2 * shift + d], 2);
1249 }
1250 }
1251 if (P[2] > 0) {
1252 // edge23
1253 shift = ii * (P[2]);
1254 bzero(&diff_edgeN23[2 * shift], sizeof(double) * 2 * (P[2]));
1255 int d = 0;
1256 for (; d != 2; ++d) {
1257 cblas_daxpy(P[2], bubble_ksi * extrude_zeta23,
1258 &diffL23[d * (p[2] + 1)], 1, &diff_edgeN23[2 * shift + d],
1259 2);
1260 cblas_daxpy(P[2],
1261 diff_bubble_ksi[d] * extrude_zeta23 +
1262 bubble_ksi * diff_extrude_zeta23[d],
1263 L23, 1, &diff_edgeN23[2 * shift + d], 2);
1264 }
1265 }
1266 if (P[3] > 0) {
1267 // edge30
1268 shift = ii * (P[3]);
1269 bzero(&diff_edgeN30[2 * shift], sizeof(double) * 2 * (P[3]));
1270 int d = 0;
1271 for (; d != 2; ++d) {
1272 cblas_daxpy(P[3], bubble_zeta * extrude_ksi30,
1273 &diffL30[d * (p[3] + 1)], 1, &diff_edgeN30[2 * shift + d],
1274 2);
1275 cblas_daxpy(P[3],
1276 diff_bubble_zeta[d] * extrude_ksi30 +
1277 bubble_zeta * diff_extrude_ksi30[d],
1278 L30, 1, &diff_edgeN30[2 * shift + d], 2);
1279 }
1280 }
1281 }
1282 }
1284}
static PetscErrorCode ierr
Definition: h1.c:25
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
const int N
Definition: speed_test.cpp:3

◆ 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 284 of file h1.c.

289 {
291
292 int P[6];
293 int ee = 0;
294 for (; ee < 6; ee++)
295 P[ee] = NBEDGE_H1(p[ee]);
296 int edges_nodes[2 * 6] = {0, 1, 1, 2, 2, 0, 0, 3, 1, 3, 2, 3};
297 double diff_ksi01[3], diff_ksi12[3], diff_ksi20[3];
298 double diff_ksi03[3], diff_ksi13[3], diff_ksi23[3];
299 double *edges_diff_ksi[6] = {diff_ksi01, diff_ksi12, diff_ksi20,
300 diff_ksi03, diff_ksi13, diff_ksi23};
301 for (ee = 0; ee < 6; ee++) {
302 int dd = 0;
303 for (; dd < 3; dd++) {
304 edges_diff_ksi[ee][dd] = (diffN[edges_nodes[2 * ee + 1] * 3 + dd] -
305 diffN[edges_nodes[2 * ee + 0] * 3 + dd]) *
306 sense[ee];
307 }
308 }
309 int ii = 0;
310 for (; ii < GDIM; ii++) {
311 int node_shift = ii * 4;
312 double edges_ksi[6];
313 for (ee = 0; ee < 6; ee++) {
314 edges_ksi[ee] = (N[node_shift + edges_nodes[2 * ee + 1]] -
315 N[node_shift + edges_nodes[2 * ee + 0]]) *
316 sense[ee];
317 }
318 for (ee = 0; ee < 6; ee++) {
319 if (P[ee] == 0)
320 continue;
321 double L[p[ee] + 1], diffL[3 * (p[ee] + 1)];
322 ierr = base_polynomials(p[ee], edges_ksi[ee], edges_diff_ksi[ee], L,
323 diffL, 3);
324 CHKERRQ(ierr);
325 double v = N[node_shift + edges_nodes[2 * ee + 0]] *
326 N[node_shift + edges_nodes[2 * ee + 1]];
327 if (edgeN != NULL)
328 if (edgeN[ee] != NULL) {
329 int shift = ii * P[ee];
330 {
331 double *edge_n_ptr = &edgeN[ee][shift];
332 double *l_ptr = L;
333 double scalar = N[node_shift + edges_nodes[2 * ee + 0]] *
334 N[node_shift + edges_nodes[2 * ee + 1]];
335 int size = P[ee];
336 for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
337 *edge_n_ptr = (*l_ptr) * scalar;
338 ++edge_n_ptr;
339 }
340 }
341 }
342 if (diff_edgeN != NULL)
343 if (diff_edgeN[ee] != NULL) {
344 int shift = ii * P[ee];
345 {
346 double *diff_edge_n_ptr = &diff_edgeN[ee][3 * shift];
347 double *diff_l_x = &diffL[0 * (p[ee] + 1)];
348 double *diff_l_y = &diffL[1 * (p[ee] + 1)];
349 double *diff_l_z = &diffL[2 * (p[ee] + 1)];
350 double *l_ptr = L;
351 double scalar_x = diffN[3 * edges_nodes[2 * ee + 0] + 0] *
352 N[node_shift + edges_nodes[2 * ee + 1]] +
353 N[node_shift + edges_nodes[2 * ee + 0]] *
354 diffN[3 * edges_nodes[2 * ee + 1] + 0];
355 double scalar_y = diffN[3 * edges_nodes[2 * ee + 0] + 1] *
356 N[node_shift + edges_nodes[2 * ee + 1]] +
357 N[node_shift + edges_nodes[2 * ee + 0]] *
358 diffN[3 * edges_nodes[2 * ee + 1] + 1];
359 double scalar_z = diffN[3 * edges_nodes[2 * ee + 0] + 2] *
360 N[node_shift + edges_nodes[2 * ee + 1]] +
361 N[node_shift + edges_nodes[2 * ee + 0]] *
362 diffN[3 * edges_nodes[2 * ee + 1] + 2];
363
364 int size = P[ee];
365 for (size_t jj = 0; jj != size;
366 ++jj, ++diff_l_x, ++diff_l_y, ++diff_l_z, ++l_ptr) {
367
368 *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
369 ++diff_edge_n_ptr;
370 *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
371 ++diff_edge_n_ptr;
372 *diff_edge_n_ptr = v * (*diff_l_z) + scalar_z * (*l_ptr);
373 ++diff_edge_n_ptr;
374
375 }
376 }
377
378 }
379 }
380 }
382}
static Index< 'L', 3 > L
double v
phase velocity of light in medium (cm/ns)
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

◆ 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 // edge 01
79 shift = ii * (P[0]);
80 {
81 double *edge_n_ptr = &edgeN01[shift];
82 double *l_ptr = L01;
83 double scalar = N[node_shift + 0] * N[node_shift + 1];
84 int size = P[0];
85 for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
86 *edge_n_ptr = (*l_ptr) * scalar;
87 ++edge_n_ptr;
88 }
89 }
90
91 // edge 12
92 shift = ii * (P[1]);
93 {
94 double *edge_n_ptr = &edgeN12[shift];
95 double *l_ptr = L12;
96 double scalar = N[node_shift + 1] * N[node_shift + 2];
97 int size = P[1];
98 for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
99 *edge_n_ptr = (*l_ptr) * scalar;
100 ++edge_n_ptr;
101 }
102 }
103
104 // edge 20
105 shift = ii * (P[2]);
106 {
107 double *edge_n_ptr = &edgeN20[shift];
108 double *l_ptr = L20;
109 double scalar = N[node_shift + 2] * N[node_shift + 0];
110 int size = P[2];
111 for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
112 *edge_n_ptr = (*l_ptr) * scalar;
113 ++edge_n_ptr;
114 }
115 }
116 }
117 if (diff_edgeN != NULL) {
118 if (P[0] > 0) {
119 // edge01
120 shift = ii * (P[0]);
121 {
122 double *diff_edge_n_ptr = &diff_edgeN01[2 * shift];
123 double *diff_l_x = &diffL01[0 * (p[0] + 1)];
124 double *diff_l_y = &diffL01[1 * (p[0] + 1)];
125 double scalar_x = diffN[2 * 0 + 0] * N[node_shift + 1] +
126 N[node_shift + 0] * diffN[2 * 1 + 0];
127 double scalar_y = diffN[2 * 0 + 1] * N[node_shift + 1] +
128 N[node_shift + 0] * diffN[2 * 1 + 1];
129 double v = N[node_shift + 0] * N[node_shift + 1];
130 double *l_ptr = L01;
131
132 int size = P[0];
133 for (size_t jj = 0; jj != size;
134 ++jj, ++diff_l_x, ++diff_l_y, ++l_ptr) {
135
136 *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
137 ++diff_edge_n_ptr;
138 *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
139 ++diff_edge_n_ptr;
140 }
141 }
142
143 }
144 if (P[1] > 0) {
145 // edge12
146 shift = ii * (P[1]);
147 {
148 double *diff_edge_n_ptr = &diff_edgeN12[2 * shift];
149 double *diff_l_x = &diffL12[0 * (p[1] + 1)];
150 double *diff_l_y = &diffL12[1 * (p[1] + 1)];
151 double scalar_x = diffN[2 * 1 + 0] * N[node_shift + 2] +
152 N[node_shift + 1] * diffN[2 * 2 + 0];
153 double scalar_y = diffN[2 * 1 + 1] * N[node_shift + 2] +
154 N[node_shift + 1] * diffN[2 * 2 + 1];
155 double v = N[node_shift + 1] * N[node_shift + 2];
156 double *l_ptr = L12;
157
158 int size = P[1];
159 for (size_t jj = 0; jj != size;
160 ++jj, ++diff_l_x, ++diff_l_y, ++l_ptr) {
161
162 *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
163 ++diff_edge_n_ptr;
164 *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
165 ++diff_edge_n_ptr;
166 }
167 }
168
169 }
170 if (P[2] > 0) {
171 // edge20
172 shift = ii * (P[2]);
173
174 {
175 double *diff_edge_n_ptr = &diff_edgeN20[2 * shift];
176 double *diff_l_x = &diffL20[0 * (p[2] + 1)];
177 double *diff_l_y = &diffL20[1 * (p[2] + 1)];
178 double scalar_x = diffN[2 * 2 + 0] * N[node_shift + 0] +
179 N[node_shift + 2] * diffN[2 * 0 + 0];
180 double scalar_y = diffN[2 * 2 + 1] * N[node_shift + 0] +
181 N[node_shift + 2] * diffN[2 * 0 + 1];
182 double v = N[node_shift + 2] * N[node_shift + 0];
183 double *l_ptr = L20;
184
185 int size = P[2];
186 for (size_t jj = 0; jj != size;
187 ++jj, ++diff_l_x, ++diff_l_y, ++l_ptr) {
188
189 *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
190 ++diff_edge_n_ptr;
191 *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
192 ++diff_edge_n_ptr;
193 }
194 }
195
196 }
197 }
198 }
200}

◆ H1_FaceGradientOfDeformation_hierarchical()

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

Definition at line 630 of file h1.c.

632 {
634 int col, row = 0;
635 for (; row < 3; row++)
636 for (col = 0; col < 3; col++)
637 F[3 * row + col] =
638 cblas_ddot(NBFACETRI_H1(p), &diffN[col], 3, &dofs[row], 3);
640}
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.

◆ H1_FaceShapeDiffMBTETinvJ()

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

Definition at line 584 of file h1.c.

586 {
588 int ff = 0, ii, gg;
589 for (; ff < 4; ff++) {
590 for (ii = 0; ii < NBFACETRI_H1(p[ff]); ii++) {
591 for (gg = 0; gg < GDIM; gg++) {
592 int shift1 = NBFACETRI_H1(base_p[ff]) * gg;
593 int shift2 = NBFACETRI_H1(p[ff]) * gg;
594 cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
595 &(face_diffN[ff])[3 * shift1 + 3 * ii], 1, 0.,
596 &(face_diffNinvJac[ff])[3 * shift2 + 3 * ii], 1);
597 }
598 }
599 }
601}

◆ 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 383 of file h1.c.

388 {
390
391 int P[4];
392 int ff = 0;
393 for (; ff < 4; ff++)
394 P[ff] = NBFACETRI_H1(p[ff]);
395 double diff_ksiL0F0[3], diff_ksiL1F0[3];
396 double diff_ksiL0F1[3], diff_ksiL1F1[3];
397 double diff_ksiL0F2[3], diff_ksiL1F2[3];
398 double diff_ksiL0F3[3], diff_ksiL1F3[3];
399 double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL1F0, diff_ksiL0F1,
400 diff_ksiL1F1, diff_ksiL0F2, diff_ksiL1F2,
401 diff_ksiL0F3, diff_ksiL1F3};
402 int dd = 0;
403 for (; dd < 3; dd++) {
404 for (ff = 0; ff < 4; ff++) {
405 diff_ksi_faces[ff * 2 + 0][dd] =
406 (diffN[faces_nodes[3 * ff + 1] * 3 + dd] -
407 diffN[faces_nodes[3 * ff + 0] * 3 + dd]);
408 diff_ksi_faces[ff * 2 + 1][dd] =
409 (diffN[faces_nodes[3 * ff + 2] * 3 + dd] -
410 diffN[faces_nodes[3 * ff + 0] * 3 + dd]);
411 }
412 }
413 int ii = 0;
414 for (; ii < GDIM; ii++) {
415 int node_shift = ii * 4;
416 double ksi_faces[8];
417 for (ff = 0; ff < 4; ff++) {
418 ksi_faces[2 * ff + 0] = N[node_shift + faces_nodes[3 * ff + 1]] -
419 N[node_shift + faces_nodes[3 * ff + 0]];
420 ksi_faces[2 * ff + 1] = N[node_shift + faces_nodes[3 * ff + 2]] -
421 N[node_shift + faces_nodes[3 * ff + 0]];
422 }
423 int shift;
424 for (ff = 0; ff < 4; ff++) {
425 if (P[ff] == 0)
426 continue;
427 double L0[p[ff] + 1], L1[p[ff] + 1];
428 double diffL0[3 * (p[ff] + 1)], diffL1[3 * (p[ff] + 1)];
429 ierr = base_polynomials(p[ff], ksi_faces[ff * 2 + 0],
430 diff_ksi_faces[ff * 2 + 0], L0, diffL0, 3);
431 CHKERRQ(ierr);
432 ierr = base_polynomials(p[ff], ksi_faces[ff * 2 + 1],
433 diff_ksi_faces[ff * 2 + 1], L1, diffL1, 3);
434 CHKERRQ(ierr);
435 double v = N[node_shift + faces_nodes[3 * ff + 0]] *
436 N[node_shift + faces_nodes[3 * ff + 1]] *
437 N[node_shift + faces_nodes[3 * ff + 2]];
438 double v2[3] = {0, 0, 0};
439 dd = 0;
440 double n1n2 = N[node_shift + faces_nodes[3 * ff + 1]] *
441 N[node_shift + faces_nodes[3 * ff + 2]];
442 double n0n2 = N[node_shift + faces_nodes[3 * ff + 0]] *
443 N[node_shift + faces_nodes[3 * ff + 2]];
444 double n0n1 = N[node_shift + faces_nodes[3 * ff + 0]] *
445 N[node_shift + faces_nodes[3 * ff + 1]];
446 for (; dd < 3; dd++) {
447 v2[dd] = diffN[3 * faces_nodes[3 * ff + 0] + dd] * n1n2 +
448 diffN[3 * faces_nodes[3 * ff + 1] + dd] * n0n2 +
449 diffN[3 * faces_nodes[3 * ff + 2] + dd] * n0n1;
450 }
451 shift = ii * P[ff];
452 int jj = 0;
453 int oo = 0;
454 for (; oo <= (p[ff] - 3); oo++) {
455 int pp0 = 0;
456 for (; pp0 <= oo; pp0++) {
457 int pp1 = oo - pp0;
458 if (pp1 >= 0) {
459 if (faceN != NULL)
460 if (faceN[ff] != NULL) {
461 faceN[ff][shift + jj] = v * L0[pp0] * L1[pp1];
462 }
463 if (diff_faceN != NULL)
464 if (diff_faceN[ff] != NULL) {
465 dd = 0;
466 double L0L1 = L0[pp0] * L1[pp1];
467 for (; dd < 3; dd++) {
468 diff_faceN[ff][3 * shift + 3 * jj + dd] =
469 (L0[pp0] * diffL1[dd * (p[ff] + 1) + pp1] +
470 diffL0[dd * (p[ff] + 1) + pp0] * L1[pp1]) *
471 v;
472 diff_faceN[ff][3 * shift + 3 * jj + dd] += L0L1 * v2[dd];
473 }
474 }
475 jj++;
476 }
477 }
478 }
479 if (jj != P[ff])
480 SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P[ff]);
481 }
482 }
484}

◆ 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 201 of file h1.c.

206 {
208
209 int P = NBFACETRI_H1(p);
210 if (P == 0)
212 double diff_ksiL0[2], diff_ksiL1[2];
213 double *diff_ksi_faces[] = {diff_ksiL0, diff_ksiL1};
214 int dd = 0;
215 if (diff_faceN != NULL) {
216 for (; dd < 2; dd++) {
217 diff_ksi_faces[0][dd] =
218 diffN[face_nodes[1] * 2 + dd] - diffN[face_nodes[0] * 2 + dd];
219 diff_ksi_faces[1][dd] =
220 diffN[face_nodes[2] * 2 + dd] - diffN[face_nodes[0] * 2 + dd];
221 }
222 }
223 int ii = 0;
224 for (; ii < GDIM; ii++) {
225 int node_shift = ii * 3;
226 double ksi_faces[2];
227 ksi_faces[0] =
228 N[node_shift + face_nodes[1]] - N[node_shift + face_nodes[0]];
229 ksi_faces[1] =
230 N[node_shift + face_nodes[2]] - N[node_shift + face_nodes[0]];
231 double L0[p + 1], L1[p + 1];
232 double diffL0[2 * (p + 1)], diffL1[2 * (p + 1)];
233 ierr = base_polynomials(p, ksi_faces[0], diff_ksi_faces[0], L0, diffL0, 2);
234 CHKERRQ(ierr);
235 ierr = base_polynomials(p, ksi_faces[1], diff_ksi_faces[1], L1, diffL1, 2);
236 CHKERRQ(ierr);
237 double v = N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[1]] *
238 N[node_shift + face_nodes[2]];
239 double v2[2] = {0, 0};
240 if (diff_faceN != NULL) {
241 double n1n2 =
242 N[node_shift + face_nodes[1]] * N[node_shift + face_nodes[2]];
243 double n0n2 =
244 N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[2]];
245 double n0n1 =
246 N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[1]];
247 dd = 0;
248 for (; dd < 2; dd++) {
249 v2[dd] = diffN[face_nodes[0] * 2 + dd] * n1n2 +
250 diffN[face_nodes[1] * 2 + dd] * n0n2 +
251 diffN[face_nodes[2] * 2 + dd] * n0n1;
252 }
253 }
254 int shift = ii * P;
255 int jj = 0;
256 int oo = 0;
257 for (; oo <= (p - 3); oo++) {
258 int pp0 = 0;
259 for (; pp0 <= oo; pp0++) {
260 int pp1 = oo - pp0;
261 if (pp1 >= 0) {
262 if (faceN != NULL) {
263 faceN[shift + jj] = v * L0[pp0] * L1[pp1];
264 }
265 if (diff_faceN != NULL) {
266 dd = 0;
267 for (; dd < 2; dd++) {
268 double *val = &diff_faceN[2 * shift + 2 * jj + dd];
269 *val = (L0[pp0] * diffL1[dd * (p + 1) + pp1] +
270 diffL0[dd * (p + 1) + pp0] * L1[pp1]) *
271 v;
272 *val += L0[pp0] * L1[pp1] * v2[dd];
273 }
274 }
275 jj++;
276 }
277 }
278 }
279 if (jj != P)
280 SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P);
281 }
283}

◆ 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 652 of file h1.c.

657 {
659 // TODO: return separately components of the tensor product between two edges
660
661 int P[3];
662 int ff = 0;
663 for (; ff < 3; ff++)
664 P[ff] = NBFACEQUAD_H1(p[ff]);
665
666 double ksi_faces[6];
667 double diff_ksiL0F0[3], diff_ksiL3F0[3];
668 double diff_ksiL0F1[3], diff_ksiL3F1[3];
669 double diff_ksiL0F2[3], diff_ksiL3F2[3];
670 double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL3F0, diff_ksiL0F1,
671 diff_ksiL3F1, diff_ksiL0F2, diff_ksiL3F2};
672 int ii = 0;
673 for (; ii < GDIM; ii++) {
674 int node_shift = ii * 6;
675 int node_diff_shift = 3 * node_shift;
676 int ff = 0;
677 for (; ff < 3; ff++) {
678 if (P[ff] == 0)
679 continue;
680 int n0 = faces_nodes[4 * ff + 0];
681 int n1 = faces_nodes[4 * ff + 1];
682 int n2 = faces_nodes[4 * ff + 2];
683 int n3 = faces_nodes[4 * ff + 3];
684 int e0 = 2 * ff + 0;
685 int e1 = 2 * ff + 1;
686 double ksi0 = N[node_shift + n0] + N[node_shift + n3];
687 double ksi1 = N[node_shift + n1] + N[node_shift + n2];
688 double eta0 = N[node_shift + n0] + N[node_shift + n1];
689 double eta1 = N[node_shift + n2] + N[node_shift + n3];
690 ksi_faces[e0] = ksi1 - ksi0;
691 ksi_faces[e1] = eta1 - eta0;
692
693 int dd = 0;
694 for (; dd < 3; dd++) {
695 double diff_ksi0 = diffN[node_diff_shift + 3 * n0 + dd] +
696 diffN[node_diff_shift + 3 * n3 + dd];
697 double diff_ksi1 = diffN[node_diff_shift + 3 * n1 + dd] +
698 diffN[node_diff_shift + 3 * n2 + dd];
699 double diff_eta0 = diffN[node_diff_shift + 3 * n0 + dd] +
700 diffN[node_diff_shift + 3 * n1 + dd];
701 double diff_eta1 = diffN[node_diff_shift + 3 * n2 + dd] +
702 diffN[node_diff_shift + 3 * n3 + dd];
703 diff_ksi_faces[e0][dd] = diff_ksi1 - diff_ksi0;
704 diff_ksi_faces[e1][dd] = diff_eta1 - diff_eta0;
705 }
706 double L0[p[ff] + 1], L1[p[ff] + 1];
707 double diffL0[3 * (p[ff] + 1)], diffL1[3 * (p[ff] + 1)];
708 ierr = base_polynomials(p[ff], ksi_faces[e0], diff_ksi_faces[e0], L0,
709 diffL0, 3);
710 CHKERRQ(ierr);
711 ierr = base_polynomials(p[ff], ksi_faces[e1], diff_ksi_faces[e1], L1,
712 diffL1, 3);
713 CHKERRQ(ierr);
714
715 double v = N[node_shift + n0] * N[node_shift + n2] +
716 N[node_shift + n1] * N[node_shift + n3];
717
718 double diff_v[3];
719 dd = 0;
720 for (; dd < 3; ++dd)
721 diff_v[dd] =
722
723 diffN[node_diff_shift + 3 * n0 + dd] * N[node_shift + n2] +
724 N[node_shift + n0] * diffN[node_diff_shift + 3 * n2 + dd] +
725
726 diffN[node_diff_shift + 3 * n1 + dd] * N[node_shift + n3] +
727 N[node_shift + n1] * diffN[node_diff_shift + 3 * n3 + dd];
728
729 int shift;
730 shift = ii * P[ff];
731
732 if (faceN != NULL) {
733 if (faceN[ff] != NULL) {
734 int jj = 0;
735 int oo = 0;
736 for (; oo <= p[ff] - 2; ++oo) {
737 int pp0 = 0;
738 for (; pp0 < oo; ++pp0) {
739 int pp1 = oo;
740 faceN[ff][shift + jj] = L0[pp0] * L1[pp1] * v;
741 ++jj;
742 faceN[ff][shift + jj] = L0[pp1] * L1[pp0] * v;
743 ++jj;
744 }
745 faceN[ff][shift + jj] = L0[oo] * L1[oo] * v;
746 ++jj;
747 }
748 if (jj != P[ff])
749 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
750 "Inconsistent implementation (bug in the code) %d != %d",
751 jj, P);
752 }
753 }
754
755 if (diff_faceN != NULL) {
756 if (diff_faceN[ff] != NULL) {
757 int jj = 0;
758 int oo = 0;
759 for (; oo <= p[ff] - 2; ++oo) {
760 int pp0 = 0;
761 for (; pp0 < oo; ++pp0) {
762 int pp1 = oo;
763 int dd;
764 for (dd = 0; dd < 3; dd++) {
765 diff_faceN[ff][3 * shift + 3 * jj + dd] =
766 (L0[pp0] * diffL1[dd * (p[ff] + 1) + pp1] +
767 diffL0[dd * (p[ff] + 1) + pp0] * L1[pp1]) *
768 v +
769 L0[pp0] * L1[pp1] * diff_v[dd];
770 }
771 ++jj;
772 for (dd = 0; dd < 3; dd++) {
773 diff_faceN[ff][3 * shift + 3 * jj + dd] =
774 (L0[pp1] * diffL1[dd * (p[ff] + 1) + pp0] +
775 diffL0[dd * (p[ff] + 1) + pp1] * L1[pp0]) *
776 v +
777 L0[pp1] * L1[pp0] * diff_v[dd];
778 }
779 ++jj;
780 }
781 for (dd = 0; dd < 3; dd++) {
782 diff_faceN[ff][3 * shift + 3 * jj + dd] =
783 (L0[oo] * diffL1[dd * (p[ff] + 1) + oo] +
784 diffL0[dd * (p[ff] + 1) + oo] * L1[oo]) *
785 v +
786 L0[oo] * L1[oo] * diff_v[dd];
787 }
788 ++jj;
789 }
790 if (jj != P[ff])
791 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
792 "Inconsistent implementation (bug in the code) %d != %d",
793 jj, P);
794 }
795 }
796 }
797 }
799}
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.

◆ 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 969 of file h1.c.

974 {
976 // TODO: return separately components of the tensor product between two edges
977
978 int P = NBFACEQUAD_H1(p);
979 if (P == 0)
981 double diff_ksiL0F0[2], diff_ksiL1F0[2];
982 double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL1F0};
983 int ii = 0;
984 for (; ii < GDIM; ii++) {
985 int node_shift = ii * 4;
986 int node_diff_shift = 2 * node_shift;
987
988 int n0 = faces_nodes[0];
989 int n1 = faces_nodes[1];
990 int n2 = faces_nodes[2];
991 int n3 = faces_nodes[3];
992 double ksi0 = N[node_shift + n0] + N[node_shift + n3];
993 double ksi1 = N[node_shift + n1] + N[node_shift + n2];
994 double eta0 = N[node_shift + n0] + N[node_shift + n1];
995 double eta1 = N[node_shift + n2] + N[node_shift + n3];
996 double ksi_faces = ksi1 - ksi0;
997 double eta_faces = eta1 - eta0;
998 double diff_ksi_faces[2];
999 double diff_eta_faces[2];
1000 int dd = 0;
1001 for (; dd < 2; dd++) {
1002 double diff_ksi0 = diffN[node_diff_shift + 2 * n0 + dd] +
1003 diffN[node_diff_shift + 2 * n3 + dd];
1004 double diff_ksi1 = diffN[node_diff_shift + 2 * n1 + dd] +
1005 diffN[node_diff_shift + 2 * n2 + dd];
1006 double diff_eta0 = diffN[node_diff_shift + 2 * n0 + dd] +
1007 diffN[node_diff_shift + 2 * n1 + dd];
1008 double diff_eta1 = diffN[node_diff_shift + 2 * n2 + dd] +
1009 diffN[node_diff_shift + 2 * n3 + dd];
1010
1011 diff_ksi_faces[dd] = diff_ksi1 - diff_ksi0;
1012 diff_eta_faces[dd] = diff_eta1 - diff_eta0;
1013 }
1014 double L0[p + 1], L1[p + 1];
1015 double diffL0[2 * (p + 1)], diffL1[2 * (p + 1)];
1016 ierr = base_polynomials(p, ksi_faces, diff_ksi_faces, L0, diffL0, 2);
1017 CHKERRQ(ierr);
1018 ierr = base_polynomials(p, eta_faces, diff_eta_faces, L1, diffL1, 2);
1019 CHKERRQ(ierr);
1020
1021 double v = N[node_shift + n0] * N[node_shift + n2] +
1022 N[node_shift + n1] * N[node_shift + n3];
1023
1024 double diff_v[2];
1025 dd = 0;
1026 for (; dd < 2; ++dd)
1027 diff_v[dd] =
1028
1029 diffN[node_diff_shift + 2 * n0 + dd] * N[node_shift + n2] +
1030 N[node_shift + n0] * diffN[node_diff_shift + 2 * n2 + dd] +
1031
1032 diffN[node_diff_shift + 2 * n1 + dd] * N[node_shift + n3] +
1033 N[node_shift + n1] * diffN[node_diff_shift + 2 * n3 + dd];
1034
1035 const int shift = ii * P;
1036
1037 if (faceN != NULL) {
1038 int jj = 0;
1039 int oo = 0;
1040 for (; oo <= p - 2; ++oo) {
1041 int pp0 = 0;
1042 for (; pp0 < oo; ++pp0) {
1043 int pp1 = oo;
1044 faceN[shift + jj] = L0[pp0] * L1[pp1] * v;
1045 ++jj;
1046 faceN[shift + jj] = L0[pp1] * L1[pp0] * v;
1047 ++jj;
1048 }
1049 faceN[shift + jj] = L0[oo] * L1[oo] * v;
1050 ++jj;
1051 }
1052 if (jj != P)
1053 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1054 "Inconsistent implementation (bug in the code) %d != %d", jj,
1055 P);
1056 }
1057
1058 if (diff_faceN != NULL) {
1059 int jj = 0;
1060 int oo = 0;
1061 for (; oo <= p - 2; ++oo) {
1062 int pp0 = 0;
1063 for (; pp0 < oo; ++pp0) {
1064 int pp1 = oo;
1065 int dd;
1066 for (dd = 0; dd < 2; dd++) {
1067 diff_faceN[2 * shift + 2 * jj + dd] =
1068 (L0[pp0] * diffL1[dd * (p + 1) + pp1] +
1069 diffL0[dd * (p + 1) + pp0] * L1[pp1]) *
1070 v +
1071 L0[pp0] * L1[pp1] * diff_v[dd];
1072 }
1073 ++jj;
1074 for (dd = 0; dd < 2; dd++) {
1075 diff_faceN[2 * shift + 2 * jj + dd] =
1076 (L0[pp1] * diffL1[dd * (p + 1) + pp0] +
1077 diffL0[dd * (p + 1) + pp1] * L1[pp0]) *
1078 v +
1079 L0[pp1] * L1[pp0] * diff_v[dd];
1080 }
1081 ++jj;
1082 }
1083 for (dd = 0; dd < 2; dd++) {
1084 diff_faceN[2 * shift + 2 * jj + dd] =
1085 (L0[oo] * diffL1[dd * (p + 1) + oo] +
1086 diffL0[dd * (p + 1) + oo] * L1[oo]) *
1087 v +
1088 L0[oo] * L1[oo] * diff_v[dd];
1089 }
1090 ++jj;
1091 }
1092 if (jj != P)
1093 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1094 "Inconsistent implementation (bug in the code) %d != %d", jj,
1095 P);
1096 }
1097 }
1099}

◆ H1_VolumeGradientOfDeformation_hierarchical()

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

Definition at line 641 of file h1.c.

643 {
645 int col, row = 0;
646 for (; row < 3; row++)
647 for (col = 0; col < 3; col++)
648 F[3 * row + col] =
649 cblas_ddot(NBVOLUMETET_H1(p), &diffN[col], 3, &dofs[row], 3);
651}
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.

◆ H1_VolumeShapeDiffMBTETinvJ()

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

Definition at line 602 of file h1.c.

605 {
607 int ii, gg;
608 for (ii = 0; ii < NBVOLUMETET_H1(p); ii++) {
609 for (gg = 0; gg < GDIM; gg++) {
610 int shift1 = NBVOLUMETET_H1(base_p) * gg;
611 int shift2 = NBVOLUMETET_H1(p) * gg;
612 cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
613 &(volume_diffN)[3 * shift1 + 3 * ii], 1, 0.,
614 &(volume_diffNinvJac)[3 * shift2 + 3 * ii], 1);
615 }
616 }
618}

◆ 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 800 of file h1.c.

805 {
807
808 int P = NBVOLUMEPRISM_H1(p);
809 if (P == 0)
811 double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
812 int ii = 0;
813 for (; ii < GDIM; ii++) {
814 int node_shift = ii * 6;
815 int node_diff_shift = ii * 18;
816
817 double n0 = N[node_shift + 0];
818 double n1 = N[node_shift + 1];
819 double n2 = N[node_shift + 2];
820 double n3 = N[node_shift + 3];
821 double n4 = N[node_shift + 4];
822 double n5 = N[node_shift + 5];
823
824 double ksiL0 = n1 + n4 - n0 - n3;
825 double ksiL1 = n2 + n5 - n0 - n3;
826 double ksiL2 = (n3 + n4 + n5) - (n0 + n1 + n2);
827
828 int dd = 0;
829 for (; dd < 3; dd++) {
830 double diff_n0 = diffN[node_diff_shift + 3 * 0 + dd];
831 double diff_n1 = diffN[node_diff_shift + 3 * 1 + dd];
832 double diff_n2 = diffN[node_diff_shift + 3 * 2 + dd];
833 double diff_n3 = diffN[node_diff_shift + 3 * 3 + dd];
834 double diff_n4 = diffN[node_diff_shift + 3 * 4 + dd];
835 double diff_n5 = diffN[node_diff_shift + 3 * 5 + dd];
836 diff_ksiL0[dd] = (diff_n1 + diff_n4) - (diff_n0 + diff_n3);
837 diff_ksiL1[dd] = (diff_n2 + diff_n5) - (diff_n0 + diff_n3);
838 diff_ksiL2[dd] =
839 (diff_n3 + diff_n4 + diff_n5) - (diff_n0 + diff_n1 + diff_n2);
840 }
841
842 double L0[p + 1], L1[p + 1], L2[p + 1];
843 double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
844 ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
845 CHKERRQ(ierr);
846 ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
847 CHKERRQ(ierr);
848 ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
849 CHKERRQ(ierr);
850
851 double v_tri = (n0 + n3) * (n1 + n4) * (n2 + n5);
852 double v_edge = (n0 + n1 + n2) * (n3 + n4 + n5);
853 double v = v_tri * v_edge;
854
855 double diff_v_tri[3];
856 double diff_v_edge[3];
857 double diff_v[3];
858 dd = 0;
859 for (; dd < 3; dd++) {
860 double diff_n0 = diffN[node_diff_shift + 3 * 0 + dd];
861 double diff_n1 = diffN[node_diff_shift + 3 * 1 + dd];
862 double diff_n2 = diffN[node_diff_shift + 3 * 2 + dd];
863 double diff_n3 = diffN[node_diff_shift + 3 * 3 + dd];
864 double diff_n4 = diffN[node_diff_shift + 3 * 4 + dd];
865 double diff_n5 = diffN[node_diff_shift + 3 * 5 + dd];
866 diff_v_tri[dd] = (diff_n0 + diff_n3) * (n1 + n4) * (n2 + n5) +
867 (diff_n1 + diff_n4) * (n0 + n3) * (n2 + n5) +
868 (diff_n2 + diff_n5) * (n0 + n3) * (n1 + n4);
869 diff_v_edge[dd] = (diff_n0 + diff_n1 + diff_n2) * (n3 + n4 + n5) +
870 (diff_n3 + diff_n4 + diff_n5) * (n0 + n1 + n2);
871 diff_v[dd] = diff_v_tri[dd] * v_edge + v_tri * diff_v_edge[dd];
872 }
873
874 int shift = ii * P;
875
876 if (volumeN != NULL) {
877 int jj = 0;
878 int oo, pp0, pp1, zz;
879 for (oo = 0; oo <= p - 3; ++oo) {
880 for (pp0 = 0; pp0 < oo; ++pp0) {
881 for (pp1 = 0; pp1 < oo; ++pp1) {
882 const int perm[3][3] = {
883 {oo, pp0, pp1}, {pp0, oo, pp1}, {pp0, pp1, oo}};
884 for (zz = 0; zz != 3; ++zz) {
885 volumeN[shift + jj] =
886 L0[perm[zz][0]] * L1[perm[zz][1]] * L2[perm[zz][2]] * v;
887 ++jj;
888 }
889 }
890 const int perm[3][3] = {{pp0, oo, oo}, {oo, pp0, oo}, {oo, oo, pp0}};
891 for (zz = 0; zz != 3; ++zz) {
892 volumeN[shift + jj] =
893 L0[perm[zz][0]] * L1[perm[zz][1]] * L2[perm[zz][2]] * v;
894 ++jj;
895 }
896 }
897 volumeN[shift + jj] = L0[oo] * L1[oo] * L2[oo] * v;
898 ++jj;
899 }
900 if (jj != P)
901 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
902 "Inconsistent implementation (bug in the code) %d != %d", jj,
903 P);
904 }
905
906 if (diff_volumeN != NULL) {
907 int jj = 0;
908 int oo, pp0, pp1, zz, dd;
909 for (oo = 0; oo <= p - 3; ++oo) {
910 for (pp0 = 0; pp0 < oo; ++pp0) {
911 for (pp1 = 0; pp1 < oo; ++pp1) {
912 const int perm[3][3] = {
913 {oo, pp0, pp1}, {pp0, oo, pp1}, {pp0, pp1, oo}};
914 for (zz = 0; zz != 3; ++zz) {
915 const int i = perm[zz][0];
916 const int j = perm[zz][1];
917 const int k = perm[zz][2];
918 for (dd = 0; dd != 3; ++dd) {
919 diff_volumeN[3 * shift + 3 * jj + dd] =
920 (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
921 L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
922 L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
923 v +
924 L0[i] * L1[j] * L2[k] * diff_v[dd];
925 }
926 ++jj;
927 }
928 }
929 const int perm[3][3] = {{pp0, oo, oo}, {oo, pp0, oo}, {oo, oo, pp0}};
930 for (zz = 0; zz != 3; ++zz) {
931 const int i = perm[zz][0];
932 const int j = perm[zz][1];
933 const int k = perm[zz][2];
934 for (dd = 0; dd != 3; ++dd) {
935 diff_volumeN[3 * shift + 3 * jj + dd] =
936 (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
937 L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
938 L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
939 v +
940 L0[i] * L1[j] * L2[k] * diff_v[dd];
941 }
942 ++jj;
943 }
944 }
945
946 const int i = oo;
947 const int j = oo;
948 const int k = oo;
949 for (dd = 0; dd != 3; ++dd) {
950 diff_volumeN[3 * shift + 3 * jj + dd] =
951 (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
952 L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
953 L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
954 v +
955 L0[i] * L1[j] * L2[k] * diff_v[dd];
956 }
957
958 ++jj;
959 }
960 if (jj != P)
961 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
962 "Inconsistent implementation (bug in the code) %d != %d", jj,
963 P);
964 }
965 }
967}
@ L2
field with C-1 continuity
Definition: definitions.h:101
#define NBVOLUMEPRISM_H1(P)
Number of base functions on prism for H1 space.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k

◆ 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 485 of file h1.c.

490 {
492
493 int P = NBVOLUMETET_H1(p);
494 if (P == 0)
496 double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
497 int dd = 0;
498 for (; dd < 3; dd++) {
499 diff_ksiL0[dd] = (diffN[1 * 3 + dd] - diffN[0 * 3 + dd]);
500 diff_ksiL1[dd] = (diffN[2 * 3 + dd] - diffN[0 * 3 + dd]);
501 diff_ksiL2[dd] = (diffN[3 * 3 + dd] - diffN[0 * 3 + dd]);
502 }
503 int ii = 0;
504 for (; ii < GDIM; ii++) {
505 int node_shift = ii * 4;
506 double ksiL0 = N[node_shift + 1] - N[node_shift + 0];
507 double ksiL1 = N[node_shift + 2] - N[node_shift + 0];
508 double ksiL2 = N[node_shift + 3] - N[node_shift + 0];
509 double L0[p + 1], L1[p + 1], L2[p + 1];
510 double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
511 ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
512 CHKERRQ(ierr);
513 ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
514 CHKERRQ(ierr);
515 ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
516 CHKERRQ(ierr);
517 double v = N[node_shift + 0] * N[node_shift + 1] * N[node_shift + 2] *
518 N[node_shift + 3];
519 double v2[3] = {0, 0, 0};
520 dd = 0;
521 for (; dd < 3; dd++) {
522 v2[dd] = diffN[3 * 0 + dd] * N[node_shift + 1] * N[node_shift + 2] *
523 N[node_shift + 3] +
524 N[node_shift + 0] * diffN[3 * 1 + dd] * N[node_shift + 2] *
525 N[node_shift + 3] +
526 N[node_shift + 0] * N[node_shift + 1] * diffN[3 * 2 + dd] *
527 N[node_shift + 3] +
528 N[node_shift + 0] * N[node_shift + 1] * N[node_shift + 2] *
529 diffN[3 * 3 + dd];
530 }
531 int shift = ii * P;
532 int jj = 0;
533 int oo = 0;
534 for (; oo <= (p - 4); oo++) {
535 int pp0 = 0;
536 for (; pp0 <= oo; pp0++) {
537 int pp1 = 0;
538 for (; (pp0 + pp1) <= oo; pp1++) {
539 int pp2 = oo - pp0 - pp1;
540 if (pp2 >= 0) {
541 if (volumeN != NULL) {
542 volumeN[shift + jj] = L0[pp0] * L1[pp1] * L2[pp2] * v;
543 }
544 if (diff_volumeN != NULL) {
545 dd = 0;
546 for (; dd < 3; dd++) {
547 diff_volumeN[3 * shift + 3 * jj + dd] =
548 (diffL0[dd * (p + 1) + pp0] * L1[pp1] * L2[pp2] +
549 L0[pp0] * diffL1[dd * (p + 1) + pp1] * L2[pp2] +
550 L0[pp0] * L1[pp1] * diffL2[dd * (p + 1) + pp2]) *
551 v;
552 diff_volumeN[3 * shift + 3 * jj + dd] +=
553 L0[pp0] * L1[pp1] * L2[pp2] * v2[dd];
554 }
555 }
556 jj++;
557 }
558 }
559 }
560 }
561 if (jj != P)
562 SETERRQ1(PETSC_COMM_SELF, 1, "wrong order %d", jj);
563 }
565}

◆ 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
diffNderivatives of barycentric coordinates, i.e. derivatives 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}
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
static PetscErrorCode ierr
Definition: l2.c:27

◆ 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
diffNderivatives of barycentric coordinates, i.e. derivatives 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}
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.