v0.8.23
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 43 of file TetPolynomialBase.cpp.

43 {}

◆ ~TetPolynomialBase()

TetPolynomialBase::~TetPolynomialBase ( )

Definition at line 42 of file TetPolynomialBase.cpp.

42 {}

Member Function Documentation

◆ getValue()

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

Reimplemented from MoFEM::BaseFunction.

Definition at line 879 of file TetPolynomialBase.cpp.

880  {
882 
884  CHKERR ctx_ptr->query_interface(IDD_TET_BASE_FUNCTION, &iface);
885  cTx = reinterpret_cast<EntPolynomialBaseCtx *>(iface);
886 
887  int nb_gauss_pts = pts.size2();
888  if (!nb_gauss_pts) {
890  }
891 
892  if (pts.size1() < 3) {
893  SETERRQ(
894  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
895  "Wrong dimension of pts, should be at least 3 rows with coordinates");
896  }
897 
898  const FieldApproximationBase base = cTx->bAse;
900  if (cTx->copyNodeBase == LASTBASE) {
901  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4, false);
903  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
904  &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
905  } else {
906  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
907  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
908  }
909  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
910  (unsigned int)nb_gauss_pts) {
911  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
912  "Base functions or nodes has wrong number of integration points "
913  "for base %s",
914  ApproximationBaseNames[base]);
915  }
916  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
917  std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
918  data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
919 
920  switch (cTx->sPace) {
921  case H1:
922  CHKERR getValueH1(pts);
923  break;
924  case HDIV:
925  CHKERR getValueHdiv(pts);
926  break;
927  case HCURL:
928  CHKERR getValueHcurl(pts);
929  break;
930  case L2:
931  CHKERR getValueL2(pts);
932  break;
933  default:
934  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space");
935  }
936 
938 }
field with continuous normal traction
Definition: definitions.h:173
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
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:477
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:120
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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:158
FieldApproximationBase
approximation base
Definition: definitions.h:143
DataForcesAndSourcesCore & dAta
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
field with continuous tangents
Definition: definitions.h:172
MoFEMErrorCode getValueH1(MatrixDouble &pts)
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
continuous field
Definition: definitions.h:171
EntPolynomialBaseCtx * cTx
static MoFEMErrorCode shapeFunMBTET(double *shape, const double *ksi, const double *eta, const double *zeta, const double nb)
Calculate shape functions on tetrahedron.
Definition: Tools.hpp:420
field with C-1 continuity
Definition: definitions.h:174

◆ getValueH1()

MoFEMErrorCode TetPolynomialBase::getValueH1 ( MatrixDouble pts)
private

Definition at line 45 of file TetPolynomialBase.cpp.

45  {
47 
49  const FieldApproximationBase base = cTx->bAse;
50  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
51  double *diffL, const int dim) =
53 
54  int nb_gauss_pts = pts.size2();
55 
56  int sense[6], order[6];
57  if (data.spacesOnEntities[MBEDGE].test(H1)) {
58  // edges
59  if (data.dataOnEntities[MBEDGE].size() != 6) {
60  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
61  }
62  double *h1_edge_n[6], *diff_h1_egde_n[6];
63  for (int ee = 0; ee != 6; ++ee) {
64  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
65  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
66  "data inconsistency");
67  }
68  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
69  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
70  int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getDataOrder());
71  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
72  false);
73  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
74  3 * nb_dofs, false);
75  h1_edge_n[ee] =
76  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
77  diff_h1_egde_n[ee] =
78  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
79  }
81  sense, order,
82  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
83  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
84  h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
85  } else {
86  for (int ee = 0; ee != 6; ++ee) {
87  data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
88  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
89  }
90  }
91 
92  if (data.spacesOnEntities[MBTRI].test(H1)) {
93  // faces
94  if (data.dataOnEntities[MBTRI].size() != 4) {
95  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
96  }
97  double *h1_face_n[4], *diff_h1_face_n[4];
98  for (int ff = 0; ff != 4; ++ff) {
99  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
100  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
101  "data inconsistency");
102  }
103  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getDataOrder());
104  order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
105  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
106  false);
107  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
108  3 * nb_dofs, false);
109  h1_face_n[ff] =
110  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
111  diff_h1_face_n[ff] =
112  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
113  }
114  if (data.facesNodes.size1() != 4) {
115  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
116  }
117  if (data.facesNodes.size2() != 3) {
118  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
119  }
121  &*data.facesNodes.data().begin(), order,
122  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
123  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
124  h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
125 
126  } else {
127  for (int ff = 0; ff != 4; ++ff) {
128  data.dataOnEntities[MBTRI][ff].getN(base).resize(0, false);
129  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, false);
130  }
131  }
132 
133  if (data.spacesOnEntities[MBTET].test(H1)) {
134  // volume
135  int order = data.dataOnEntities[MBTET][0].getDataOrder();
136  int nb_vol_dofs = NBVOLUMETET_H1(order);
137  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
138  false);
139  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
140  3 * nb_vol_dofs, false);
142  data.dataOnEntities[MBTET][0].getDataOrder(),
143  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
144  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
145  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
146  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
147  nb_gauss_pts, base_polynomials);
148  } else {
149  data.dataOnEntities[MBTET][0].getN(base).resize(0, 0, false);
150  data.dataOnEntities[MBTET][0].getDiffN(base).resize(0, 0, false);
151  }
152 
154 }
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 for H1 space.
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
FieldApproximationBase
approximation base
Definition: definitions.h:143
DataForcesAndSourcesCore & dAta
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
continuous field
Definition: definitions.h:171
EntPolynomialBaseCtx * cTx
const int order
approximation order

◆ getValueHcurl()

MoFEMErrorCode TetPolynomialBase::getValueHcurl ( MatrixDouble pts)
private

Definition at line 860 of file TetPolynomialBase.cpp.

860  {
862 
863  switch (cTx->bAse) {
867  break;
870  break;
871  default:
872  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
873  }
874 
876 }
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:477
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:145
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
EntPolynomialBaseCtx * cTx
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)

◆ getValueHcurlAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlAinsworthBase ( MatrixDouble pts)
private

Definition at line 610 of file TetPolynomialBase.cpp.

610  {
612 
614  const FieldApproximationBase base = cTx->bAse;
615  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
616  double *diffL, const int dim) =
618 
619  int nb_gauss_pts = pts.size2();
620 
621  // edges
622  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
623  int sense[6], order[6];
624  if (data.dataOnEntities[MBEDGE].size() != 6) {
625  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
626  }
627  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
628  for (int ee = 0; ee != 6; ee++) {
629  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
630  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
631  "data inconsistency");
632  }
633  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
634  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
635  int nb_dofs = NBEDGE_AINSWORTH_HCURL(
636  data.dataOnEntities[MBEDGE][ee].getDataOrder());
637  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
638  3 * nb_dofs, false);
639  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
640  9 * nb_dofs, false);
641  hcurl_edge_n[ee] =
642  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
643  diff_hcurl_edge_n[ee] =
644  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
645  }
647  sense, order,
648  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
649  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
650  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
651  } else {
652  for (int ee = 0; ee != 6; ee++) {
653  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
654  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
655  false);
656  }
657  }
658 
659  // triangles
660  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
661  int order[4];
662  // faces
663  if (data.dataOnEntities[MBTRI].size() != 4) {
664  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
665  }
666  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
667  for (int ff = 0; ff != 4; ff++) {
668  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
669  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
670  "data inconsistency");
671  }
672  order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
673  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
674  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
675  3 * nb_dofs, false);
676  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
677  9 * nb_dofs, false);
678  hcurl_base_n[ff] =
679  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
680  diff_hcurl_base_n[ff] =
681  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
682  }
683  if (data.facesNodes.size1() != 4) {
684  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
685  }
686  if (data.facesNodes.size2() != 3) {
687  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
688  }
690  &*data.facesNodes.data().begin(), order,
691  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
692  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
693  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
694  } else {
695  for (int ff = 0; ff != 4; ff++) {
696  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
697  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
698  false);
699  }
700  }
701 
702  if (data.spacesOnEntities[MBTET].test(HCURL)) {
703 
704  // volume
705  int order = data.dataOnEntities[MBTET][0].getDataOrder();
706  int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
707  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
708  3 * nb_vol_dofs, false);
709  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
710  9 * nb_vol_dofs, false);
712  data.dataOnEntities[MBTET][0].getDataOrder(),
713  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
714  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
715  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
716  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
717  nb_gauss_pts, base_polynomials);
718 
719  } else {
720  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
721  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
722  }
723 
725 }
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:477
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
FieldApproximationBase
approximation base
Definition: definitions.h:143
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:1052
field with continuous tangents
Definition: definitions.h:172
#define CHKERR
Inline error check.
Definition: definitions.h:596
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:1403
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:16
#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:407
EntPolynomialBaseCtx * cTx
const int order
approximation order

◆ getValueHcurlDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlDemkowiczBase ( MatrixDouble pts)
private

Definition at line 728 of file TetPolynomialBase.cpp.

728  {
730 
732  const FieldApproximationBase base = cTx->bAse;
733  if (base != DEMKOWICZ_JACOBI_BASE) {
734  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
735  "This should be used only with DEMKOWICZ_JACOBI_BASE "
736  "but base is %s",
737  ApproximationBaseNames[base]);
738  }
739 
740  int nb_gauss_pts = pts.size2();
741 
742  // edges
743  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
744  int sense[6], order[6];
745  if (data.dataOnEntities[MBEDGE].size() != 6) {
746  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
747  "wrong size of data structure, expected space for six edges "
748  "but is %d",
749  data.dataOnEntities[MBEDGE].size());
750  }
751  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
752  for (int ee = 0; ee != 6; ee++) {
753  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
754  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
755  "orintation of edges is not set");
756  }
757  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
758  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
759  int nb_dofs = NBEDGE_DEMKOWICZ_HCURL(
760  data.dataOnEntities[MBEDGE][ee].getDataOrder());
761  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
762  3 * nb_dofs, false);
763  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
764  9 * nb_dofs, false);
765  hcurl_edge_n[ee] =
766  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
767  diff_hcurl_edge_n[ee] =
768  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
769  }
771  sense, order,
772  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
773  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
774  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
775  } else {
776  // No DOFs on edges, resize base function matrices, indicating that no
777  // dofs on them.
778  for (int ee = 0; ee != 6; ee++) {
779  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
780  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
781  false);
782  }
783  }
784 
785  // triangles
786  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
787  int order[4];
788  // faces
789  if (data.dataOnEntities[MBTRI].size() != 4) {
790  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
791  "data structure for storing face h-curl base have wrong size "
792  "should be four but is %d",
793  data.dataOnEntities[MBTRI].size());
794  }
795  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
796  for (int ff = 0; ff != 4; ff++) {
797  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
798  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
799  "orintation of face is not set");
800  }
801  order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
802  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order[ff]);
803  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
804  3 * nb_dofs, false);
805  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
806  9 * nb_dofs, false);
807  hcurl_base_n[ff] =
808  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
809  diff_hcurl_base_n[ff] =
810  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
811  }
812  if (data.facesNodes.size1() != 4) {
813  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
814  "data inconsistency, should be four faces");
815  }
816  if (data.facesNodes.size2() != 3) {
817  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
818  "data inconsistency, should be three nodes on face");
819  }
821  &*data.facesNodes.data().begin(), order,
822  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
823  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
824  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
825  } else {
826  // No DOFs on faces, resize base function matrices, indicating that no
827  // dofs on them.
828  for (int ff = 0; ff != 4; ff++) {
829  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
830  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
831  false);
832  }
833  }
834 
835  if (data.spacesOnEntities[MBTET].test(HCURL)) {
836  // volume
837  int order = data.dataOnEntities[MBTET][0].getDataOrder();
838  int nb_vol_dofs = NBVOLUMETET_DEMKOWICZ_HCURL(order);
839  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
840  3 * nb_vol_dofs, false);
841  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
842  9 * nb_vol_dofs, false);
844  data.dataOnEntities[MBTET][0].getDataOrder(),
845  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
846  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
847  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
848  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
849  nb_gauss_pts);
850  } else {
851  // No DOFs on faces, resize base function matrices, indicating that no
852  // dofs on them.
853  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
854  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
855  }
856 
858 }
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:2467
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
static const char *const ApproximationBaseNames[]
Definition: definitions.h:158
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:2394
FieldApproximationBase
approximation base
Definition: definitions.h:143
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:2070
DataForcesAndSourcesCore & dAta
field with continuous tangents
Definition: definitions.h:172
#define CHKERR
Inline error check.
Definition: definitions.h:596
#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:407
EntPolynomialBaseCtx * cTx
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
const int order
approximation order

◆ getValueHdiv()

MoFEMErrorCode TetPolynomialBase::getValueHdiv ( MatrixDouble pts)
private

Definition at line 593 of file TetPolynomialBase.cpp.

593  {
595 
596  switch (cTx->bAse) {
599  return getValueHdivAinsworthBase(pts);
601  return getValueHdivDemkowiczBase(pts);
602  default:
603  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
604  }
605 
607 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
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:508
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:145
EntPolynomialBaseCtx * cTx

◆ getValueHdivAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBase ( MatrixDouble pts)
private

Definition at line 185 of file TetPolynomialBase.cpp.

185  {
187 
189  const FieldApproximationBase base = cTx->bAse;
190  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
191  double *diffL, const int dim) =
193 
194  int nb_gauss_pts = pts.size2();
195 
196  // face shape functions
197 
198  double *phi_f_e[4][3];
199  double *phi_f[4];
200  double *diff_phi_f_e[4][3];
201  double *diff_phi_f[4];
202 
203  N_face_edge.resize(4, 3, false);
204  N_face_bubble.resize(4, false);
205  diffN_face_edge.resize(4, 3, false);
206  diffN_face_bubble.resize(4, false);
207 
208  int faces_order[4];
209  for (int ff = 0; ff != 4; ++ff) {
210  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
211  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
212  }
213  faces_order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
214  // three edges on face
215  for (int ee = 0; ee < 3; ee++) {
216  N_face_edge(ff, ee).resize(
217  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
218  false);
219  diffN_face_edge(ff, ee).resize(
220  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
221  false);
222  phi_f_e[ff][ee] = &((N_face_edge(ff, ee))(0, 0));
223  diff_phi_f_e[ff][ee] = &((diffN_face_edge(ff, ee))(0, 0));
224  }
225  N_face_bubble[ff].resize(nb_gauss_pts,
226  3 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
227  false);
228  diffN_face_bubble[ff].resize(
229  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
230  false);
231  phi_f[ff] = &*(N_face_bubble[ff].data().begin());
232  diff_phi_f[ff] = &*(diffN_face_bubble[ff].data().begin());
233  }
234 
236  &data.facesNodes(0, 0), faces_order,
237  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
238  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_f_e,
239  diff_phi_f_e, nb_gauss_pts, base_polynomials);
240 
242  &data.facesNodes(0, 0), faces_order,
243  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
244  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_f, diff_phi_f,
245  nb_gauss_pts, base_polynomials);
246 
247  // volume shape functions
248 
249  double *phi_v_e[6];
250  double *phi_v_f[4];
251  double *phi_v;
252  double *diff_phi_v_e[6];
253  double *diff_phi_v_f[4];
254  double *diff_phi_v;
255 
256  int volume_order = data.dataOnEntities[MBTET][0].getDataOrder();
257 
258  N_volume_edge.resize(6, false);
259  diffN_volume_edge.resize(6, false);
260  for (int ee = 0; ee != 6; ++ee) {
261  N_volume_edge[ee].resize(
262  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
263  diffN_volume_edge[ee].resize(
264  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
265  phi_v_e[ee] = &*(N_volume_edge[ee].data().begin());
266  diff_phi_v_e[ee] = &*(diffN_volume_edge[ee].data().begin());
267  }
269  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
270  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_e,
271  diff_phi_v_e, nb_gauss_pts, base_polynomials);
272 
273  N_volume_face.resize(4, false);
274  diffN_volume_face.resize(4, false);
275  for (int ff = 0; ff != 4; ++ff) {
276  N_volume_face[ff].resize(
277  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
278  diffN_volume_face[ff].resize(
279  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
280  phi_v_f[ff] = &*(N_volume_face[ff].data().begin());
281  diff_phi_v_f[ff] = &*(diffN_volume_face[ff].data().begin());
282  }
284  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
285  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_f,
286  diff_phi_v_f, nb_gauss_pts, base_polynomials);
287 
288  N_volume_bubble.resize(
289  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
290  diffN_volume_bubble.resize(
291  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
292  phi_v = &*(N_volume_bubble.data().begin());
293  diff_phi_v = &*(diffN_volume_bubble.data().begin());
295  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
296  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v, diff_phi_v,
297  nb_gauss_pts, base_polynomials);
298 
299  // Set shape functions into data structure Shape functions hast to be put
300  // in arrays in order which guarantee hierarchical series of degrees of
301  // freedom, i.e. in other words dofs form sub-entities has to be group
302  // by order.
303 
306 
307  // faces
308  if (data.dataOnEntities[MBTRI].size() != 4) {
309  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
310  }
311  for (int ff = 0; ff != 4; ff++) {
312  data.dataOnEntities[MBTRI][ff].getN(base).resize(
313  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
314  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
315  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
316  if (NBFACETRI_AINSWORTH_HDIV(faces_order[ff]) == 0)
317  continue;
318  // face
319  double *base_ptr =
320  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
321  FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
322  &base_ptr[HVEC2], 3);
323  double *diff_base_ptr =
324  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
326  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
327  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
328  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
329  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
330  &diff_base_ptr[HVEC2_2], 9);
331  // face-face
332  boost::shared_ptr<FTensor::Tensor1<double *, 3> > t_base_f;
333  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> > t_diff_base_f;
334  if (NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]) > 0) {
335  base_ptr = phi_f[ff];
336  t_base_f = boost::shared_ptr<FTensor::Tensor1<double *, 3> >(
337  new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
338  &base_ptr[HVEC2], 3));
339  diff_base_ptr = diff_phi_f[ff];
340  t_diff_base_f = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> >(
342  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
343  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
344  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
345  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
346  &diff_base_ptr[HVEC2_2], 9));
347  }
348  // edge-face
349  base_ptr = phi_f_e[ff][0];
350  FTensor::Tensor1<double *, 3> t_base_f_e0(base_ptr, &base_ptr[HVEC1],
351  &base_ptr[HVEC2], 3);
352  diff_base_ptr = diff_phi_f_e[ff][0];
353  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e0(
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  base_ptr = phi_f_e[ff][1];
360  FTensor::Tensor1<double *, 3> t_base_f_e1(base_ptr, &base_ptr[HVEC1],
361  &base_ptr[HVEC2], 3);
362  diff_base_ptr = diff_phi_f_e[ff][1];
363  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e1(
364  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
365  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
366  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
367  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
368  &diff_base_ptr[HVEC2_2], 9);
369  base_ptr = phi_f_e[ff][2];
370  FTensor::Tensor1<double *, 3> t_base_f_e2(base_ptr, &base_ptr[HVEC1],
371  &base_ptr[HVEC2], 3);
372  diff_base_ptr = diff_phi_f_e[ff][2];
373  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e2(
374  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
375  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
376  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
377  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
378  &diff_base_ptr[HVEC2_2], 9);
379  for (int gg = 0; gg != nb_gauss_pts; gg++) {
380  for (int oo = 0; oo != faces_order[ff]; oo++) {
381  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
382  dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
383  t_base(i) = t_base_f_e0(i);
384  ++t_base;
385  ++t_base_f_e0;
386  t_diff_base(i, j) = t_diff_base_f_e0(i, j);
387  ++t_diff_base;
388  ++t_diff_base_f_e0;
389  t_base(i) = t_base_f_e1(i);
390  ++t_base;
391  ++t_base_f_e1;
392  t_diff_base(i, j) = t_diff_base_f_e1(i, j);
393  ++t_diff_base;
394  ++t_diff_base_f_e1;
395  t_base(i) = t_base_f_e2(i);
396  ++t_base;
397  ++t_base_f_e2;
398  t_diff_base(i, j) = t_diff_base_f_e2(i, j);
399  ++t_diff_base;
400  ++t_diff_base_f_e2;
401  }
402  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
403  dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
404  t_base(i) = (*t_base_f)(i);
405  ++t_base;
406  ++(*t_base_f);
407  t_diff_base(i, j) = (*t_diff_base_f)(i, j);
408  ++t_diff_base;
409  ++(*t_diff_base_f);
410  }
411  }
412  }
413  }
414 
415  // volume
416  data.dataOnEntities[MBTET][0].getN(base).resize(
417  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
418  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
419  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
420  if (NBVOLUMETET_AINSWORTH_HDIV(volume_order) > 0) {
421  double *base_ptr =
422  &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
423  FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
424  &base_ptr[HVEC2], 3);
425  double *diff_base_ptr =
426  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
428  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
429  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
430  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
431  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
432  &diff_base_ptr[HVEC2_2], 9);
433  // edges
434  std::vector<FTensor::Tensor1<double *, 3> > t_base_v_e;
435  t_base_v_e.reserve(6);
436  std::vector<FTensor::Tensor2<double *, 3, 3> > t_diff_base_v_e;
437  t_diff_base_v_e.reserve(6);
438  for (int ee = 0; ee != 6; ee++) {
439  base_ptr = phi_v_e[ee];
440  diff_base_ptr = diff_phi_v_e[ee];
441  t_base_v_e.push_back(FTensor::Tensor1<double *, 3>(
442  base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
443  t_diff_base_v_e.push_back(FTensor::Tensor2<double *, 3, 3>(
444  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
445  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
446  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
447  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
448  &diff_base_ptr[HVEC2_2], 9));
449  }
450  // faces
451  std::vector<FTensor::Tensor1<double *, 3> > t_base_v_f;
452  t_base_v_f.reserve(4);
453  std::vector<FTensor::Tensor2<double *, 3, 3> > t_diff_base_v_f;
454  t_diff_base_v_f.reserve(4);
455  if (NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order) > 0) {
456  for (int ff = 0; ff != 4; ff++) {
457  base_ptr = phi_v_f[ff];
458  diff_base_ptr = diff_phi_v_f[ff];
459  t_base_v_f.push_back(FTensor::Tensor1<double *, 3>(
460  base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
461  t_diff_base_v_f.push_back(FTensor::Tensor2<double *, 3, 3>(
462  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
463  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
464  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
465  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
466  &diff_base_ptr[HVEC2_2], 9));
467  }
468  }
469  boost::shared_ptr<FTensor::Tensor1<double *, 3> > t_base_v;
470  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> > t_diff_base_v;
471  if (NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order) > 0) {
472  base_ptr = phi_v;
473  t_base_v = boost::shared_ptr<FTensor::Tensor1<double *, 3> >(
474  new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
475  &base_ptr[HVEC2], 3));
476  diff_base_ptr = diff_phi_v;
477  t_diff_base_v = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> >(
479  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
480  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
481  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
482  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
483  &diff_base_ptr[HVEC2_2], 9));
484  }
485  for (int gg = 0; gg != nb_gauss_pts; gg++) {
486  for (int oo = 0; oo < volume_order; oo++) {
487  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
488  dd < NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
489  for (int ee = 0; ee < 6; ee++) {
490  t_base(i) = t_base_v_e[ee](i);
491  ++t_base;
492  ++t_base_v_e[ee];
493  t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
494  ++t_diff_base;
495  ++t_diff_base_v_e[ee];
496  }
497  }
498  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
499  dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
500  for (int ff = 0; ff < 4; ff++) {
501  t_base(i) = t_base_v_f[ff](i);
502  ++t_base;
503  ++t_base_v_f[ff];
504  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
505  ++t_diff_base;
506  ++t_diff_base_v_f[ff];
507  }
508  }
509  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
510  dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
511  t_base(i) = (*t_base_v)(i);
512  ++t_base;
513  ++(*t_base_v);
514  t_diff_base(i, j) = (*t_diff_base_v)(i, j);
515  ++t_diff_base;
516  ++(*t_diff_base_v);
517  }
518  }
519  }
520  }
521 
523 }
ublas::matrix< int > facesNodes
nodes on finite element faces
ublas::matrix< MatrixDouble > diffN_face_edge
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
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:477
#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:484
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:368
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
FieldApproximationBase
approximation base
Definition: definitions.h:143
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:596
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:280
#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:13
#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:407
#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:137
EntPolynomialBaseCtx * cTx
ublas::vector< MatrixDouble > N_face_bubble
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)

◆ getValueHdivDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBase ( MatrixDouble pts)
private

Definition at line 525 of file TetPolynomialBase.cpp.

525  {
527 
529  const FieldApproximationBase base = cTx->bAse;
530  if (base != DEMKOWICZ_JACOBI_BASE) {
531  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
532  "This should be used only with DEMKOWICZ_JACOBI_BASE "
533  "but base is %s",
534  ApproximationBaseNames[base]);
535  }
536  int nb_gauss_pts = pts.size2();
537 
538  int volume_order = data.dataOnEntities[MBTET][0].getDataOrder();
539 
540  int p_f[4];
541  double *phi_f[4];
542  double *diff_phi_f[4];
543 
544  // Calculate base function on tet faces
545  for (int ff = 0; ff != 4; ff++) {
546  int face_order = data.dataOnEntities[MBTRI][ff].getDataOrder();
547  int order = volume_order > face_order ? volume_order : face_order;
548  data.dataOnEntities[MBTRI][ff].getN(base).resize(
549  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
550  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
551  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
552  p_f[ff] = order;
553  phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
554  diff_phi_f[ff] =
555  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
557  continue;
559  &data.facesNodes(ff, 0), order,
560  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
561  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
562  phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
563  }
564 
565  // Calculate base functions in tet interior
566  if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
567  data.dataOnEntities[MBTET][0].getN(base).resize(
568  nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
569  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
570  nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
571  double *phi_v =
572  &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
573  double *diff_phi_v =
574  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
576  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
577  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
578  diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
579  }
580 
581  // Set size of face base correctly
582  for (int ff = 0; ff != 4; ff++) {
583  int face_order = data.dataOnEntities[MBTRI][ff].getDataOrder();
584  data.dataOnEntities[MBTRI][ff].getN(base).resize(
585  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
586  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
587  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
588  }
589 
591 }
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:617
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
static const char *const ApproximationBaseNames[]
Definition: definitions.h:158
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:762
#define NBFACETRI_DEMKOWICZ_HDIV(P)
FieldApproximationBase
approximation base
Definition: definitions.h:143
DataForcesAndSourcesCore & dAta
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
EntPolynomialBaseCtx * cTx
const int order
approximation order

◆ getValueL2()

MoFEMErrorCode TetPolynomialBase::getValueL2 ( MatrixDouble pts)
private

Definition at line 156 of file TetPolynomialBase.cpp.

156  {
158 
160  const FieldApproximationBase base = cTx->bAse;
161  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
162  double *diffL, const int dim) =
164 
165  int nb_gauss_pts = pts.size2();
166 
167  data.dataOnEntities[MBTET][0].getN(base).resize(
168  nb_gauss_pts,
169  NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getDataOrder()), false);
170  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
171  nb_gauss_pts,
172  3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getDataOrder()), false);
173 
175  data.dataOnEntities[MBTET][0].getDataOrder(),
176  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
177  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
178  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
179  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
180  nb_gauss_pts, base_polynomials);
181 
183 }
#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:477
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
FieldApproximationBase
approximation base
Definition: definitions.h:143
DataForcesAndSourcesCore & dAta
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
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 26 of file TetPolynomialBase.cpp.

27  {
28 
30  *iface = NULL;
31  if (uuid == IDD_TET_BASE_FUNCTION) {
32  *iface = const_cast<TetPolynomialBase *>(this);
34  } else {
35  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
36  }
37  ierr = BaseFunction::query_interface(uuid, iface);
38  CHKERRG(ierr);
40 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static const MOFEMuuid IDD_TET_BASE_FUNCTION
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: