v0.8.16
Public Member Functions | Private Member Functions | Private Attributes | List of all members
MoFEM::TetPolynomialBase Struct Reference

Calculate base functions on tetrahedral. More...

#include <src/approximation/TetPolynomialBase.hpp>

Inheritance diagram for MoFEM::TetPolynomialBase:
[legend]
Collaboration diagram for MoFEM::TetPolynomialBase:
[legend]

Public Member Functions

MoFEMErrorCode query_interface (const MOFEMuuid &uuid, UnknownInterface **iface) const
 
 TetPolynomialBase ()
 
 ~TetPolynomialBase ()
 
MoFEMErrorCode getValue (MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
- Public Member Functions inherited from MoFEM::BaseFunction
 BaseFunction ()
 
 ~BaseFunction ()
 
virtual MoFEMErrorCode getValue (MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (const MOFEMuuid &uuid, bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE , bool VERIFY = false>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, IFACE *&iface) const
 Get interface by uuid and return reference to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()
 
virtual MoFEMErrorCode getLibVersion (Version &version) const
 Get library version. More...
 
virtual const MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version) const
 Get database major version. More...
 
virtual MoFEMErrorCode getInterfaceVersion (Version &version) const
 Get database major version. More...
 
template<>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, UnknownInterface *&iface) const
 

Private Member Functions

MoFEMErrorCode getValueH1 (MatrixDouble &pts)
 
MoFEMErrorCode getValueL2 (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdiv (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurl (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivAinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurlAinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivDemkowiczBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurlDemkowiczBase (MatrixDouble &pts)
 

Private Attributes

EntPolynomialBaseCtxcTx
 
ublas::matrix< MatrixDoubleN_face_edge
 
ublas::vector< MatrixDoubleN_face_bubble
 
ublas::vector< MatrixDoubleN_volume_edge
 
ublas::vector< MatrixDoubleN_volume_face
 
MatrixDouble N_volume_bubble
 
ublas::matrix< MatrixDoublediffN_face_edge
 
ublas::vector< MatrixDoublediffN_face_bubble
 
ublas::vector< MatrixDoublediffN_volume_edge
 
ublas::vector< MatrixDoublediffN_volume_face
 
MatrixDouble diffN_volume_bubble
 

Additional Inherited Members

- Protected Member Functions inherited from MoFEM::UnknownInterface
boost::typeindex::type_index getClassIdx (const MOFEMuuid &uid) const
 Get type name for interface Id. More...
 
MOFEMuuid getUId (const boost::typeindex::type_index &class_idx) const
 Get interface Id for class name. More...
 

Detailed Description

Calculate base functions on tetrahedral.

Definition at line 31 of file TetPolynomialBase.hpp.

Constructor & Destructor Documentation

◆ TetPolynomialBase()

TetPolynomialBase::TetPolynomialBase ( )

Definition at line 71 of file TetPolynomialBase.cpp.

71 {}

◆ ~TetPolynomialBase()

TetPolynomialBase::~TetPolynomialBase ( )

Definition at line 70 of file TetPolynomialBase.cpp.

70 {}

Member Function Documentation

◆ getValue()

MoFEMErrorCode TetPolynomialBase::getValue ( MatrixDouble pts,
boost::shared_ptr< BaseFunctionCtx ctx_ptr 
)
virtual

Reimplemented from MoFEM::BaseFunction.

Definition at line 907 of file TetPolynomialBase.cpp.

908  {
910 
912  CHKERR ctx_ptr->query_interface(IDD_TET_BASE_FUNCTION, &iface);
913  cTx = reinterpret_cast<EntPolynomialBaseCtx *>(iface);
914 
915  int nb_gauss_pts = pts.size2();
916  if (!nb_gauss_pts) {
918  }
919 
920  if (pts.size1() < 3) {
921  SETERRQ(
922  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
923  "Wrong dimension of pts, should be at least 3 rows with coordinates");
924  }
925 
926  const FieldApproximationBase base = cTx->bAse;
928  if (cTx->copyNodeBase == LASTBASE) {
929  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4, false);
931  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
932  &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
933  } else {
934  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
935  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
936  }
937  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
938  (unsigned int)nb_gauss_pts) {
939  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
940  "Base functions or nodes has wrong number of integration points "
941  "for base %s",
942  ApproximationBaseNames[base]);
943  }
944  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
946  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
947 
948  switch (cTx->sPace) {
949  case H1:
950  CHKERR getValueH1(pts);
951  break;
952  case HDIV:
953  CHKERR getValueHdiv(pts);
954  break;
955  case HCURL:
956  CHKERR getValueHcurl(pts);
957  break;
958  case L2:
959  CHKERR getValueL2(pts);
960  break;
961  default:
962  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space");
963  }
964 
966 }
field with continuous normal traction
Definition: definitions.h:170
data structure for finite element entityIt keeps that about indices of degrees of freedom...
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:291
Class used to pass element data to calculate base functions on tet,triangle,edge. ...
const FieldApproximationBase bAse
base class for all interface classes
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
MoFEMErrorCode getValueL2(MatrixDouble &pts)
static const MOFEMuuid IDD_TET_BASE_FUNCTION
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
const FieldApproximationBase copyNodeBase
static const char *const ApproximationBaseNames[]
Definition: definitions.h:155
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:303
FieldApproximationBase
approximation base
Definition: definitions.h:140
DataForcesAndSourcesCore & dAta
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
field with continuous tangents
Definition: definitions.h:169
MoFEMErrorCode getValueH1(MatrixDouble &pts)
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
continuous field
Definition: definitions.h:168
EntPolynomialBaseCtx * cTx
field with C-1 continuity
Definition: definitions.h:171

◆ getValueH1()

MoFEMErrorCode TetPolynomialBase::getValueH1 ( MatrixDouble pts)
private

Definition at line 73 of file TetPolynomialBase.cpp.

73  {
75 
77  const FieldApproximationBase base = cTx->bAse;
78  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
79  double *diffL, const int dim) =
81 
82  int nb_gauss_pts = pts.size2();
83 
84  int sense[6], order[6];
85  if (data.spacesOnEntities[MBEDGE].test(H1)) {
86  // edges
87  if (data.dataOnEntities[MBEDGE].size() != 6) {
88  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
89  }
90  double *h1_edge_n[6], *diff_h1_egde_n[6];
91  for (int ee = 0; ee != 6; ++ee) {
92  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
93  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
94  "data inconsistency");
95  }
96  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
97  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
98  int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getDataOrder());
99  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
100  false);
101  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
102  3 * nb_dofs, false);
103  h1_edge_n[ee] =
104  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
105  diff_h1_egde_n[ee] =
106  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
107  }
109  sense, order,
110  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
111  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
112  h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
113  } else {
114  for (int ee = 0; ee != 6; ++ee) {
115  data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
116  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
117  }
118  }
119 
120  if (data.spacesOnEntities[MBTRI].test(H1)) {
121  // faces
122  if (data.dataOnEntities[MBTRI].size() != 4) {
123  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
124  }
125  double *h1_face_n[4], *diff_h1_face_n[4];
126  for (int ff = 0; ff != 4; ++ff) {
127  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
128  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
129  "data inconsistency");
130  }
131  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getDataOrder());
132  order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
133  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
134  false);
135  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
136  3 * nb_dofs, false);
137  h1_face_n[ff] =
138  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
139  diff_h1_face_n[ff] =
140  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
141  }
142  if (data.facesNodes.size1() != 4) {
143  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
144  }
145  if (data.facesNodes.size2() != 3) {
146  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
147  }
149  &*data.facesNodes.data().begin(), order,
150  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
151  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
152  h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
153 
154  } else {
155  for (int ff = 0; ff != 4; ++ff) {
156  data.dataOnEntities[MBTRI][ff].getN(base).resize(0, false);
157  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, false);
158  }
159  }
160 
161  if (data.spacesOnEntities[MBTET].test(H1)) {
162  // volume
163  int order = data.dataOnEntities[MBTET][0].getDataOrder();
164  int nb_vol_dofs = NBVOLUMETET_H1(order);
165  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
166  false);
167  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
168  3 * nb_vol_dofs, false);
170  data.dataOnEntities[MBTET][0].getDataOrder(),
171  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
172  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
173  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
174  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
175  nb_gauss_pts, base_polynomials);
176  } else {
177  data.dataOnEntities[MBTET][0].getN(base).resize(0, 0, false);
178  data.dataOnEntities[MBTET][0].getDiffN(base).resize(0, 0, false);
179  }
180 
182 }
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))
Definition: h1.c:317
std::bitset< LASTSPACE > spacesOnEntities[MBMAXTYPE]
spaces on entity types
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))
Definition: h1.c:419
ublas::matrix< int > facesNodes
nodes on finite element faces
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
data structure for finite element entityIt keeps that about indices of degrees of freedom...
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))
Definition: h1.c:245
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron fro H1 space.
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
FieldApproximationBase
approximation base
Definition: definitions.h:140
DataForcesAndSourcesCore & dAta
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
continuous field
Definition: definitions.h:168
EntPolynomialBaseCtx * cTx

◆ getValueHcurl()

MoFEMErrorCode TetPolynomialBase::getValueHcurl ( MatrixDouble pts)
private

Definition at line 888 of file TetPolynomialBase.cpp.

888  {
890 
891  switch (cTx->bAse) {
895  break;
898  break;
899  default:
900  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
901  }
902 
904 }
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:142
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
EntPolynomialBaseCtx * cTx
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)

◆ getValueHcurlAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlAinsworthBase ( MatrixDouble pts)
private

Definition at line 638 of file TetPolynomialBase.cpp.

638  {
640 
642  const FieldApproximationBase base = cTx->bAse;
643  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
644  double *diffL, const int dim) =
646 
647  int nb_gauss_pts = pts.size2();
648 
649  // edges
650  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
651  int sense[6], order[6];
652  if (data.dataOnEntities[MBEDGE].size() != 6) {
653  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
654  }
655  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
656  for (int ee = 0; ee != 6; ee++) {
657  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
658  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
659  "data inconsistency");
660  }
661  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
662  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
663  int nb_dofs = NBEDGE_AINSWORTH_HCURL(
664  data.dataOnEntities[MBEDGE][ee].getDataOrder());
665  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
666  3 * nb_dofs, false);
667  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
668  9 * nb_dofs, false);
669  hcurl_edge_n[ee] =
670  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
671  diff_hcurl_edge_n[ee] =
672  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
673  }
675  sense, order,
676  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
677  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
678  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
679  } else {
680  for (int ee = 0; ee != 6; ee++) {
681  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
682  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
683  false);
684  }
685  }
686 
687  // triangles
688  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
689  int order[4];
690  // faces
691  if (data.dataOnEntities[MBTRI].size() != 4) {
692  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
693  }
694  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
695  for (int ff = 0; ff != 4; ff++) {
696  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
697  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
698  "data inconsistency");
699  }
700  order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
701  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
702  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
703  3 * nb_dofs, false);
704  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
705  9 * nb_dofs, false);
706  hcurl_base_n[ff] =
707  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
708  diff_hcurl_base_n[ff] =
709  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
710  }
711  if (data.facesNodes.size1() != 4) {
712  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
713  }
714  if (data.facesNodes.size2() != 3) {
715  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
716  }
718  &*data.facesNodes.data().begin(), order,
719  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
720  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
721  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
722  } else {
723  for (int ff = 0; ff != 4; ff++) {
724  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
725  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
726  false);
727  }
728  }
729 
730  if (data.spacesOnEntities[MBTET].test(HCURL)) {
731 
732  // volume
733  int order = data.dataOnEntities[MBTET][0].getDataOrder();
734  int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
735  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
736  3 * nb_vol_dofs, false);
737  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
738  9 * nb_vol_dofs, false);
740  data.dataOnEntities[MBTET][0].getDataOrder(),
741  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
742  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
743  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
744  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
745  nb_gauss_pts, base_polynomials);
746 
747  } else {
748  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
749  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
750  }
751 
753 }
std::bitset< LASTSPACE > spacesOnEntities[MBMAXTYPE]
spaces on entity types
ublas::matrix< int > facesNodes
nodes on finite element faces
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
data structure for finite element entityIt keeps that about indices of degrees of freedom...
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
FieldApproximationBase
approximation base
Definition: definitions.h:140
DataForcesAndSourcesCore & dAta
#define NBEDGE_AINSWORTH_HCURL(P)
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET(int *face_nodes, int *p, double *N, double *diffN, double *phi_f[4], double *diff_phi_f[4], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face H-curl functions.
Definition: Hcurl.cpp:1059
field with continuous tangents
Definition: definitions.h:169
#define CHKERR
Inline error check.
Definition: definitions.h:578
MoFEMErrorCode Hcurl_Ainsworth_VolumeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H-curl volume base functions.
Definition: Hcurl.cpp:1413
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Edge based H-curl base functions on tetrahedral.
Definition: Hcurl.cpp:22
#define NBVOLUMETET_AINSWORTH_HCURL(P)
#define NBFACETRI_AINSWORTH_HCURL(P)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
EntPolynomialBaseCtx * cTx

◆ getValueHcurlDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlDemkowiczBase ( MatrixDouble pts)
private

Definition at line 756 of file TetPolynomialBase.cpp.

756  {
758 
760  const FieldApproximationBase base = cTx->bAse;
761  if (base != DEMKOWICZ_JACOBI_BASE) {
762  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
763  "This should be used only with DEMKOWICZ_JACOBI_BASE "
764  "but base is %s",
765  ApproximationBaseNames[base]);
766  }
767 
768  int nb_gauss_pts = pts.size2();
769 
770  // edges
771  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
772  int sense[6], order[6];
773  if (data.dataOnEntities[MBEDGE].size() != 6) {
774  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
775  "wrong size of data structure, expected space for six edges "
776  "but is %d",
777  data.dataOnEntities[MBEDGE].size());
778  }
779  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
780  for (int ee = 0; ee != 6; ee++) {
781  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
782  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
783  "orintation of edges is not set");
784  }
785  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
786  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
787  int nb_dofs = NBEDGE_DEMKOWICZ_HCURL(
788  data.dataOnEntities[MBEDGE][ee].getDataOrder());
789  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
790  3 * nb_dofs, false);
791  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
792  9 * nb_dofs, false);
793  hcurl_edge_n[ee] =
794  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
795  diff_hcurl_edge_n[ee] =
796  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
797  }
799  sense, order,
800  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
801  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
802  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
803  } else {
804  // No DOFs on edges, resize base function matrices, indicating that no
805  // dofs on them.
806  for (int ee = 0; ee != 6; ee++) {
807  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
808  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
809  false);
810  }
811  }
812 
813  // triangles
814  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
815  int order[4];
816  // faces
817  if (data.dataOnEntities[MBTRI].size() != 4) {
818  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
819  "data structure for storing face h-curl base have wrong size "
820  "should be four but is %d",
821  data.dataOnEntities[MBTRI].size());
822  }
823  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
824  for (int ff = 0; ff != 4; ff++) {
825  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
826  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
827  "orintation of face is not set");
828  }
829  order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
830  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order[ff]);
831  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
832  3 * nb_dofs, false);
833  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
834  9 * nb_dofs, false);
835  hcurl_base_n[ff] =
836  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
837  diff_hcurl_base_n[ff] =
838  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
839  }
840  if (data.facesNodes.size1() != 4) {
841  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
842  "data inconsistency, should be four faces");
843  }
844  if (data.facesNodes.size2() != 3) {
845  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
846  "data inconsistency, should be three nodes on face");
847  }
849  &*data.facesNodes.data().begin(), order,
850  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
851  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
852  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
853  } else {
854  // No DOFs on faces, resize base function matrices, indicating that no
855  // dofs on them.
856  for (int ff = 0; ff != 4; ff++) {
857  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
858  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
859  false);
860  }
861  }
862 
863  if (data.spacesOnEntities[MBTET].test(HCURL)) {
864  // volume
865  int order = data.dataOnEntities[MBTET][0].getDataOrder();
866  int nb_vol_dofs = NBVOLUMETET_DEMKOWICZ_HCURL(order);
867  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
868  3 * nb_vol_dofs, false);
869  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
870  9 * nb_vol_dofs, false);
872  data.dataOnEntities[MBTET][0].getDataOrder(),
873  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
874  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
875  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
876  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
877  nb_gauss_pts);
878  } else {
879  // No DOFs on faces, resize base function matrices, indicating that no
880  // dofs on them.
881  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
882  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
883  }
884 
886 }
std::bitset< LASTSPACE > spacesOnEntities[MBMAXTYPE]
spaces on entity types
ublas::matrix< int > facesNodes
nodes on finite element faces
data structure for finite element entityIt keeps that about indices of degrees of freedom...
MoFEMErrorCode Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Volume base interior function.
Definition: Hcurl.cpp:2812
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
static const char *const ApproximationBaseNames[]
Definition: definitions.h:155
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTET(int *faces_nodes, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Face base interior function.
Definition: Hcurl.cpp:2723
FieldApproximationBase
approximation base
Definition: definitions.h:140
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(int *sense, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Edge based H-curl base functions on tetrahedral.
Definition: Hcurl.cpp:2197
DataForcesAndSourcesCore & dAta
field with continuous tangents
Definition: definitions.h:169
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define NBFACETRI_DEMKOWICZ_HCURL(P)
#define NBEDGE_DEMKOWICZ_HCURL(P)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
EntPolynomialBaseCtx * cTx
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)

◆ getValueHdiv()

MoFEMErrorCode TetPolynomialBase::getValueHdiv ( MatrixDouble pts)
private

Definition at line 621 of file TetPolynomialBase.cpp.

621  {
623 
624  switch (cTx->bAse) {
627  return getValueHdivAinsworthBase(pts);
629  return getValueHdivDemkowiczBase(pts);
630  default:
631  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
632  }
633 
635 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
const FieldApproximationBase bAse
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:142
EntPolynomialBaseCtx * cTx

◆ getValueHdivAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBase ( MatrixDouble pts)
private

Definition at line 213 of file TetPolynomialBase.cpp.

213  {
215 
217  const FieldApproximationBase base = cTx->bAse;
218  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
219  double *diffL, const int dim) =
221 
222  int nb_gauss_pts = pts.size2();
223 
224  // face shape functions
225 
226  double *phi_f_e[4][3];
227  double *phi_f[4];
228  double *diff_phi_f_e[4][3];
229  double *diff_phi_f[4];
230 
231  N_face_edge.resize(4, 3, false);
232  N_face_bubble.resize(4, false);
233  diffN_face_edge.resize(4, 3, false);
234  diffN_face_bubble.resize(4, false);
235 
236  int faces_order[4];
237  for (int ff = 0; ff != 4; ++ff) {
238  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
239  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
240  }
241  faces_order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
242  // three edges on face
243  for (int ee = 0; ee < 3; ee++) {
244  N_face_edge(ff, ee).resize(
245  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
246  false);
247  diffN_face_edge(ff, ee).resize(
248  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
249  false);
250  phi_f_e[ff][ee] = &((N_face_edge(ff, ee))(0, 0));
251  diff_phi_f_e[ff][ee] = &((diffN_face_edge(ff, ee))(0, 0));
252  }
253  N_face_bubble[ff].resize(nb_gauss_pts,
254  3 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
255  false);
256  diffN_face_bubble[ff].resize(
257  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
258  false);
259  phi_f[ff] = &*(N_face_bubble[ff].data().begin());
260  diff_phi_f[ff] = &*(diffN_face_bubble[ff].data().begin());
261  }
262 
264  &data.facesNodes(0, 0), faces_order,
265  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
266  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_f_e,
267  diff_phi_f_e, nb_gauss_pts, base_polynomials);
268 
270  &data.facesNodes(0, 0), faces_order,
271  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
272  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_f, diff_phi_f,
273  nb_gauss_pts, base_polynomials);
274 
275  // volume shape functions
276 
277  double *phi_v_e[6];
278  double *phi_v_f[4];
279  double *phi_v;
280  double *diff_phi_v_e[6];
281  double *diff_phi_v_f[4];
282  double *diff_phi_v;
283 
284  int volume_order = data.dataOnEntities[MBTET][0].getDataOrder();
285 
286  N_volume_edge.resize(6, false);
287  diffN_volume_edge.resize(6, false);
288  for (int ee = 0; ee != 6; ++ee) {
289  N_volume_edge[ee].resize(
290  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
291  diffN_volume_edge[ee].resize(
292  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
293  phi_v_e[ee] = &*(N_volume_edge[ee].data().begin());
294  diff_phi_v_e[ee] = &*(diffN_volume_edge[ee].data().begin());
295  }
297  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
298  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_e,
299  diff_phi_v_e, nb_gauss_pts, base_polynomials);
300 
301  N_volume_face.resize(4, false);
302  diffN_volume_face.resize(4, false);
303  for (int ff = 0; ff != 4; ++ff) {
304  N_volume_face[ff].resize(
305  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
306  diffN_volume_face[ff].resize(
307  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
308  phi_v_f[ff] = &*(N_volume_face[ff].data().begin());
309  diff_phi_v_f[ff] = &*(diffN_volume_face[ff].data().begin());
310  }
312  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
313  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_f,
314  diff_phi_v_f, nb_gauss_pts, base_polynomials);
315 
316  N_volume_bubble.resize(
317  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
318  diffN_volume_bubble.resize(
319  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
320  phi_v = &*(N_volume_bubble.data().begin());
321  diff_phi_v = &*(diffN_volume_bubble.data().begin());
323  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
324  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v, diff_phi_v,
325  nb_gauss_pts, base_polynomials);
326 
327  // Set shape functions into data structure Shape functions hast to be put
328  // in arrays in order which guarantee hierarchical series of degrees of
329  // freedom, i.e. in other words dofs form sub-entities has to be group
330  // by order.
331 
334 
335  // faces
336  if (data.dataOnEntities[MBTRI].size() != 4) {
337  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
338  }
339  for (int ff = 0; ff != 4; ff++) {
340  data.dataOnEntities[MBTRI][ff].getN(base).resize(
341  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
342  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
343  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
344  if (NBFACETRI_AINSWORTH_HDIV(faces_order[ff]) == 0)
345  continue;
346  // face
347  double *base_ptr =
348  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
349  FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
350  &base_ptr[HVEC2], 3);
351  double *diff_base_ptr =
352  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
354  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
355  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
356  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
357  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
358  &diff_base_ptr[HVEC2_2], 9);
359  // face-face
360  boost::shared_ptr<FTensor::Tensor1<double *, 3> > t_base_f;
361  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> > t_diff_base_f;
362  if (NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]) > 0) {
363  base_ptr = phi_f[ff];
364  t_base_f = boost::shared_ptr<FTensor::Tensor1<double *, 3> >(
365  new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
366  &base_ptr[HVEC2], 3));
367  diff_base_ptr = diff_phi_f[ff];
368  t_diff_base_f = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> >(
370  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
371  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
372  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
373  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
374  &diff_base_ptr[HVEC2_2], 9));
375  }
376  // edge-face
377  base_ptr = phi_f_e[ff][0];
378  FTensor::Tensor1<double *, 3> t_base_f_e0(base_ptr, &base_ptr[HVEC1],
379  &base_ptr[HVEC2], 3);
380  diff_base_ptr = diff_phi_f_e[ff][0];
381  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e0(
382  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
383  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
384  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
385  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
386  &diff_base_ptr[HVEC2_2], 9);
387  base_ptr = phi_f_e[ff][1];
388  FTensor::Tensor1<double *, 3> t_base_f_e1(base_ptr, &base_ptr[HVEC1],
389  &base_ptr[HVEC2], 3);
390  diff_base_ptr = diff_phi_f_e[ff][1];
391  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e1(
392  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
393  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
394  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
395  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
396  &diff_base_ptr[HVEC2_2], 9);
397  base_ptr = phi_f_e[ff][2];
398  FTensor::Tensor1<double *, 3> t_base_f_e2(base_ptr, &base_ptr[HVEC1],
399  &base_ptr[HVEC2], 3);
400  diff_base_ptr = diff_phi_f_e[ff][2];
401  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e2(
402  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
403  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
404  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
405  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
406  &diff_base_ptr[HVEC2_2], 9);
407  for (int gg = 0; gg != nb_gauss_pts; gg++) {
408  for (int oo = 0; oo != faces_order[ff]; oo++) {
409  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
410  dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
411  t_base(i) = t_base_f_e0(i);
412  ++t_base;
413  ++t_base_f_e0;
414  t_diff_base(i, j) = t_diff_base_f_e0(i, j);
415  ++t_diff_base;
416  ++t_diff_base_f_e0;
417  t_base(i) = t_base_f_e1(i);
418  ++t_base;
419  ++t_base_f_e1;
420  t_diff_base(i, j) = t_diff_base_f_e1(i, j);
421  ++t_diff_base;
422  ++t_diff_base_f_e1;
423  t_base(i) = t_base_f_e2(i);
424  ++t_base;
425  ++t_base_f_e2;
426  t_diff_base(i, j) = t_diff_base_f_e2(i, j);
427  ++t_diff_base;
428  ++t_diff_base_f_e2;
429  }
430  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
431  dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
432  t_base(i) = (*t_base_f)(i);
433  ++t_base;
434  ++(*t_base_f);
435  t_diff_base(i, j) = (*t_diff_base_f)(i, j);
436  ++t_diff_base;
437  ++(*t_diff_base_f);
438  }
439  }
440  }
441  }
442 
443  // volume
444  data.dataOnEntities[MBTET][0].getN(base).resize(
445  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
446  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
447  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
448  if (NBVOLUMETET_AINSWORTH_HDIV(volume_order) > 0) {
449  double *base_ptr =
450  &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
451  FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
452  &base_ptr[HVEC2], 3);
453  double *diff_base_ptr =
454  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
456  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
457  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
458  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
459  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
460  &diff_base_ptr[HVEC2_2], 9);
461  // edges
462  std::vector<FTensor::Tensor1<double *, 3> > t_base_v_e;
463  t_base_v_e.reserve(6);
464  std::vector<FTensor::Tensor2<double *, 3, 3> > t_diff_base_v_e;
465  t_diff_base_v_e.reserve(6);
466  for (int ee = 0; ee != 6; ee++) {
467  base_ptr = phi_v_e[ee];
468  diff_base_ptr = diff_phi_v_e[ee];
469  t_base_v_e.push_back(FTensor::Tensor1<double *, 3>(
470  base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
471  t_diff_base_v_e.push_back(FTensor::Tensor2<double *, 3, 3>(
472  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
473  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
474  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
475  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
476  &diff_base_ptr[HVEC2_2], 9));
477  }
478  // faces
479  std::vector<FTensor::Tensor1<double *, 3> > t_base_v_f;
480  t_base_v_f.reserve(4);
481  std::vector<FTensor::Tensor2<double *, 3, 3> > t_diff_base_v_f;
482  t_diff_base_v_f.reserve(4);
483  if (NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order) > 0) {
484  for (int ff = 0; ff != 4; ff++) {
485  base_ptr = phi_v_f[ff];
486  diff_base_ptr = diff_phi_v_f[ff];
487  t_base_v_f.push_back(FTensor::Tensor1<double *, 3>(
488  base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
489  t_diff_base_v_f.push_back(FTensor::Tensor2<double *, 3, 3>(
490  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
491  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
492  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
493  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
494  &diff_base_ptr[HVEC2_2], 9));
495  }
496  }
497  boost::shared_ptr<FTensor::Tensor1<double *, 3> > t_base_v;
498  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> > t_diff_base_v;
499  if (NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order) > 0) {
500  base_ptr = phi_v;
501  t_base_v = boost::shared_ptr<FTensor::Tensor1<double *, 3> >(
502  new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
503  &base_ptr[HVEC2], 3));
504  diff_base_ptr = diff_phi_v;
505  t_diff_base_v = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> >(
507  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
508  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
509  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
510  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
511  &diff_base_ptr[HVEC2_2], 9));
512  }
513  for (int gg = 0; gg != nb_gauss_pts; gg++) {
514  for (int oo = 0; oo < volume_order; oo++) {
515  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
516  dd < NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
517  for (int ee = 0; ee < 6; ee++) {
518  t_base(i) = t_base_v_e[ee](i);
519  ++t_base;
520  ++t_base_v_e[ee];
521  t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
522  ++t_diff_base;
523  ++t_diff_base_v_e[ee];
524  }
525  }
526  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
527  dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
528  for (int ff = 0; ff < 4; ff++) {
529  t_base(i) = t_base_v_f[ff](i);
530  ++t_base;
531  ++t_base_v_f[ff];
532  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
533  ++t_diff_base;
534  ++t_diff_base_v_f[ff];
535  }
536  }
537  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
538  dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
539  t_base(i) = (*t_base_v)(i);
540  ++t_base;
541  ++(*t_base_v);
542  t_diff_base(i, j) = (*t_diff_base_v)(i, j);
543  ++t_diff_base;
544  ++(*t_diff_base_v);
545  }
546  }
547  }
548  }
549 
551 }
ublas::matrix< int > facesNodes
nodes on finite element faces
ublas::matrix< MatrixDouble > diffN_face_edge
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
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
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
data structure for finite element entityIt keeps that about indices of degrees of freedom...
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
ublas::vector< MatrixDouble > N_volume_face
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
MoFEMErrorCode Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Interior bubble functions by Ainsworth .
Definition: Hdiv.cpp:488
ublas::matrix< MatrixDouble > N_face_edge
MoFEMErrorCode Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: Hdiv.cpp:372
FieldApproximationBase
approximation base
Definition: definitions.h:140
DataForcesAndSourcesCore & dAta
ublas::vector< MatrixDouble > diffN_volume_edge
ublas::vector< MatrixDouble > diffN_face_bubble
ublas::vector< MatrixDouble > N_volume_edge
#define CHKERR
Inline error check.
Definition: definitions.h:578
ublas::vector< MatrixDouble > diffN_volume_face
MoFEMErrorCode Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_e[6], double *diff_phi_v_e[6], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base function, Edge-based interior (volume) functions by Ainsworth .
Definition: Hdiv.cpp:284
#define NBVOLUMETET_AINSWORTH_HDIV(P)
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3], double *diff_phi_f_e[4][3], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
Definition: Hdiv.cpp:19
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
#define NBFACETRI_AINSWORTH_HDIV(P)
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[], double *diff_phi_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
Definition: Hdiv.cpp:141
EntPolynomialBaseCtx * cTx
ublas::vector< MatrixDouble > N_face_bubble
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)

◆ getValueHdivDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBase ( MatrixDouble pts)
private

Definition at line 553 of file TetPolynomialBase.cpp.

553  {
555 
557  const FieldApproximationBase base = cTx->bAse;
558  if (base != DEMKOWICZ_JACOBI_BASE) {
559  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
560  "This should be used only with DEMKOWICZ_JACOBI_BASE "
561  "but base is %s",
562  ApproximationBaseNames[base]);
563  }
564  int nb_gauss_pts = pts.size2();
565 
566  int volume_order = data.dataOnEntities[MBTET][0].getDataOrder();
567 
568  int p_f[4];
569  double *phi_f[4];
570  double *diff_phi_f[4];
571 
572  // Calculate base function on tet faces
573  for (int ff = 0; ff != 4; ff++) {
574  int face_order = data.dataOnEntities[MBTRI][ff].getDataOrder();
575  int order = volume_order > face_order ? volume_order : face_order;
576  data.dataOnEntities[MBTRI][ff].getN(base).resize(
577  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
578  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
579  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
580  p_f[ff] = order;
581  phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
582  diff_phi_f[ff] =
583  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
584  if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
585  continue;
587  &data.facesNodes(ff, 0), order,
588  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
589  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
590  phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
591  }
592 
593  // Calculate base functions in tet interior
594  if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
595  data.dataOnEntities[MBTET][0].getN(base).resize(
596  nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
597  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
598  nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
599  double *phi_v =
600  &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
601  double *diff_phi_v =
602  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
604  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
605  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
606  diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
607  }
608 
609  // Set size of face base correctly
610  for (int ff = 0; ff != 4; ff++) {
611  int face_order = data.dataOnEntities[MBTRI][ff].getDataOrder();
612  data.dataOnEntities[MBTRI][ff].getN(base).resize(
613  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
614  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
615  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
616  }
617 
619 }
ublas::matrix< int > facesNodes
nodes on finite element faces
data structure for finite element entityIt keeps that about indices of degrees of freedom...
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition: Hdiv.cpp:621
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
static const char *const ApproximationBaseNames[]
Definition: definitions.h:155
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
Definition: Hdiv.cpp:766
#define NBFACETRI_DEMKOWICZ_HDIV(P)
FieldApproximationBase
approximation base
Definition: definitions.h:140
DataForcesAndSourcesCore & dAta
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
EntPolynomialBaseCtx * cTx

◆ getValueL2()

MoFEMErrorCode TetPolynomialBase::getValueL2 ( MatrixDouble pts)
private

Definition at line 184 of file TetPolynomialBase.cpp.

184  {
186 
188  const FieldApproximationBase base = cTx->bAse;
189  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
190  double *diffL, const int dim) =
192 
193  int nb_gauss_pts = pts.size2();
194 
195  data.dataOnEntities[MBTET][0].getN(base).resize(
196  nb_gauss_pts,
197  NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getDataOrder()), false);
198  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
199  nb_gauss_pts,
200  3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getDataOrder()), false);
201 
203  data.dataOnEntities[MBTET][0].getDataOrder(),
204  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
205  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
206  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
207  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
208  nb_gauss_pts, base_polynomials);
209 
211 }
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
data structure for finite element entityIt keeps that about indices of degrees of freedom...
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
FieldApproximationBase
approximation base
Definition: definitions.h:140
DataForcesAndSourcesCore & dAta
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403
EntPolynomialBaseCtx * cTx
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

◆ query_interface()

MoFEMErrorCode TetPolynomialBase::query_interface ( const MOFEMuuid uuid,
MoFEM::UnknownInterface **  iface 
) const
virtual

Reimplemented from MoFEM::BaseFunction.

Definition at line 54 of file TetPolynomialBase.cpp.

55  {
56 
58  *iface = NULL;
59  if (uuid == IDD_TET_BASE_FUNCTION) {
60  *iface = const_cast<TetPolynomialBase *>(this);
62  } else {
63  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
64  }
65  ierr = BaseFunction::query_interface(uuid, iface);
66  CHKERRG(ierr);
68 }
Calculate base functions on tetrahedral.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:526
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
static const MOFEMuuid IDD_TET_BASE_FUNCTION
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, MoFEM::UnknownInterface **iface) const

Member Data Documentation

◆ cTx

EntPolynomialBaseCtx* MoFEM::TetPolynomialBase::cTx
private

Definition at line 43 of file TetPolynomialBase.hpp.

◆ diffN_face_bubble

ublas::vector<MatrixDouble> MoFEM::TetPolynomialBase::diffN_face_bubble
private

Definition at line 56 of file TetPolynomialBase.hpp.

◆ diffN_face_edge

ublas::matrix<MatrixDouble> MoFEM::TetPolynomialBase::diffN_face_edge
private

Definition at line 55 of file TetPolynomialBase.hpp.

◆ diffN_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::diffN_volume_bubble
private

Definition at line 59 of file TetPolynomialBase.hpp.

◆ diffN_volume_edge

ublas::vector<MatrixDouble> MoFEM::TetPolynomialBase::diffN_volume_edge
private

Definition at line 57 of file TetPolynomialBase.hpp.

◆ diffN_volume_face

ublas::vector<MatrixDouble> MoFEM::TetPolynomialBase::diffN_volume_face
private

Definition at line 58 of file TetPolynomialBase.hpp.

◆ N_face_bubble

ublas::vector<MatrixDouble> MoFEM::TetPolynomialBase::N_face_bubble
private

Definition at line 50 of file TetPolynomialBase.hpp.

◆ N_face_edge

ublas::matrix<MatrixDouble> MoFEM::TetPolynomialBase::N_face_edge
private

Definition at line 49 of file TetPolynomialBase.hpp.

◆ N_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::N_volume_bubble
private

Definition at line 53 of file TetPolynomialBase.hpp.

◆ N_volume_edge

ublas::vector<MatrixDouble> MoFEM::TetPolynomialBase::N_volume_edge
private

Definition at line 51 of file TetPolynomialBase.hpp.

◆ N_volume_face

ublas::vector<MatrixDouble> MoFEM::TetPolynomialBase::N_volume_face
private

Definition at line 52 of file TetPolynomialBase.hpp.


The documentation for this struct was generated from the following files: