v0.14.0
Public Member Functions | Static Public Member Functions | Protected Member Functions | Static Protected Member Functions | Protected 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 (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 
 TetPolynomialBase (const void *ptr=nullptr)
 
virtual ~TetPolynomialBase ()
 
MoFEMErrorCode getValue (MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
- Public Member Functions inherited from MoFEM::BaseFunction
virtual MoFEMErrorCode getValue (MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
- Public Member Functions inherited from MoFEM::BaseFunctionUnknownInterface
virtual ~BaseFunctionUnknownInterface ()=default
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface reference 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 ()=default
 

Static Public Member Functions

template<int SPACE>
static bool switchCacheBaseFace (FieldApproximationBase base, void *ptr)
 
template<int SPACE>
static bool switchCacheBaseInterior (FieldApproximationBase base, void *ptr)
 
template<int SPACE>
static bool switchCacheBrokenBaseInterior (FieldApproximationBase base, void *ptr)
 
template<int SPACE>
static void switchCacheBaseOn (FieldApproximationBase base, std::vector< void * > v)
 
template<int SPACE>
static void switchCacheBaseOff (FieldApproximationBase base, std::vector< void * > v)
 
template<int SPACE>
static void switchCacheBaseOn (std::vector< void * > v)
 
template<int SPACE>
static void switchCacheBaseOff (std::vector< void * > v)
 
static MoFEMErrorCode setDofsSideMap (const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &)
 Set map of dof to side number. More...
 
template<>
bool switchCacheBaseFace (FieldApproximationBase base, void *ptr)
 
template<>
bool switchCacheBaseInterior (FieldApproximationBase base, void *ptr)
 
template<>
bool switchCacheBrokenBaseInterior (FieldApproximationBase base, void *ptr)
 
template<>
void switchCacheBaseOn (FieldApproximationBase base, std::vector< void * > v)
 
template<>
void switchCacheBaseOff (FieldApproximationBase base, std::vector< void * > v)
 
template<>
void switchCacheBaseOn (std::vector< void * > v)
 
template<>
void switchCacheBaseOff (std::vector< void * > v)
 
template<>
bool switchCacheBaseInterior (FieldApproximationBase base, void *ptr)
 
template<>
void switchCacheBaseOn (FieldApproximationBase base, std::vector< void * > v)
 
template<>
void switchCacheBaseOff (FieldApproximationBase base, std::vector< void * > v)
 
template<>
void switchCacheBaseOn (std::vector< void * > v)
 
template<>
void switchCacheBaseOff (std::vector< void * > v)
 
template<>
bool switchCacheBaseFace (FieldApproximationBase base, void *ptr)
 
template<>
bool switchCacheBaseInterior (FieldApproximationBase base, void *ptr)
 
template<>
bool switchCacheBrokenBaseInterior (FieldApproximationBase base, void *ptr)
 
template<>
void switchCacheBaseOn (FieldApproximationBase base, std::vector< void * > v)
 
template<>
void switchCacheBaseOff (FieldApproximationBase base, std::vector< void * > v)
 
template<>
void switchCacheBaseOn (std::vector< void * > v)
 
template<>
void switchCacheBaseOff (std::vector< void * > v)
 
template<>
bool switchCacheBaseInterior (FieldApproximationBase base, void *ptr)
 
template<>
void switchCacheBaseOn (FieldApproximationBase base, std::vector< void * > v)
 
template<>
void switchCacheBaseOff (FieldApproximationBase base, std::vector< void * > v)
 
template<>
void switchCacheBaseOn (std::vector< void * > v)
 
template<>
void switchCacheBaseOff (std::vector< void * > v)
 
- Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode getLibVersion (Version &version)
 Get library version. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

Protected Member Functions

MoFEMErrorCode getValueH1 (MatrixDouble &pts)
 Get base functions for H1 space. More...
 
MoFEMErrorCode getValueL2 (MatrixDouble &pts)
 Get base functions for L2 space. More...
 
MoFEMErrorCode getValueHdiv (MatrixDouble &pts)
 Get base functions for Hdiv space. More...
 
MoFEMErrorCode getValueHcurl (MatrixDouble &pts)
 Get base functions for Hcurl space. More...
 
MoFEMErrorCode getValueH1AinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueH1BernsteinBezierBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueL2AinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueL2BernsteinBezierBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivAinsworthBaseImpl (MatrixDouble &pts, MatrixDouble &shape_functions, MatrixDouble &diff_shape_functions, int volume_order, std::array< int, 4 > &faces_order, std::array< int, 3 *4 > &faces_nodes, boost::function< int(int)> broken_nbfacetri_edge_hdiv, boost::function< int(int)> broken_nbfacetri_face_hdiv, boost::function< int(int)> broken_nbvolumetet_edge_hdiv, boost::function< int(int)> broken_nbvolumetet_face_hdiv, boost::function< int(int)> broken_nbvolumetet_volume_hdiv)
 
MoFEMErrorCode getValueHdivAinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivAinsworthBrokenBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurlAinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivDemkowiczBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivDemkowiczBrokenBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurlDemkowiczBase (MatrixDouble &pts)
 

Static Protected Member Functions

static MoFEMErrorCode setDofsSideMapHdiv (const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &dofs_side_map)
 Set the Dofs Side Map Hdiv object. More...
 

Protected Attributes

const void * vPtr
 
EntPolynomialBaseCtx * cTx
 
MatrixInt senseFaceAlpha
 
ublas::matrix< MatrixDouble > N_face_edge
 
ublas::vector< MatrixDouble > N_face_bubble
 
ublas::vector< MatrixDouble > N_volume_edge
 
ublas::vector< MatrixDouble > N_volume_face
 
MatrixDouble N_volume_bubble
 
ublas::matrix< MatrixDouble > diffN_face_edge
 
ublas::vector< MatrixDouble > diffN_face_bubble
 
ublas::vector< MatrixDouble > diffN_volume_edge
 
ublas::vector< MatrixDouble > diffN_volume_face
 
MatrixDouble diffN_volume_bubble
 

Additional Inherited Members

- Public Types inherited from MoFEM::BaseFunction
using DofsSideMap = multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > >>, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > >
 Map entity stype and side to element/entity dof index. More...
 

Detailed Description

Calculate base functions on tetrahedral.

Examples
EshelbianPlasticity.cpp.

Definition at line 17 of file TetPolynomialBase.hpp.

Constructor & Destructor Documentation

◆ TetPolynomialBase()

TetPolynomialBase::TetPolynomialBase ( const void *  ptr = nullptr)

Definition at line 95 of file TetPolynomialBase.cpp.

95 : vPtr(ptr) {}

◆ ~TetPolynomialBase()

TetPolynomialBase::~TetPolynomialBase ( )
virtual

Definition at line 97 of file TetPolynomialBase.cpp.

97  {
98  if (vPtr) {
99 
100  auto erase = [&](auto cache) {
101  if (cache.find(vPtr) != cache.end())
102  cache.erase(vPtr);
103  };
104 
105  for (auto b = 0; b != LASTBASE; ++b) {
106  erase(TetBaseCache::hdivBaseInterior[b]);
107  erase(TetBaseCache::hdivBrokenBaseInterior[b]);
108  erase(TetBaseCache::hDivBaseFace[b]);
109  }
110  }
111 }

Member Function Documentation

◆ getValue()

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

Reimplemented from MoFEM::BaseFunction.

Definition at line 2048 of file TetPolynomialBase.cpp.

2049  {
2050  MoFEMFunctionBegin;
2051 
2052  cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
2053 
2054  int nb_gauss_pts = pts.size2();
2055  if (!nb_gauss_pts)
2056  MoFEMFunctionReturnHot(0);
2057 
2058  if (pts.size1() < 3)
2059  SETERRQ(
2060  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2061  "Wrong dimension of pts, should be at least 3 rows with coordinates");
2062 
2063  if (cTx->bAse != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
2064  const FieldApproximationBase base = cTx->bAse;
2065  EntitiesFieldData &data = cTx->dAta;
2066  if (cTx->copyNodeBase == LASTBASE) {
2067  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
2068  false);
2069  CHKERR Tools::shapeFunMBTET(
2070  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
2071  &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
2072  } else {
2073  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
2074  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
2075  }
2076  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
2077  (unsigned int)nb_gauss_pts) {
2078  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2079  "Base functions or nodes has wrong number of integration points "
2080  "for base %s",
2081  ApproximationBaseNames[base]);
2082  }
2083  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
2084  std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
2085  data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
2086  }
2087 
2088  switch (cTx->spaceContinuity) {
2089  case CONTINUOUS:
2090 
2091  switch (cTx->sPace) {
2092  case H1:
2093  CHKERR getValueH1(pts);
2094  break;
2095  case HDIV:
2096  CHKERR getValueHdiv(pts);
2097  break;
2098  case HCURL:
2099  CHKERR getValueHcurl(pts);
2100  break;
2101  case L2:
2102  CHKERR getValueL2(pts);
2103  break;
2104  default:
2105  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2106  FieldSpaceNames[cTx->sPace]);
2107  }
2108  break;
2109 
2110  case DISCONTINUOUS:
2111 
2112  switch (cTx->sPace) {
2113  case HDIV:
2114  CHKERR getValueHdiv(pts);
2115  break;
2116  default:
2117  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2118  FieldSpaceNames[cTx->sPace]);
2119  }
2120  break;
2121 
2122  default:
2123  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
2124  }
2125 
2126  MoFEMFunctionReturn(0);
2127 }

◆ getValueH1()

MoFEMErrorCode TetPolynomialBase::getValueH1 ( MatrixDouble &  pts)
protected

Get base functions for H1 space.

Parameters
ptsmatrix of integration pts
Returns
MoFEMErrorCode
Note
matrix of integration points on rows has local coordinates of finite element on columns are integration pts.

Definition at line 113 of file TetPolynomialBase.cpp.

113  {
114  MoFEMFunctionBegin;
115 
116  switch (cTx->bAse) {
117  case AINSWORTH_LEGENDRE_BASE:
118  case AINSWORTH_LOBATTO_BASE:
119  CHKERR getValueH1AinsworthBase(pts);
120  break;
121  case AINSWORTH_BERNSTEIN_BEZIER_BASE:
122  CHKERR getValueH1BernsteinBezierBase(pts);
123  break;
124  default:
125  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
126  }
127 
128  MoFEMFunctionReturn(0);
129 }

◆ getValueH1AinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueH1AinsworthBase ( MatrixDouble &  pts)
protected

Definition at line 131 of file TetPolynomialBase.cpp.

131  {
132  MoFEMFunctionBegin;
133 
134  EntitiesFieldData &data = cTx->dAta;
135  const FieldApproximationBase base = cTx->bAse;
136  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
137  double *diffL, const int dim) =
138  cTx->basePolynomialsType0;
139 
140  int nb_gauss_pts = pts.size2();
141 
142  int sense[6], order[6];
143  if (data.spacesOnEntities[MBEDGE].test(H1)) {
144  // edges
145  if (data.dataOnEntities[MBEDGE].size() != 6) {
146  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
147  }
148  double *h1_edge_n[6], *diff_h1_egde_n[6];
149  for (int ee = 0; ee != 6; ++ee) {
150  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
151  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
152  "data inconsistency");
153  }
154  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
155  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
156  int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
157  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
158  false);
159  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
160  3 * nb_dofs, false);
161  h1_edge_n[ee] =
162  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
163  diff_h1_egde_n[ee] =
164  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
165  }
166  CHKERR H1_EdgeShapeFunctions_MBTET(
167  sense, order,
168  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
169  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
170  h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
171  } else {
172  for (int ee = 0; ee != 6; ++ee) {
173  data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
174  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
175  }
176  }
177 
178  if (data.spacesOnEntities[MBTRI].test(H1)) {
179  // faces
180  if (data.dataOnEntities[MBTRI].size() != 4) {
181  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
182  }
183  double *h1_face_n[4], *diff_h1_face_n[4];
184  for (int ff = 0; ff != 4; ++ff) {
185  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
186  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
187  "data inconsistency");
188  }
189  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getOrder());
190  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
191  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
192  false);
193  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
194  3 * nb_dofs, false);
195  h1_face_n[ff] =
196  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
197  diff_h1_face_n[ff] =
198  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
199  }
200  if (data.facesNodes.size1() != 4) {
201  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
202  }
203  if (data.facesNodes.size2() != 3) {
204  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
205  }
206  CHKERR H1_FaceShapeFunctions_MBTET(
207  &*data.facesNodes.data().begin(), order,
208  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
209  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
210  h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
211 
212  } else {
213  for (int ff = 0; ff != 4; ++ff) {
214  data.dataOnEntities[MBTRI][ff].getN(base).resize(0, false);
215  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, false);
216  }
217  }
218 
219  if (data.spacesOnEntities[MBTET].test(H1)) {
220  // volume
221  int order = data.dataOnEntities[MBTET][0].getOrder();
222  int nb_vol_dofs = NBVOLUMETET_H1(order);
223  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
224  false);
225  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
226  3 * nb_vol_dofs, false);
227  CHKERR H1_VolumeShapeFunctions_MBTET(
228  data.dataOnEntities[MBTET][0].getOrder(),
229  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
230  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
231  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
232  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
233  nb_gauss_pts, base_polynomials);
234  } else {
235  data.dataOnEntities[MBTET][0].getN(base).resize(0, 0, false);
236  data.dataOnEntities[MBTET][0].getDiffN(base).resize(0, 0, false);
237  }
238 
239  MoFEMFunctionReturn(0);
240 }

◆ getValueH1BernsteinBezierBase()

MoFEMErrorCode TetPolynomialBase::getValueH1BernsteinBezierBase ( MatrixDouble &  pts)
protected

Definition at line 243 of file TetPolynomialBase.cpp.

243  {
244  MoFEMFunctionBegin;
245 
246  EntitiesFieldData &data = cTx->dAta;
247  const std::string field_name = cTx->fieldName;
248  const int nb_gauss_pts = pts.size2();
249 
250  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
251  (unsigned int)nb_gauss_pts)
252  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
253  "Base functions or nodes has wrong number of integration points "
254  "for base %s",
255  ApproximationBaseNames[NOBASE]);
256  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
257 
258  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
259  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
260  if (!ptr)
261  ptr.reset(new MatrixInt());
262  return *ptr;
263  };
264 
265  auto get_base = [field_name](auto &data) -> MatrixDouble & {
266  auto &ptr = data.getBBNSharedPtr(field_name);
267  if (!ptr)
268  ptr.reset(new MatrixDouble());
269  return *ptr;
270  };
271 
272  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
273  auto &ptr = data.getBBDiffNSharedPtr(field_name);
274  if (!ptr)
275  ptr.reset(new MatrixDouble());
276  return *ptr;
277  };
278 
279  auto get_alpha_by_name_ptr =
280  [](auto &data,
281  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
282  return data.getBBAlphaIndicesSharedPtr(field_name);
283  };
284 
285  auto get_base_by_name_ptr =
286  [](auto &data,
287  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
288  return data.getBBNSharedPtr(field_name);
289  };
290 
291  auto get_diff_base_by_name_ptr =
292  [](auto &data,
293  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
294  return data.getBBDiffNSharedPtr(field_name);
295  };
296 
297  auto get_alpha_by_order_ptr =
298  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
299  return data.getBBAlphaIndicesByOrderSharedPtr(o);
300  };
301 
302  auto get_base_by_order_ptr =
303  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
304  return data.getBBNByOrderSharedPtr(o);
305  };
306 
307  auto get_diff_base_by_order_ptr =
308  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
309  return data.getBBDiffNByOrderSharedPtr(o);
310  };
311 
312  auto &vert_ent_data = data.dataOnEntities[MBVERTEX][0];
313  auto &vertex_alpha = get_alpha(vert_ent_data);
314  vertex_alpha.resize(4, 4, false);
315  vertex_alpha.clear();
316  for (int n = 0; n != 4; ++n)
317  vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
318 
319  auto &vert_get_n = get_base(vert_ent_data);
320  auto &vert_get_diff_n = get_diff_base(vert_ent_data);
321  vert_get_n.resize(nb_gauss_pts, 4, false);
322  vert_get_diff_n.resize(nb_gauss_pts, 12, false);
323  CHKERR BernsteinBezier::baseFunctionsTet(
324  1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
325  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &vert_get_n(0, 0),
326  &vert_get_diff_n(0, 0));
327  for (int n = 0; n != 4; ++n) {
328  const double f = boost::math::factorial<double>(
329  data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
330  for (int g = 0; g != nb_gauss_pts; ++g) {
331  vert_get_n(g, n) *= f;
332  for (int d = 0; d != 3; ++d)
333  vert_get_diff_n(g, 3 * n + d) *= f;
334  }
335  }
336 
337  // edges
338  if (data.spacesOnEntities[MBEDGE].test(H1)) {
339  if (data.dataOnEntities[MBEDGE].size() != 6)
340  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
341  "Wrong size of ent data");
342 
343  constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
344  {0, 3}, {1, 3}, {2, 3}};
345  for (int ee = 0; ee != 6; ++ee) {
346  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
347  const int sense = ent_data.getSense();
348  if (sense == 0)
349  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
350  "Sense of the edge unknown");
351  const int order = ent_data.getOrder();
352  const int nb_dofs = NBEDGE_H1(order);
353 
354  if (nb_dofs) {
355  if (get_alpha_by_order_ptr(ent_data, order)) {
356  get_alpha_by_name_ptr(ent_data, field_name) =
357  get_alpha_by_order_ptr(ent_data, order);
358  get_base_by_name_ptr(ent_data, field_name) =
359  get_base_by_order_ptr(ent_data, order);
360  get_diff_base_by_name_ptr(ent_data, field_name) =
361  get_diff_base_by_order_ptr(ent_data, order);
362  } else {
363  auto &get_n = get_base(ent_data);
364  auto &get_diff_n = get_diff_base(ent_data);
365  get_n.resize(nb_gauss_pts, nb_dofs, false);
366  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
367 
368  auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
369  edge_alpha.resize(nb_dofs, 4, false);
370  CHKERR BernsteinBezier::generateIndicesEdgeTet(ee, order,
371  &edge_alpha(0, 0));
372  if (sense == -1) {
373  for (int i = 0; i != edge_alpha.size1(); ++i) {
374  int a = edge_alpha(i, edges_nodes[ee][0]);
375  edge_alpha(i, edges_nodes[ee][0]) =
376  edge_alpha(i, edges_nodes[ee][1]);
377  edge_alpha(i, edges_nodes[ee][1]) = a;
378  }
379  }
380  CHKERR BernsteinBezier::baseFunctionsTet(
381  order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
382  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
383  &get_diff_n(0, 0));
384 
385  get_alpha_by_order_ptr(ent_data, order) =
386  get_alpha_by_name_ptr(ent_data, field_name);
387  get_base_by_order_ptr(ent_data, order) =
388  get_base_by_name_ptr(ent_data, field_name);
389  get_diff_base_by_order_ptr(ent_data, order) =
390  get_diff_base_by_name_ptr(ent_data, field_name);
391  }
392  }
393  }
394  } else {
395  for (int ee = 0; ee != 6; ++ee) {
396  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
397  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
398  auto &get_n = get_base(ent_data);
399  auto &get_diff_n = get_diff_base(ent_data);
400  get_n.resize(nb_gauss_pts, 0, false);
401  get_diff_n.resize(nb_gauss_pts, 0, false);
402  }
403  }
404 
405  // face
406  if (data.spacesOnEntities[MBTRI].test(H1)) {
407  if (data.dataOnEntities[MBTRI].size() != 4)
408  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
409  "Wrong size of ent data");
410  if (data.facesNodes.size1() != 4)
411  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
412  if (data.facesNodes.size2() != 3)
413  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
414 
415  for (int ff = 0; ff != 4; ++ff) {
416  auto &ent_data = data.dataOnEntities[MBTRI][ff];
417  const int order = ent_data.getOrder();
418  const int nb_dofs = NBFACETRI_H1(order);
419 
420  if (nb_dofs) {
421  if (get_alpha_by_order_ptr(ent_data, order)) {
422  get_alpha_by_name_ptr(ent_data, field_name) =
423  get_alpha_by_order_ptr(ent_data, order);
424  get_base_by_name_ptr(ent_data, field_name) =
425  get_base_by_order_ptr(ent_data, order);
426  get_diff_base_by_name_ptr(ent_data, field_name) =
427  get_diff_base_by_order_ptr(ent_data, order);
428  } else {
429 
430  auto &get_n = get_base(ent_data);
431  auto &get_diff_n = get_diff_base(ent_data);
432  get_n.resize(nb_gauss_pts, nb_dofs, false);
433  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
434 
435  auto &face_alpha = get_alpha(ent_data);
436  face_alpha.resize(nb_dofs, 4, false);
437 
438  CHKERR BernsteinBezier::generateIndicesTriTet(ff, order,
439  &face_alpha(0, 0));
440  senseFaceAlpha.resize(face_alpha.size1(), face_alpha.size2(), false);
441  senseFaceAlpha.clear();
442  constexpr int tri_nodes[4][3] = {
443  {0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
444  for (int d = 0; d != nb_dofs; ++d)
445  for (int n = 0; n != 3; ++n)
446  senseFaceAlpha(d, data.facesNodes(ff, n)) =
447  face_alpha(d, tri_nodes[ff][n]);
448  face_alpha.swap(senseFaceAlpha);
449  CHKERR BernsteinBezier::baseFunctionsTet(
450  order, lambda.size1(), face_alpha.size1(), &face_alpha(0, 0),
451  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
452  &get_diff_n(0, 0));
453 
454  get_alpha_by_order_ptr(ent_data, order) =
455  get_alpha_by_name_ptr(ent_data, field_name);
456  get_base_by_order_ptr(ent_data, order) =
457  get_base_by_name_ptr(ent_data, field_name);
458  get_diff_base_by_order_ptr(ent_data, order) =
459  get_diff_base_by_name_ptr(ent_data, field_name);
460  }
461  }
462  }
463  } else {
464  for (int ff = 0; ff != 4; ++ff) {
465  auto &ent_data = data.dataOnEntities[MBTRI][ff];
466  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
467  auto &get_n = get_base(ent_data);
468  auto &get_diff_n = get_diff_base(ent_data);
469  get_n.resize(nb_gauss_pts, 0, false);
470  get_diff_n.resize(nb_gauss_pts, 0, false);
471  }
472  }
473 
474  if (data.spacesOnEntities[MBTET].test(H1)) {
475  if (data.dataOnEntities[MBTET].size() != 1)
476  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
477  "Wrong size ent of ent data");
478 
479  auto &ent_data = data.dataOnEntities[MBTET][0];
480  const int order = ent_data.getOrder();
481  const int nb_dofs = NBVOLUMETET_H1(order);
482  if (get_alpha_by_order_ptr(ent_data, order)) {
483  get_alpha_by_name_ptr(ent_data, field_name) =
484  get_alpha_by_order_ptr(ent_data, order);
485  get_base_by_name_ptr(ent_data, field_name) =
486  get_base_by_order_ptr(ent_data, order);
487  get_diff_base_by_name_ptr(ent_data, field_name) =
488  get_diff_base_by_order_ptr(ent_data, order);
489  } else {
490 
491  auto &get_n = get_base(ent_data);
492  auto &get_diff_n = get_diff_base(ent_data);
493  get_n.resize(nb_gauss_pts, nb_dofs, false);
494  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
495  if (nb_dofs) {
496  auto &tet_alpha = get_alpha(ent_data);
497  tet_alpha.resize(nb_dofs, 4, false);
498 
499  CHKERR BernsteinBezier::generateIndicesTetTet(order, &tet_alpha(0, 0));
500  CHKERR BernsteinBezier::baseFunctionsTet(
501  order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
502  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
503  &get_diff_n(0, 0));
504 
505  get_alpha_by_order_ptr(ent_data, order) =
506  get_alpha_by_name_ptr(ent_data, field_name);
507  get_base_by_order_ptr(ent_data, order) =
508  get_base_by_name_ptr(ent_data, field_name);
509  get_diff_base_by_order_ptr(ent_data, order) =
510  get_diff_base_by_name_ptr(ent_data, field_name);
511  }
512  }
513  } else {
514  auto &ent_data = data.dataOnEntities[MBTET][0];
515  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
516  auto &get_n = get_base(ent_data);
517  auto &get_diff_n = get_diff_base(ent_data);
518  get_n.resize(nb_gauss_pts, 0, false);
519  get_diff_n.resize(nb_gauss_pts, 0, false);
520  }
521 
522  MoFEMFunctionReturn(0);
523 }

◆ getValueHcurl()

MoFEMErrorCode TetPolynomialBase::getValueHcurl ( MatrixDouble &  pts)
protected

Get base functions for Hcurl space.

Parameters
ptsmatrix of integration pts
Returns
MoFEMErrorCode
Note
matrix of integration points on rows has local coordinates of finite element on columns are integration pts.

Definition at line 2029 of file TetPolynomialBase.cpp.

2029  {
2030  MoFEMFunctionBegin;
2031 
2032  switch (cTx->bAse) {
2033  case AINSWORTH_LEGENDRE_BASE:
2034  case AINSWORTH_LOBATTO_BASE:
2035  CHKERR getValueHcurlAinsworthBase(pts);
2036  break;
2037  case DEMKOWICZ_JACOBI_BASE:
2038  CHKERR getValueHcurlDemkowiczBase(pts);
2039  break;
2040  default:
2041  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
2042  }
2043 
2044  MoFEMFunctionReturn(0);
2045 }

◆ getValueHcurlAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlAinsworthBase ( MatrixDouble &  pts)
protected

Definition at line 1779 of file TetPolynomialBase.cpp.

1779  {
1780  MoFEMFunctionBegin;
1781 
1782  EntitiesFieldData &data = cTx->dAta;
1783  const FieldApproximationBase base = cTx->bAse;
1784  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
1785  double *diffL, const int dim) =
1786  cTx->basePolynomialsType0;
1787 
1788  int nb_gauss_pts = pts.size2();
1789 
1790  // edges
1791  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1792  int sense[6], order[6];
1793  if (data.dataOnEntities[MBEDGE].size() != 6) {
1794  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1795  }
1796  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1797  for (int ee = 0; ee != 6; ee++) {
1798  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1799  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1800  "data inconsistency");
1801  }
1802  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1803  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1804  int nb_dofs =
1805  NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1806  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1807  3 * nb_dofs, false);
1808  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1809  9 * nb_dofs, false);
1810  hcurl_edge_n[ee] =
1811  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1812  diff_hcurl_edge_n[ee] =
1813  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1814  }
1815  CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(
1816  sense, order,
1817  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1818  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1819  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
1820  } else {
1821  for (int ee = 0; ee != 6; ee++) {
1822  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1823  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1824  false);
1825  }
1826  }
1827 
1828  // triangles
1829  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1830  int order[4];
1831  // faces
1832  if (data.dataOnEntities[MBTRI].size() != 4) {
1833  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1834  }
1835  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1836  for (int ff = 0; ff != 4; ff++) {
1837  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1838  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1839  "data inconsistency");
1840  }
1841  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1842  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
1843  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1844  3 * nb_dofs, false);
1845  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1846  9 * nb_dofs, false);
1847  hcurl_base_n[ff] =
1848  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1849  diff_hcurl_base_n[ff] =
1850  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1851  }
1852  if (data.facesNodes.size1() != 4) {
1853  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1854  }
1855  if (data.facesNodes.size2() != 3) {
1856  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1857  }
1858  CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET(
1859  &*data.facesNodes.data().begin(), order,
1860  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1861  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1862  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
1863  } else {
1864  for (int ff = 0; ff != 4; ff++) {
1865  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1866  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1867  false);
1868  }
1869  }
1870 
1871  if (data.spacesOnEntities[MBTET].test(HCURL)) {
1872 
1873  // volume
1874  int order = data.dataOnEntities[MBTET][0].getOrder();
1875  int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
1876  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1877  3 * nb_vol_dofs, false);
1878  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1879  9 * nb_vol_dofs, false);
1880  CHKERR Hcurl_Ainsworth_VolumeFunctions_MBTET(
1881  data.dataOnEntities[MBTET][0].getOrder(),
1882  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1883  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1884  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1885  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1886  nb_gauss_pts, base_polynomials);
1887 
1888  } else {
1889  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1890  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1891  }
1892 
1893  MoFEMFunctionReturn(0);
1894 }

◆ getValueHcurlDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlDemkowiczBase ( MatrixDouble &  pts)
protected

Definition at line 1897 of file TetPolynomialBase.cpp.

1897  {
1898  MoFEMFunctionBegin;
1899 
1900  EntitiesFieldData &data = cTx->dAta;
1901  const FieldApproximationBase base = cTx->bAse;
1902  if (base != DEMKOWICZ_JACOBI_BASE) {
1903  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1904  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1905  "but base is %s",
1906  ApproximationBaseNames[base]);
1907  }
1908 
1909  int nb_gauss_pts = pts.size2();
1910 
1911  // edges
1912  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1913  int sense[6], order[6];
1914  if (data.dataOnEntities[MBEDGE].size() != 6) {
1915  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1916  "wrong size of data structure, expected space for six edges "
1917  "but is %d",
1918  data.dataOnEntities[MBEDGE].size());
1919  }
1920  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1921  for (int ee = 0; ee != 6; ee++) {
1922  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1923  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1924  "orintation of edges is not set");
1925  }
1926  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1927  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1928  int nb_dofs =
1929  NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1930  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1931  3 * nb_dofs, false);
1932  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1933  9 * nb_dofs, false);
1934  hcurl_edge_n[ee] =
1935  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1936  diff_hcurl_edge_n[ee] =
1937  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1938  }
1939  CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(
1940  sense, order,
1941  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1942  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1943  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
1944  } else {
1945  // No DOFs on edges, resize base function matrices, indicating that no
1946  // dofs on them.
1947  for (int ee = 0; ee != 6; ee++) {
1948  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1949  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1950  false);
1951  }
1952  }
1953 
1954  // triangles
1955  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1956  int order[4];
1957  // faces
1958  if (data.dataOnEntities[MBTRI].size() != 4) {
1959  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1960  "data structure for storing face h-curl base have wrong size "
1961  "should be four but is %d",
1962  data.dataOnEntities[MBTRI].size());
1963  }
1964  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1965  for (int ff = 0; ff != 4; ff++) {
1966  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1967  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1968  "orintation of face is not set");
1969  }
1970  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1971  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order[ff]);
1972  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1973  3 * nb_dofs, false);
1974  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1975  9 * nb_dofs, false);
1976  hcurl_base_n[ff] =
1977  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1978  diff_hcurl_base_n[ff] =
1979  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1980  }
1981  if (data.facesNodes.size1() != 4) {
1982  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1983  "data inconsistency, should be four faces");
1984  }
1985  if (data.facesNodes.size2() != 3) {
1986  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1987  "data inconsistency, should be three nodes on face");
1988  }
1989  CHKERR Hcurl_Demkowicz_FaceBaseFunctions_MBTET(
1990  &*data.facesNodes.data().begin(), order,
1991  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1992  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1993  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
1994  } else {
1995  // No DOFs on faces, resize base function matrices, indicating that no
1996  // dofs on them.
1997  for (int ff = 0; ff != 4; ff++) {
1998  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1999  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
2000  false);
2001  }
2002  }
2003 
2004  if (data.spacesOnEntities[MBTET].test(HCURL)) {
2005  // volume
2006  int order = data.dataOnEntities[MBTET][0].getOrder();
2007  int nb_vol_dofs = NBVOLUMETET_DEMKOWICZ_HCURL(order);
2008  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
2009  3 * nb_vol_dofs, false);
2010  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
2011  9 * nb_vol_dofs, false);
2012  CHKERR Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(
2013  data.dataOnEntities[MBTET][0].getOrder(),
2014  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
2015  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
2016  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
2017  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
2018  nb_gauss_pts);
2019  } else {
2020  // No DOFs on faces, resize base function matrices, indicating that no
2021  // dofs on them.
2022  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
2023  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
2024  }
2025 
2026  MoFEMFunctionReturn(0);
2027 }

◆ getValueHdiv()

MoFEMErrorCode TetPolynomialBase::getValueHdiv ( MatrixDouble &  pts)
protected

Get base functions for Hdiv space.

Parameters
ptsmatrix of integration pts
Returns
MoFEMErrorCode
Note
matrix of integration points on rows has local coordinates of finite element on columns are integration pts.

Definition at line 1744 of file TetPolynomialBase.cpp.

1744  {
1745  MoFEMFunctionBeginHot;
1746 
1747  switch (cTx->spaceContinuity) {
1748  case CONTINUOUS:
1749  switch (cTx->bAse) {
1750  case AINSWORTH_LEGENDRE_BASE:
1751  case AINSWORTH_LOBATTO_BASE:
1752  return getValueHdivAinsworthBase(pts);
1753  case DEMKOWICZ_JACOBI_BASE:
1754  return getValueHdivDemkowiczBase(pts);
1755  default:
1756  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1757  }
1758  break;
1759  case DISCONTINUOUS:
1760  switch (cTx->bAse) {
1761  case AINSWORTH_LEGENDRE_BASE:
1762  case AINSWORTH_LOBATTO_BASE:
1763  return getValueHdivAinsworthBrokenBase(pts);
1764  case DEMKOWICZ_JACOBI_BASE:
1765  return getValueHdivDemkowiczBrokenBase(pts);
1766  default:
1767  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1768  }
1769  break;
1770 
1771  default:
1772  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1773  }
1774 
1775  MoFEMFunctionReturnHot(0);
1776 }

◆ getValueHdivAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBase ( MatrixDouble &  pts)
protected

Definition at line 909 of file TetPolynomialBase.cpp.

909  {
910  MoFEMFunctionBegin;
911 
912  std::array<int, 4> faces_order;
913  std::array<int, 4 * 3> faces_nodes;
914 
915  FieldApproximationBase base = cTx->bAse;
916  EntitiesFieldData &data = cTx->dAta;
917 
918  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
919  std::copy(data.facesNodes.data().begin(), data.facesNodes.data().end(),
920  faces_nodes.begin());
921  for (int ff = 0; ff != 4; ff++) {
922  faces_order[ff] = cTx->dAta.dataOnEntities[MBTRI][ff].getOrder();
923  }
924 
925  CHKERR getValueHdivAinsworthBaseImpl(
926  pts, data.dataOnEntities[MBVERTEX][0].getN(base),
927  data.dataOnEntities[MBVERTEX][0].getDiffN(base), volume_order,
928  faces_order, faces_nodes,
929 
930  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv,
931  AinsworthOrderHooks::broken_nbfacetri_face_hdiv,
932  AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv,
933  AinsworthOrderHooks::broken_nbvolumetet_face_hdiv,
934  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv
935 
936  );
937 
938  // Set shape functions into data structure Shape functions hast to be put
939  // in arrays in order which guarantee hierarchical series of degrees of
940  // freedom, i.e. in other words dofs form sub-entities has to be group
941  // by order.
942 
943  FTENSOR_INDEX(3, i);
944  FTENSOR_INDEX(3, j);
945 
946  int nb_gauss_pts = pts.size2();
947 
948  // faces
949  if (data.dataOnEntities[MBTRI].size() != 4) {
950  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
951  }
952 
953  // face-face
954  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_f_f[] = {
955  getFTensor1FromPtr<3>(&*(N_face_bubble[0].data().begin())),
956  getFTensor1FromPtr<3>(&*(N_face_bubble[1].data().begin())),
957  getFTensor1FromPtr<3>(&*(N_face_bubble[2].data().begin())),
958  getFTensor1FromPtr<3>(&*(N_face_bubble[3].data().begin()))};
959  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_f[] = {
960  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[0].data().begin())),
961  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[1].data().begin())),
962  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[2].data().begin())),
963  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[3].data().begin()))};
964  // face-edge
965  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_f_e[] = {
966  getFTensor1FromPtr<3>(&*(N_face_edge(0, 0).data().begin())),
967  getFTensor1FromPtr<3>(&*(N_face_edge(0, 1).data().begin())),
968  getFTensor1FromPtr<3>(&*(N_face_edge(0, 2).data().begin())),
969  getFTensor1FromPtr<3>(&*(N_face_edge(1, 0).data().begin())),
970  getFTensor1FromPtr<3>(&*(N_face_edge(1, 1).data().begin())),
971  getFTensor1FromPtr<3>(&*(N_face_edge(1, 2).data().begin())),
972  getFTensor1FromPtr<3>(&*(N_face_edge(2, 0).data().begin())),
973  getFTensor1FromPtr<3>(&*(N_face_edge(2, 1).data().begin())),
974  getFTensor1FromPtr<3>(&*(N_face_edge(2, 2).data().begin())),
975  getFTensor1FromPtr<3>(&*(N_face_edge(3, 0).data().begin())),
976  getFTensor1FromPtr<3>(&*(N_face_edge(3, 1).data().begin())),
977  getFTensor1FromPtr<3>(&*(N_face_edge(3, 2).data().begin()))};
978  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_e[] = {
979  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 0).data().begin())),
980  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 1).data().begin())),
981  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 2).data().begin())),
982  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 0).data().begin())),
983  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 1).data().begin())),
984  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 2).data().begin())),
985  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 0).data().begin())),
986  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 1).data().begin())),
987  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 2).data().begin())),
988  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 0).data().begin())),
989  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 1).data().begin())),
990  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 2).data().begin()))};
991 
992  for (int ff = 0; ff != 4; ff++) {
993  int face_order = faces_order[ff];
994  auto face_dofs =
995  3 * NBFACETRI_AINSWORTH_EDGE_HDIV(
996  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order)) +
997  NBFACETRI_AINSWORTH_FACE_HDIV(
998  AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order));
999  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1000  3 * face_dofs, false);
1001  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1002  9 * face_dofs, false);
1003  if (face_dofs) {
1004  double *base_ptr =
1005  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1006  double *diff_base_ptr =
1007  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1008  auto t_base = getFTensor1FromPtr<3>(base_ptr);
1009  auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1010 
1011  auto max_face_order =
1012  std::max(face_order,
1013  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order));
1014  max_face_order =
1015  std::max(max_face_order,
1016  AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order));
1017 
1018  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1019  for (int oo = 0; oo != max_face_order; oo++) {
1020 
1021  // face-edge
1022  if (oo < AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order))
1023  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
1024  dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1025  for (int ee = 0; ee != 3; ++ee) {
1026  t_base(i) = t_base_f_e[ff * 3 + ee](i);
1027  ++t_base;
1028  ++t_base_f_e[ff * 3 + ee];
1029  }
1030  for (int ee = 0; ee != 3; ++ee) {
1031  t_diff_base(i, j) = t_diff_base_f_e[ff * 3 + ee](i, j);
1032  ++t_diff_base;
1033  ++t_diff_base_f_e[ff * 3 + ee];
1034  }
1035  }
1036 
1037  // face-face
1038  if (oo < AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order))
1039  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
1040  dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1041  t_base(i) = t_base_f_f[ff](i);
1042  ++t_base;
1043  ++t_base_f_f[ff];
1044  t_diff_base(i, j) = t_diff_base_f_f[ff](i, j);
1045  ++t_diff_base;
1046  ++t_diff_base_f_f[ff];
1047  }
1048  }
1049  }
1050  }
1051  }
1052 
1053  // volume
1054  int volume_dofs =
1055  6 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(
1056  AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order)) +
1057  4 * NBVOLUMETET_AINSWORTH_FACE_HDIV(
1058  AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order)) +
1059  NBVOLUMETET_AINSWORTH_VOLUME_HDIV(
1060  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order));
1061  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * volume_dofs,
1062  false);
1063  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1064  9 * volume_dofs, false);
1065  if (volume_dofs) {
1066  double *base_ptr =
1067  &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1068  double *diff_base_ptr =
1069  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1070  auto t_base = getFTensor1FromPtr<3>(base_ptr);
1071  auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1072 
1073  // volume-edge
1074  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_e[] = {
1075  getFTensor1FromPtr<3>(&*N_volume_edge[0].data().begin()),
1076  getFTensor1FromPtr<3>(&*N_volume_edge[1].data().begin()),
1077  getFTensor1FromPtr<3>(&*N_volume_edge[2].data().begin()),
1078  getFTensor1FromPtr<3>(&*N_volume_edge[3].data().begin()),
1079  getFTensor1FromPtr<3>(&*N_volume_edge[4].data().begin()),
1080  getFTensor1FromPtr<3>(&*N_volume_edge[5].data().begin())};
1081  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_e[] = {
1082  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[0].data().begin()),
1083  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[1].data().begin()),
1084  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[2].data().begin()),
1085  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[3].data().begin()),
1086  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[4].data().begin()),
1087  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[5].data().begin())};
1088 
1089  // volume-faces
1090  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_f[] = {
1091  getFTensor1FromPtr<3>(&*N_volume_face[0].data().begin()),
1092  getFTensor1FromPtr<3>(&*N_volume_face[1].data().begin()),
1093  getFTensor1FromPtr<3>(&*N_volume_face[2].data().begin()),
1094  getFTensor1FromPtr<3>(&*N_volume_face[3].data().begin())};
1095  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1096  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[0].data().begin()),
1097  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[1].data().begin()),
1098  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[2].data().begin()),
1099  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[3].data().begin())};
1100 
1101  // volume-bubble
1102  base_ptr = &*(N_volume_bubble.data().begin());
1103  diff_base_ptr = &*(diffN_volume_bubble.data().begin());
1104  auto t_base_v = getFTensor1FromPtr<3>(base_ptr);
1105  auto t_diff_base_v = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1106 
1107  auto max_volume_order = std::max(
1108  volume_order,
1109  AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order));
1110  max_volume_order = std::max(
1111  max_volume_order,
1112  AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order));
1113  max_volume_order = std::max(
1114  max_volume_order,
1115  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order));
1116 
1117  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1118  for (int oo = 0; oo < max_volume_order; oo++) {
1119 
1120  // volume-edge
1121  if (oo <
1122  AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order))
1123  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
1124  dd != NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1125  for (int ee = 0; ee < 6; ee++) {
1126  t_base(i) = t_base_v_e[ee](i);
1127  ++t_base;
1128  ++t_base_v_e[ee];
1129  t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
1130  ++t_diff_base;
1131  ++t_diff_base_v_e[ee];
1132  }
1133  }
1134 
1135  // volume-face
1136  if (oo <
1137  AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order))
1138  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
1139  dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1140  for (int ff = 0; ff < 4; ff++) {
1141  t_base(i) = t_base_v_f[ff](i);
1142  ++t_base;
1143  ++t_base_v_f[ff];
1144  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1145  ++t_diff_base;
1146  ++t_diff_base_v_f[ff];
1147  }
1148  }
1149 
1150  // volume-bubble
1151  if (oo <
1152  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order))
1153  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
1154  dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
1155  t_base(i) = t_base_v(i);
1156  ++t_base;
1157  ++t_base_v;
1158  t_diff_base(i, j) = t_diff_base_v(i, j);
1159  ++t_diff_base;
1160  ++t_diff_base_v;
1161  }
1162  }
1163  }
1164  }
1165 
1166  MoFEMFunctionReturn(0);
1167 }

◆ getValueHdivAinsworthBaseImpl()

MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBaseImpl ( MatrixDouble &  pts,
MatrixDouble &  shape_functions,
MatrixDouble &  diff_shape_functions,
int  volume_order,
std::array< int, 4 > &  faces_order,
std::array< int, 3 *4 > &  faces_nodes,
boost::function< int(int)>  broken_nbfacetri_edge_hdiv,
boost::function< int(int)>  broken_nbfacetri_face_hdiv,
boost::function< int(int)>  broken_nbvolumetet_edge_hdiv,
boost::function< int(int)>  broken_nbvolumetet_face_hdiv,
boost::function< int(int)>  broken_nbvolumetet_volume_hdiv 
)
protected

Definition at line 781 of file TetPolynomialBase.cpp.

796  {
797  MoFEMFunctionBegin;
798 
799  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
800  double *diffL, const int dim) =
801  cTx->basePolynomialsType0;
802 
803  int nb_gauss_pts = pts.size2();
804 
805  // face shape functions
806 
807  double *phi_f_e[4][3];
808  double *phi_f[4];
809  double *diff_phi_f_e[4][3];
810  double *diff_phi_f[4];
811 
812  N_face_edge.resize(4, 3, false);
813  N_face_bubble.resize(4, false);
814  diffN_face_edge.resize(4, 3, false);
815  diffN_face_bubble.resize(4, false);
816 
817  for (int ff = 0; ff != 4; ++ff) {
818  const auto face_edge_dofs = NBFACETRI_AINSWORTH_EDGE_HDIV(
819  broken_nbfacetri_edge_hdiv(faces_order[ff]));
820  // three edges on face
821  for (int ee = 0; ee < 3; ee++) {
822  N_face_edge(ff, ee).resize(nb_gauss_pts, 3 * face_edge_dofs, false);
823  diffN_face_edge(ff, ee).resize(nb_gauss_pts, 9 * face_edge_dofs, false);
824  phi_f_e[ff][ee] = &*N_face_edge(ff, ee).data().begin();
825  diff_phi_f_e[ff][ee] = &*diffN_face_edge(ff, ee).data().begin();
826  }
827  auto face_bubble_dofs = NBFACETRI_AINSWORTH_FACE_HDIV(
828  broken_nbfacetri_face_hdiv(faces_order[ff]));
829  N_face_bubble[ff].resize(nb_gauss_pts, 3 * face_bubble_dofs, false);
830  diffN_face_bubble[ff].resize(nb_gauss_pts, 9 * face_bubble_dofs, false);
831  phi_f[ff] = &*(N_face_bubble[ff].data().begin());
832  diff_phi_f[ff] = &*(diffN_face_bubble[ff].data().begin());
833  }
834 
835  constexpr int nb_nodes_on_tet = 4;
836 
837  for (int ff = 0; ff < 4; ff++) {
838  CHKERR Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
839  &faces_nodes[3 * ff], broken_nbfacetri_edge_hdiv(faces_order[ff]),
840  &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
841  phi_f_e[ff], diff_phi_f_e[ff], nb_gauss_pts, nb_nodes_on_tet,
842  base_polynomials);
843  }
844 
845  for (int ff = 0; ff < 4; ff++) {
846  CHKERR Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
847  &faces_nodes[3 * ff], broken_nbfacetri_face_hdiv(faces_order[ff]),
848  &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
849  phi_f[ff], diff_phi_f[ff], nb_gauss_pts, nb_nodes_on_tet,
850  base_polynomials);
851  }
852 
853  // volume shape functions
854 
855  double *phi_v_e[6];
856  double *phi_v_f[4];
857  double *phi_v;
858  double *diff_phi_v_e[6];
859  double *diff_phi_v_f[4];
860  double *diff_phi_v;
861 
862  const auto volume_edge_dofs = NBVOLUMETET_AINSWORTH_EDGE_HDIV(
863  broken_nbvolumetet_edge_hdiv(volume_order));
864  N_volume_edge.resize(6, false);
865  diffN_volume_edge.resize(6, false);
866  for (int ee = 0; ee != 6; ++ee) {
867  N_volume_edge[ee].resize(nb_gauss_pts, 3 * volume_edge_dofs, false);
868  diffN_volume_edge[ee].resize(nb_gauss_pts, 9 * volume_edge_dofs, false);
869  phi_v_e[ee] = &*(N_volume_edge[ee].data().begin());
870  diff_phi_v_e[ee] = &*(diffN_volume_edge[ee].data().begin());
871  }
872  if (volume_edge_dofs)
873  CHKERR Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(
874  broken_nbvolumetet_edge_hdiv(volume_order),
875  &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
876  phi_v_e, diff_phi_v_e, nb_gauss_pts, base_polynomials);
877 
878  const auto volume_face_dofs = NBVOLUMETET_AINSWORTH_FACE_HDIV(
879  broken_nbvolumetet_face_hdiv(volume_order));
880  N_volume_face.resize(4, false);
881  diffN_volume_face.resize(4, false);
882  for (int ff = 0; ff != 4; ++ff) {
883  N_volume_face[ff].resize(nb_gauss_pts, 3 * volume_face_dofs, false);
884  diffN_volume_face[ff].resize(nb_gauss_pts, 9 * volume_face_dofs, false);
885  phi_v_f[ff] = &*(N_volume_face[ff].data().begin());
886  diff_phi_v_f[ff] = &*(diffN_volume_face[ff].data().begin());
887  }
888  if (volume_face_dofs)
889  CHKERR Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(
890  broken_nbvolumetet_face_hdiv(volume_order),
891  &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
892  phi_v_f, diff_phi_v_f, nb_gauss_pts, base_polynomials);
893 
894  auto volume_bubble_dofs = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(
895  broken_nbvolumetet_volume_hdiv(volume_order));
896  N_volume_bubble.resize(nb_gauss_pts, 3 * volume_bubble_dofs, false);
897  diffN_volume_bubble.resize(nb_gauss_pts, 9 * volume_bubble_dofs, false);
898  phi_v = &*(N_volume_bubble.data().begin());
899  diff_phi_v = &*(diffN_volume_bubble.data().begin());
900  if (volume_bubble_dofs)
901  CHKERR Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(
902  broken_nbvolumetet_volume_hdiv(volume_order),
903  &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
904  phi_v, diff_phi_v, nb_gauss_pts, base_polynomials);
905 
906  MoFEMFunctionReturn(0);
907 }

◆ getValueHdivAinsworthBrokenBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBrokenBase ( MatrixDouble &  pts)
protected

Definition at line 1170 of file TetPolynomialBase.cpp.

1170  {
1171  MoFEMFunctionBegin;
1172 
1173  // Set shape functions into data structure Shape functions has to be put
1174  // in arrays in order which guarantee hierarchical series of degrees of
1175  // freedom, i.e. in other words dofs form sub-entities has to be group
1176  // by order.
1177 
1178  FieldApproximationBase base = cTx->bAse;
1179  EntitiesFieldData &data = cTx->dAta;
1180 
1181  FTENSOR_INDEX(3, i);
1182  FTENSOR_INDEX(3, j);
1183 
1184  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1185  int nb_gauss_pts = pts.size2();
1186  int nb_dofs_face =
1187  3 * NBFACETRI_AINSWORTH_EDGE_HDIV(
1188  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(volume_order)) +
1189  NBFACETRI_AINSWORTH_FACE_HDIV(
1190  AinsworthOrderHooks::broken_nbfacetri_face_hdiv(volume_order));
1191  int nb_dofs_volume =
1192  6 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(
1193  AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order)) +
1194  4 * NBVOLUMETET_AINSWORTH_FACE_HDIV(
1195  AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order)) +
1196  NBVOLUMETET_AINSWORTH_VOLUME_HDIV(
1197  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order));
1198 
1199  int nb_dofs = 4 * nb_dofs_face + nb_dofs_volume;
1200  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
1201  false);
1202  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 9 * nb_dofs,
1203  false);
1204  if (nb_dofs == 0)
1205  MoFEMFunctionReturnHot(0);
1206 
1207  auto get_interior_cache = [this](auto base) -> TetBaseCache::BaseCacheMI * {
1208  if (vPtr) {
1209  auto it = TetBaseCache::hdivBrokenBaseInterior[base].find(vPtr);
1210  if (it != TetBaseCache::hdivBrokenBaseInterior[base].end()) {
1211  return &it->second;
1212  }
1213  }
1214  return nullptr;
1215  };
1216 
1217  auto interior_cache_ptr = get_interior_cache(base);
1218 
1219  if (interior_cache_ptr) {
1220  auto it =
1221  interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
1222  if (it != interior_cache_ptr->end()) {
1223  noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1224  noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1225  MoFEMFunctionReturnHot(0);
1226  }
1227  }
1228 
1229  std::array<int, 4 * 3> faces_nodes = {0, 1, 3, 1, 2, 3, 0, 3, 2, 0, 2, 1};
1230  std::array<int, 4> faces_order{volume_order, volume_order, volume_order,
1231  volume_order};
1232  CHKERR getValueHdivAinsworthBaseImpl(
1233  pts, data.dataOnEntities[MBVERTEX][0].getN(base),
1234  data.dataOnEntities[MBVERTEX][0].getDiffN(base), volume_order,
1235  faces_order, faces_nodes,
1236 
1237  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv,
1238  AinsworthOrderHooks::broken_nbfacetri_face_hdiv,
1239  AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv,
1240  AinsworthOrderHooks::broken_nbvolumetet_face_hdiv,
1241  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv
1242 
1243  );
1244 
1245  auto *base_ptr = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1246  auto *diff_base_ptr =
1247  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1248  auto t_base = getFTensor1FromPtr<3>(base_ptr);
1249  auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1250 
1251  // face-edge
1252  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_f_e[] = {
1253  getFTensor1FromPtr<3>(&*(N_face_edge(0, 0).data().begin())),
1254  getFTensor1FromPtr<3>(&*(N_face_edge(0, 1).data().begin())),
1255  getFTensor1FromPtr<3>(&*(N_face_edge(0, 2).data().begin())),
1256  getFTensor1FromPtr<3>(&*(N_face_edge(1, 0).data().begin())),
1257  getFTensor1FromPtr<3>(&*(N_face_edge(1, 1).data().begin())),
1258  getFTensor1FromPtr<3>(&*(N_face_edge(1, 2).data().begin())),
1259  getFTensor1FromPtr<3>(&*(N_face_edge(2, 0).data().begin())),
1260  getFTensor1FromPtr<3>(&*(N_face_edge(2, 1).data().begin())),
1261  getFTensor1FromPtr<3>(&*(N_face_edge(2, 2).data().begin())),
1262  getFTensor1FromPtr<3>(&*(N_face_edge(3, 0).data().begin())),
1263  getFTensor1FromPtr<3>(&*(N_face_edge(3, 1).data().begin())),
1264  getFTensor1FromPtr<3>(&*(N_face_edge(3, 2).data().begin()))};
1265  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_e[] = {
1266  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 0).data().begin())),
1267  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 1).data().begin())),
1268  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 2).data().begin())),
1269  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 0).data().begin())),
1270  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 1).data().begin())),
1271  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 2).data().begin())),
1272  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 0).data().begin())),
1273  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 1).data().begin())),
1274  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 2).data().begin())),
1275  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 0).data().begin())),
1276  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 1).data().begin())),
1277  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 2).data().begin()))};
1278 
1279  // face-face
1280  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_f_f[] = {
1281  getFTensor1FromPtr<3>(&*(N_face_bubble[0].data().begin())),
1282  getFTensor1FromPtr<3>(&*(N_face_bubble[1].data().begin())),
1283  getFTensor1FromPtr<3>(&*(N_face_bubble[2].data().begin())),
1284  getFTensor1FromPtr<3>(&*(N_face_bubble[3].data().begin()))};
1285  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_f[] = {
1286  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[0].data().begin())),
1287  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[1].data().begin())),
1288  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[2].data().begin())),
1289  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[3].data().begin()))};
1290 
1291  // volume-edge
1292  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_e[] = {
1293  getFTensor1FromPtr<3>(&*N_volume_edge[0].data().begin()),
1294  getFTensor1FromPtr<3>(&*N_volume_edge[1].data().begin()),
1295  getFTensor1FromPtr<3>(&*N_volume_edge[2].data().begin()),
1296  getFTensor1FromPtr<3>(&*N_volume_edge[3].data().begin()),
1297  getFTensor1FromPtr<3>(&*N_volume_edge[4].data().begin()),
1298  getFTensor1FromPtr<3>(&*N_volume_edge[5].data().begin())};
1299  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_e[] = {
1300  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[0].data().begin()),
1301  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[1].data().begin()),
1302  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[2].data().begin()),
1303  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[3].data().begin()),
1304  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[4].data().begin()),
1305  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[5].data().begin())};
1306 
1307  // volume-faces
1308  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_f[] = {
1309  getFTensor1FromPtr<3>(&*N_volume_face[0].data().begin()),
1310  getFTensor1FromPtr<3>(&*N_volume_face[1].data().begin()),
1311  getFTensor1FromPtr<3>(&*N_volume_face[2].data().begin()),
1312  getFTensor1FromPtr<3>(&*N_volume_face[3].data().begin())};
1313  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1314  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[0].data().begin()),
1315  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[1].data().begin()),
1316  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[2].data().begin()),
1317  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[3].data().begin())};
1318 
1319  // volume-bubble
1320  auto *base_vol_ptr = &*(N_volume_bubble.data().begin());
1321  auto *diff_base_vol_ptr = &*(diffN_volume_bubble.data().begin());
1322  auto t_base_v = getFTensor1FromPtr<3>(base_vol_ptr);
1323  auto t_diff_base_v = getFTensor2HVecFromPtr<3, 3>(diff_base_vol_ptr);
1324 
1325  int count_dofs = 0;
1326  int count_dofs_face = 0;
1327  int count_dofs_volume = 0;
1328 
1329  auto max_volume_order =
1330  std::max(volume_order,
1331  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(volume_order));
1332  max_volume_order =
1333  std::max(max_volume_order,
1334  AinsworthOrderHooks::broken_nbfacetri_face_hdiv(volume_order));
1335  max_volume_order =
1336  std::max(max_volume_order,
1337  AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order));
1338  max_volume_order =
1339  std::max(max_volume_order,
1340  AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order));
1341  max_volume_order = std::max(
1342  max_volume_order,
1343  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order));
1344 
1345  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1346  for (int oo = 0; oo < max_volume_order; oo++) {
1347 
1348  // faces-edge (((P) > 0) ? (P) : 0)
1349  if (oo < AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(volume_order))
1350  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
1351  dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1352  for (auto ff = 0; ff != 4; ++ff) {
1353  for (int ee = 0; ee != 3; ++ee) {
1354  t_base(i) = t_base_f_e[ff * 3 + ee](i);
1355  ++t_base;
1356  ++t_base_f_e[ff * 3 + ee];
1357  ++count_dofs;
1358  ++count_dofs_face;
1359  }
1360  for (int ee = 0; ee != 3; ++ee) {
1361  t_diff_base(i, j) = t_diff_base_f_e[ff * 3 + ee](i, j);
1362  ++t_diff_base;
1363  ++t_diff_base_f_e[ff * 3 + ee];
1364  }
1365  }
1366  }
1367 
1368  // face-face (P - 1) * (P - 2) / 2
1369  if (oo < AinsworthOrderHooks::broken_nbfacetri_face_hdiv(volume_order))
1370  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
1371  dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1372  for (auto ff = 0; ff != 4; ++ff) {
1373  t_base(i) = t_base_f_f[ff](i);
1374  ++t_base;
1375  ++t_base_f_f[ff];
1376  t_diff_base(i, j) = t_diff_base_f_f[ff](i, j);
1377  ++t_diff_base;
1378  ++t_diff_base_f_f[ff];
1379  ++count_dofs;
1380  ++count_dofs_face;
1381  }
1382  }
1383 
1384  // volume-edge (P - 1)
1385  if (oo < AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order))
1386  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
1387  dd != NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1388  for (int ee = 0; ee < 6; ++ee) {
1389  t_base(i) = t_base_v_e[ee](i);
1390  ++t_base;
1391  ++t_base_v_e[ee];
1392  t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
1393  ++t_diff_base;
1394  ++t_diff_base_v_e[ee];
1395  ++count_dofs;
1396  ++count_dofs_volume;
1397  }
1398  }
1399 
1400  // volume-face (P - 1) * (P - 2)
1401  if (oo < AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order))
1402  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
1403  dd != NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1404  for (int ff = 0; ff < 4; ff++) {
1405  t_base(i) = t_base_v_f[ff](i);
1406  ++t_base;
1407  ++t_base_v_f[ff];
1408  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1409  ++t_diff_base;
1410  ++t_diff_base_v_f[ff];
1411  ++count_dofs;
1412  ++count_dofs_volume;
1413  }
1414  }
1415 
1416  // volume-bubble (P - 3) * (P - 2) * (P - 1) / 2
1417  if (oo <
1418  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order))
1419  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
1420  dd != NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
1421  t_base(i) = t_base_v(i);
1422  ++t_base;
1423  ++t_base_v;
1424  t_diff_base(i, j) = t_diff_base_v(i, j);
1425  ++t_diff_base;
1426  ++t_diff_base_v;
1427  ++count_dofs;
1428  ++count_dofs_volume;
1429  }
1430  }
1431  }
1432 
1433 #ifndef NDEBUG
1434  if (nb_dofs != count_dofs / nb_gauss_pts) {
1435  MOFEM_LOG_CHANNEL("SELF");
1436  MOFEM_LOG("SELF", Sev::error) << "Nb dofs face: " << 4 * nb_dofs_face
1437  << " -> " << count_dofs_face / nb_gauss_pts;
1438  MOFEM_LOG("SELF", Sev::error) << "Nb dofs volume: " << nb_dofs_volume
1439  << " -> " << count_dofs_volume / nb_gauss_pts;
1440  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1441  "Number of dofs %d is different than expected %d",
1442  count_dofs / nb_gauss_pts, nb_dofs);
1443  }
1444 #endif // NDEBUG
1445 
1446  if (interior_cache_ptr) {
1447  auto p = interior_cache_ptr->emplace(
1448  TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1449  p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1450  p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1451  }
1452 
1453  MoFEMFunctionReturn(0);
1454 }

◆ getValueHdivDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBase ( MatrixDouble &  pts)
protected

Definition at line 1456 of file TetPolynomialBase.cpp.

1456  {
1457  MoFEMFunctionBegin;
1458 
1459  EntitiesFieldData &data = cTx->dAta;
1460  const FieldApproximationBase base = cTx->bAse;
1461  if (base != DEMKOWICZ_JACOBI_BASE) {
1462  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1463  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1464  "but base is %s",
1465  ApproximationBaseNames[base]);
1466  }
1467  int nb_gauss_pts = pts.size2();
1468 
1469  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1470 
1471  int p_f[4];
1472  double *phi_f[4];
1473  double *diff_phi_f[4];
1474 
1475  auto get_face_cache_ptr = [this]() -> TetBaseCache::HDivBaseFaceCacheMI * {
1476  if (vPtr) {
1477  auto it = TetBaseCache::hDivBaseFace[DEMKOWICZ_JACOBI_BASE].find(vPtr);
1478  if (it != TetBaseCache::hDivBaseFace[DEMKOWICZ_JACOBI_BASE].end()) {
1479  return &it->second;
1480  }
1481  }
1482  return nullptr;
1483  };
1484 
1485  auto face_cache_ptr = get_face_cache_ptr();
1486 
1487  // Calculate base function on tet faces
1488  for (int ff = 0; ff != 4; ff++) {
1489  int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1490  int order = volume_order > face_order ? volume_order : face_order;
1491  data.dataOnEntities[MBTRI][ff].getN(base).resize(
1492  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1493  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1494  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1495  if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
1496  continue;
1497 
1498  if (face_cache_ptr) {
1499  auto it = face_cache_ptr->find(boost::make_tuple(
1500 
1501  face_order, nb_gauss_pts,
1502 
1503  data.facesNodes(ff, 0), data.facesNodes(ff, 1), data.facesNodes(ff, 2)
1504 
1505  ));
1506  if (it != face_cache_ptr->end()) {
1507  noalias(data.dataOnEntities[MBTRI][ff].getN(base)) = it->N;
1508  noalias(data.dataOnEntities[MBTRI][ff].getDiffN(base)) = it->diffN;
1509  continue;
1510  }
1511  }
1512 
1513  p_f[ff] = order;
1514  phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1515  diff_phi_f[ff] =
1516  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1517 
1518  CHKERR Hdiv_Demkowicz_Face_MBTET_ON_FACE(
1519  &data.facesNodes(ff, 0), order,
1520  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1521  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1522  phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1523  if (face_cache_ptr) {
1524  auto p = face_cache_ptr->emplace(TetBaseCache::HDivBaseCacheItem{
1525  face_order, nb_gauss_pts, data.facesNodes(ff, 0),
1526  data.facesNodes(ff, 1), data.facesNodes(ff, 2)});
1527  p.first->N = data.dataOnEntities[MBTRI][ff].getN(base);
1528  p.first->diffN = data.dataOnEntities[MBTRI][ff].getDiffN(base);
1529  }
1530  }
1531 
1532  auto get_interior_cache = [this]() -> TetBaseCache::BaseCacheMI * {
1533  if (vPtr) {
1534  auto it =
1535  TetBaseCache::hdivBaseInterior[DEMKOWICZ_JACOBI_BASE].find(vPtr);
1536  if (it != TetBaseCache::hdivBaseInterior[DEMKOWICZ_JACOBI_BASE].end()) {
1537  return &it->second;
1538  }
1539  }
1540  return nullptr;
1541  };
1542 
1543  auto interior_cache_ptr = get_interior_cache();
1544 
1545  // Calculate base functions in tet interior
1546  if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
1547  data.dataOnEntities[MBTET][0].getN(base).resize(
1548  nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1549  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
1550  nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1551 
1552  for (int v = 0; v != 1; ++v) {
1553  if (interior_cache_ptr) {
1554  auto it = interior_cache_ptr->find(
1555  boost::make_tuple(volume_order, nb_gauss_pts));
1556  if (it != interior_cache_ptr->end()) {
1557  noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1558  noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1559  continue;
1560  }
1561  }
1562 
1563  double *phi_v = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1564  double *diff_phi_v =
1565  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1566 
1567  CHKERR Hdiv_Demkowicz_Interior_MBTET(
1568  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1569  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
1570  diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
1571  if (interior_cache_ptr) {
1572  auto p = interior_cache_ptr->emplace(
1573  TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1574  p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1575  p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1576  }
1577  }
1578  }
1579 
1580  // Set size of face base correctly
1581  for (int ff = 0; ff != 4; ff++) {
1582  int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1583  data.dataOnEntities[MBTRI][ff].getN(base).resize(
1584  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1585  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1586  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1587  }
1588 
1589  MoFEMFunctionReturn(0);
1590 }

◆ getValueHdivDemkowiczBrokenBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBrokenBase ( MatrixDouble &  pts)
protected

Definition at line 1593 of file TetPolynomialBase.cpp.

1593  {
1594  MoFEMFunctionBegin;
1595 
1596  EntitiesFieldData &data = cTx->dAta;
1597  const FieldApproximationBase base = cTx->bAse;
1598  if (base != DEMKOWICZ_JACOBI_BASE) {
1599  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1600  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1601  "but base is %s",
1602  ApproximationBaseNames[base]);
1603  }
1604  int nb_gauss_pts = pts.size2();
1605 
1606  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1607  int nb_dofs_face = NBFACETRI_DEMKOWICZ_HDIV(volume_order);
1608  int nb_dofs_volume = NBVOLUMETET_DEMKOWICZ_HDIV(volume_order);
1609  int nb_dofs = 4 * nb_dofs_face + nb_dofs_volume;
1610  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
1611  false);
1612  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 9 * nb_dofs,
1613  false);
1614  if (nb_dofs == 0)
1615  MoFEMFunctionReturnHot(0);
1616 
1617  auto get_interior_cache = [this]() -> TetBaseCache::BaseCacheMI * {
1618  if (vPtr) {
1619  auto it =
1620  TetBaseCache::hdivBrokenBaseInterior[DEMKOWICZ_JACOBI_BASE].find(
1621  vPtr);
1622  if (it !=
1623  TetBaseCache::hdivBrokenBaseInterior[DEMKOWICZ_JACOBI_BASE].end()) {
1624  return &it->second;
1625  }
1626  }
1627  return nullptr;
1628  };
1629 
1630  auto interior_cache_ptr = get_interior_cache();
1631 
1632  if (interior_cache_ptr) {
1633  auto it =
1634  interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
1635  if (it != interior_cache_ptr->end()) {
1636  noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1637  noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1638  MoFEMFunctionReturnHot(0);
1639  }
1640  }
1641 
1642  std::array<MatrixDouble, 4> face_base_fun{
1643  MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1644  MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1645  MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1646  MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face)};
1647  std::array<MatrixDouble, 4> face_diff_base{
1648  MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1649  MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1650  MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1651  MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face)};
1652 
1653  int faces_nodes[4][3] = {{0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
1654 
1655  std::array<int, 4> p_f{volume_order, volume_order, volume_order,
1656  volume_order};
1657  std::array<double *, 4> phi_f{
1658  &*face_base_fun[0].data().begin(), &*face_base_fun[1].data().begin(),
1659  &*face_base_fun[2].data().begin(), &*face_base_fun[3].data().begin()};
1660  std::array<double *, 4> diff_phi_f{
1661  &*face_diff_base[0].data().begin(), &*face_diff_base[1].data().begin(),
1662  &*face_diff_base[2].data().begin(), &*face_diff_base[3].data().begin()};
1663 
1664  // Calculate base function on tet faces
1665  for (int ff = 0; ff != 4; ff++) {
1666  CHKERR Hdiv_Demkowicz_Face_MBTET_ON_FACE(
1667  // &data.facesNodes(ff, 0)
1668  faces_nodes[ff], p_f[ff],
1669  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1670  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1671  phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1672  }
1673 
1674  MatrixDouble vol_bases(nb_gauss_pts, 3 * nb_dofs_volume);
1675  MatrixDouble vol_diff_bases(nb_gauss_pts, 9 * nb_dofs_volume);
1676  auto *phi_v = &*vol_bases.data().begin();
1677  auto *diff_phi_v = &*vol_diff_bases.data().begin();
1678  CHKERR Hdiv_Demkowicz_Interior_MBTET(
1679  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1680  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f.data(),
1681  phi_f.data(), diff_phi_f.data(), phi_v, diff_phi_v, nb_gauss_pts);
1682 
1683  // faces
1684  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_f[] = {
1685  getFTensor1FromPtr<3>(phi_f[0]), getFTensor1FromPtr<3>(phi_f[1]),
1686  getFTensor1FromPtr<3>(phi_f[2]), getFTensor1FromPtr<3>(phi_f[3])};
1687  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1688  getFTensor2HVecFromPtr<3, 3>(diff_phi_f[0]),
1689  getFTensor2HVecFromPtr<3, 3>(diff_phi_f[1]),
1690  getFTensor2HVecFromPtr<3, 3>(diff_phi_f[2]),
1691  getFTensor2HVecFromPtr<3, 3>(diff_phi_f[3])};
1692 
1693  // volumes
1694  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v =
1695  getFTensor1FromPtr<3>(&*vol_bases.data().begin());
1696  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v =
1697  getFTensor2HVecFromPtr<3, 3>(&*vol_diff_bases.data().begin());
1698 
1699  auto t_base = getFTensor1FromPtr<3>(
1700  &*data.dataOnEntities[MBTET][0].getN(base).data().begin());
1701  auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(
1702  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin());
1703 
1704  FTENSOR_INDEX(3, i);
1705  FTENSOR_INDEX(3, j);
1706 
1707  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1708  for (int oo = 0; oo < volume_order; oo++) {
1709  // face
1710  for (auto dd = NBFACETRI_DEMKOWICZ_HDIV(oo);
1711  dd != NBFACETRI_DEMKOWICZ_HDIV(oo + 1); ++dd) {
1712  for (auto ff = 0; ff != 4; ++ff) {
1713  t_base(i) = t_base_v_f[ff](i);
1714  ++t_base;
1715  ++t_base_v_f[ff];
1716  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1717  ++t_diff_base;
1718  ++t_diff_base_v_f[ff];
1719  }
1720  }
1721  // volume
1722  for (auto dd = NBVOLUMETET_DEMKOWICZ_HDIV(oo);
1723  dd != NBVOLUMETET_DEMKOWICZ_HDIV(oo + 1); ++dd) {
1724  t_base(i) = t_base_v(i);
1725  ++t_base;
1726  ++t_base_v;
1727  t_diff_base(i, j) = t_diff_base_v(i, j);
1728  ++t_diff_base;
1729  ++t_diff_base_v;
1730  }
1731  }
1732  }
1733 
1734  if (interior_cache_ptr) {
1735  auto p = interior_cache_ptr->emplace(
1736  TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1737  p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1738  p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1739  }
1740 
1741  MoFEMFunctionReturn(0);
1742 }

◆ getValueL2()

MoFEMErrorCode TetPolynomialBase::getValueL2 ( MatrixDouble &  pts)
protected

Get base functions for L2 space.

Parameters
ptsmatrix of integration pts
Returns
MoFEMErrorCode
Note
matrix of integration points on rows has local coordinates of finite element on columns are integration pts.

Definition at line 525 of file TetPolynomialBase.cpp.

525  {
526  MoFEMFunctionBegin;
527 
528  switch (cTx->bAse) {
529  case AINSWORTH_LEGENDRE_BASE:
530  case AINSWORTH_LOBATTO_BASE:
531  case DEMKOWICZ_JACOBI_BASE:
532  CHKERR getValueL2AinsworthBase(pts);
533  break;
534  case AINSWORTH_BERNSTEIN_BEZIER_BASE:
535  CHKERR getValueL2BernsteinBezierBase(pts);
536  break;
537  default:
538  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
539  }
540 
541  MoFEMFunctionReturn(0);
542 }

◆ getValueL2AinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueL2AinsworthBase ( MatrixDouble &  pts)
protected

Definition at line 544 of file TetPolynomialBase.cpp.

544  {
545  MoFEMFunctionBegin;
546 
547  EntitiesFieldData &data = cTx->dAta;
548  const FieldApproximationBase base = cTx->bAse;
549  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
550  double *diffL, const int dim) =
551  cTx->basePolynomialsType0;
552 
553  int nb_gauss_pts = pts.size2();
554  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
555  int nb_dofs = NBVOLUMETET_L2(volume_order);
556 
557  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_dofs, false);
558  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 3 * nb_dofs,
559  false);
560 
561  if (!nb_dofs)
562  MoFEMFunctionReturnHot(0);
563 
564  auto get_interior_cache = [this](auto base) -> TetBaseCache::BaseCacheMI * {
565  if (vPtr) {
566  auto it = TetBaseCache::l2BaseInterior[base].find(vPtr);
567  if (it != TetBaseCache::l2BaseInterior[base].end()) {
568  return &it->second;
569  }
570  }
571  return nullptr;
572  };
573 
574  auto interior_cache_ptr = get_interior_cache(base);
575 
576  if (interior_cache_ptr) {
577  auto it =
578  interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
579  if (it != interior_cache_ptr->end()) {
580  noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
581  noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
582  MoFEMFunctionReturnHot(0);
583  }
584  }
585 
586  CHKERR L2_Ainsworth_ShapeFunctions_MBTET(
587  data.dataOnEntities[MBTET][0].getOrder(),
588  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
589  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
590  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
591  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
592  nb_gauss_pts, base_polynomials);
593 
594  if (interior_cache_ptr) {
595  auto p = interior_cache_ptr->emplace(
596  TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
597  p.first->N = data.dataOnEntities[MBTET][0].getN(base);
598  p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
599  }
600 
601  MoFEMFunctionReturn(0);
602 }

◆ getValueL2BernsteinBezierBase()

MoFEMErrorCode TetPolynomialBase::getValueL2BernsteinBezierBase ( MatrixDouble &  pts)
protected

Definition at line 605 of file TetPolynomialBase.cpp.

605  {
606  MoFEMFunctionBegin;
607 
608  EntitiesFieldData &data = cTx->dAta;
609  const std::string field_name = cTx->fieldName;
610  const int nb_gauss_pts = pts.size2();
611 
612  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
613  (unsigned int)nb_gauss_pts)
614  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
615  "Base functions or nodes has wrong number of integration points "
616  "for base %s",
617  ApproximationBaseNames[NOBASE]);
618  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
619 
620  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
621  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
622  if (!ptr)
623  ptr.reset(new MatrixInt());
624  return *ptr;
625  };
626 
627  auto get_base = [field_name](auto &data) -> MatrixDouble & {
628  auto &ptr = data.getBBNSharedPtr(field_name);
629  if (!ptr)
630  ptr.reset(new MatrixDouble());
631  return *ptr;
632  };
633 
634  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
635  auto &ptr = data.getBBDiffNSharedPtr(field_name);
636  if (!ptr)
637  ptr.reset(new MatrixDouble());
638  return *ptr;
639  };
640 
641  auto get_alpha_by_name_ptr =
642  [](auto &data,
643  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
644  return data.getBBAlphaIndicesSharedPtr(field_name);
645  };
646 
647  auto get_base_by_name_ptr =
648  [](auto &data,
649  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
650  return data.getBBNSharedPtr(field_name);
651  };
652 
653  auto get_diff_base_by_name_ptr =
654  [](auto &data,
655  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
656  return data.getBBDiffNSharedPtr(field_name);
657  };
658 
659  auto get_alpha_by_order_ptr =
660  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
661  return data.getBBAlphaIndicesByOrderSharedPtr(o);
662  };
663 
664  auto get_base_by_order_ptr =
665  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
666  return data.getBBNByOrderSharedPtr(o);
667  };
668 
669  auto get_diff_base_by_order_ptr =
670  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
671  return data.getBBDiffNByOrderSharedPtr(o);
672  };
673 
674  if (data.spacesOnEntities[MBTET].test(L2)) {
675  if (data.dataOnEntities[MBTET].size() != 1)
676  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
677  "Wrong size ent of ent data");
678 
679  auto &ent_data = data.dataOnEntities[MBTET][0];
680  const int order = ent_data.getOrder();
681  const int nb_dofs = NBVOLUMETET_L2(order);
682 
683  if (get_alpha_by_order_ptr(ent_data, order)) {
684  get_alpha_by_name_ptr(ent_data, field_name) =
685  get_alpha_by_order_ptr(ent_data, order);
686  get_base_by_name_ptr(ent_data, field_name) =
687  get_base_by_order_ptr(ent_data, order);
688  get_diff_base_by_name_ptr(ent_data, field_name) =
689  get_diff_base_by_order_ptr(ent_data, order);
690  } else {
691 
692  auto &get_n = get_base(ent_data);
693  auto &get_diff_n = get_diff_base(ent_data);
694  get_n.resize(nb_gauss_pts, nb_dofs, false);
695  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
696 
697  if (nb_dofs) {
698 
699  if (order == 0) {
700 
701  if (nb_dofs != 1)
702  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
703  "Inconsistent number of DOFs");
704 
705  auto &tri_alpha = get_alpha(ent_data);
706  tri_alpha.clear();
707  get_n(0, 0) = 1;
708  get_diff_n.clear();
709 
710  } else {
711 
712  if (nb_dofs != 4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order) +
713  NBVOLUMETET_H1(order) ||
714  nb_dofs != NBVOLUMETET_L2(order))
715  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
716  "Inconsistent number of DOFs");
717 
718  auto &tet_alpha = get_alpha(ent_data);
719  tet_alpha.resize(nb_dofs, 4, false);
720 
721  CHKERR BernsteinBezier::generateIndicesVertexTet(order,
722  &tet_alpha(0, 0));
723  if (order > 1) {
724  std::array<int, 6> edge_n{order, order, order, order, order, order};
725  std::array<int *, 6> tet_edge_ptr{
726  &tet_alpha(4, 0),
727  &tet_alpha(4 + 1 * NBEDGE_H1(order), 0),
728  &tet_alpha(4 + 2 * NBEDGE_H1(order), 0),
729  &tet_alpha(4 + 3 * NBEDGE_H1(order), 0),
730  &tet_alpha(4 + 4 * NBEDGE_H1(order), 0),
731  &tet_alpha(4 + 5 * NBEDGE_H1(order), 0)};
732  CHKERR BernsteinBezier::generateIndicesEdgeTet(edge_n.data(),
733  tet_edge_ptr.data());
734  if (order > 2) {
735  std::array<int, 6> face_n{order, order, order, order};
736  std::array<int *, 6> tet_face_ptr{
737  &tet_alpha(4 + 6 * NBEDGE_H1(order), 0),
738  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 1 * NBFACETRI_H1(order),
739  0),
740  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 2 * NBFACETRI_H1(order),
741  0),
742  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 3 * NBFACETRI_H1(order),
743  0),
744  };
745  CHKERR BernsteinBezier::generateIndicesTriTet(
746  face_n.data(), tet_face_ptr.data());
747  if (order > 3)
748  CHKERR BernsteinBezier::generateIndicesTetTet(
749  order,
750  &tet_alpha(
751  4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order), 0));
752  }
753  }
754 
755  CHKERR BernsteinBezier::baseFunctionsTet(
756  order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
757  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
758  &get_diff_n(0, 0));
759 
760  get_alpha_by_order_ptr(ent_data, order) =
761  get_alpha_by_name_ptr(ent_data, field_name);
762  get_base_by_order_ptr(ent_data, order) =
763  get_base_by_name_ptr(ent_data, field_name);
764  get_diff_base_by_order_ptr(ent_data, order) =
765  get_diff_base_by_name_ptr(ent_data, field_name);
766  }
767  }
768  }
769  } else {
770  auto &ent_data = data.dataOnEntities[MBTET][0];
771  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
772  auto &get_n = get_base(ent_data);
773  auto &get_diff_n = get_diff_base(ent_data);
774  get_n.resize(nb_gauss_pts, 0, false);
775  get_diff_n.resize(nb_gauss_pts, 0, false);
776  }
777 
778  MoFEMFunctionReturn(0);
779 }

◆ query_interface()

MoFEMErrorCode TetPolynomialBase::query_interface ( boost::typeindex::type_index  type_index,
UnknownInterface **  iface 
) const
virtual

Reimplemented from MoFEM::BaseFunction.

Definition at line 87 of file TetPolynomialBase.cpp.

88  {
89 
90  MoFEMFunctionBeginHot;
91  *iface = const_cast<TetPolynomialBase *>(this);
92  MoFEMFunctionReturnHot(0);
93 }

◆ setDofsSideMap()

MoFEMErrorCode TetPolynomialBase::setDofsSideMap ( const FieldSpace  space,
const FieldContinuity  continuity,
const FieldApproximationBase  base,
DofsSideMap &  dofs_side_map 
)
static

Set map of dof to side number.

That is used for broken space to establish connection between dofs in the interior of element/entity and side of element/entity to which that dof is associated. That depends on implementation of the base for given space, and has to be implemented while implementing base function for given space.

Parameters
space
continuity
base
DofsSideMap
Returns
MoFEMErrorCode

Definition at line 2129 of file TetPolynomialBase.cpp.

2134  {
2135  MoFEMFunctionBegin;
2136 
2137  switch (continuity) {
2138  case DISCONTINUOUS:
2139 
2140  switch (space) {
2141  case HDIV:
2142  CHKERR setDofsSideMapHdiv(continuity, base, dofs_side_map);
2143  break;
2144  default:
2145  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2146  FieldSpaceNames[space]);
2147  }
2148  break;
2149 
2150  default:
2151  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2152  "Unknown (or not implemented) continuity");
2153  }
2154 
2155  MoFEMFunctionReturn(0);
2156 }

◆ setDofsSideMapHdiv()

MoFEMErrorCode TetPolynomialBase::setDofsSideMapHdiv ( const FieldContinuity  continuity,
const FieldApproximationBase  base,
DofsSideMap &  dofs_side_map 
)
staticprotected

Set the Dofs Side Map Hdiv object.

Parameters
space
continuity
base
dofs_side_map
Returns
MoFEMErrorCode

Definition at line 2159 of file TetPolynomialBase.cpp.

2161  {
2162  MoFEMFunctionBegin;
2163 
2164  // That has to be consistent with implementation of getValueHdiv for
2165  // particular base functions.
2166 
2167  auto set_ainsworth = [&dofs_side_map]() {
2168  MoFEMFunctionBegin;
2169 
2170  dofs_side_map.clear();
2171 
2172  int dof = 0;
2173  for (int oo = 0; oo < Field::maxBrokenDofsOrder; oo++) {
2174 
2175  // faces-edge (((P) > 0) ? (P) : 0)
2176  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
2177  dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
2178  for (auto ff = 0; ff != 4; ++ff) {
2179  for (int ee = 0; ee != 3; ++ee) {
2180  dofs_side_map.insert(DofsSideMapData{MBTRI, ff, dof});
2181  ++dof;
2182  }
2183  }
2184  }
2185 
2186  // face-face (P - 1) * (P - 2) / 2
2187  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
2188  dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
2189  for (auto ff = 0; ff != 4; ++ff) {
2190  dofs_side_map.insert(DofsSideMapData{MBTRI, ff, dof});
2191  ++dof;
2192  }
2193  }
2194 
2195  // volume-edge (P - 1)
2196  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
2197  dd != NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
2198  for (int ee = 0; ee < 6; ee++) {
2199  dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2200  ++dof;
2201  }
2202  }
2203  // volume-face (P - 1) * (P - 2)
2204  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
2205  dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
2206  for (int ff = 0; ff < 4; ff++) {
2207  dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2208  ++dof;
2209  }
2210  }
2211  // volume-bubble (P - 3) * (P - 2) * (P - 1) / 2
2212  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
2213  dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
2214  dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2215  ++dof;
2216  }
2217  }
2218 
2219  MoFEMFunctionReturn(0);
2220  };
2221 
2222  auto set_demkowicz = [&dofs_side_map]() {
2223  MoFEMFunctionBegin;
2224 
2225  dofs_side_map.clear();
2226 
2227  int dof = 0;
2228  for (int oo = 0; oo < Field::maxBrokenDofsOrder; oo++) {
2229 
2230  // face
2231  for (auto dd = NBFACETRI_DEMKOWICZ_HDIV(oo);
2232  dd != NBFACETRI_DEMKOWICZ_HDIV(oo + 1); ++dd) {
2233  for (auto ff = 0; ff != 4; ++ff) {
2234  dofs_side_map.insert(DofsSideMapData{MBTRI, ff, dof});
2235  ++dof;
2236  }
2237  }
2238  // volume
2239  for (auto dd = NBVOLUMETET_DEMKOWICZ_HDIV(oo);
2240  dd != NBVOLUMETET_DEMKOWICZ_HDIV(oo + 1); ++dd) {
2241  dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2242  ++dof;
2243  }
2244  }
2245 
2246  MoFEMFunctionReturn(0);
2247  };
2248 
2249  switch (continuity) {
2250  case DISCONTINUOUS:
2251  switch (base) {
2252  case AINSWORTH_LEGENDRE_BASE:
2253  case AINSWORTH_LOBATTO_BASE:
2254  return set_ainsworth();
2255  case DEMKOWICZ_JACOBI_BASE:
2256  return set_demkowicz();
2257  default:
2258  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
2259  }
2260  break;
2261 
2262  default:
2263  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2264  "Unknown (or not implemented) continuity");
2265  }
2266 
2267  MoFEMFunctionReturn(0);
2268 }

◆ switchCacheBaseFace() [1/3]

template<int SPACE>
static bool MoFEM::TetPolynomialBase::switchCacheBaseFace ( FieldApproximationBase  base,
void *  ptr 
)
static

◆ switchCacheBaseFace() [2/3]

template<>
bool MoFEM::TetPolynomialBase::switchCacheBaseFace ( FieldApproximationBase  base,
void *  ptr 
)
static

◆ switchCacheBaseFace() [3/3]

template<>
bool MoFEM::TetPolynomialBase::switchCacheBaseFace ( FieldApproximationBase  base,
void *  ptr 
)
static

Definition at line 2289 of file TetPolynomialBase.cpp.

2290  {
2291  return tetCacheSwitch(ptr, TetBaseCache::hDivBaseFace[base],
2292  std::string("hDivBaseFace") +
2293  ApproximationBaseNames[base]);
2294 }

◆ switchCacheBaseInterior() [1/5]

template<int SPACE>
static bool MoFEM::TetPolynomialBase::switchCacheBaseInterior ( FieldApproximationBase  base,
void *  ptr 
)
static

◆ switchCacheBaseInterior() [2/5]

template<>
bool MoFEM::TetPolynomialBase::switchCacheBaseInterior ( FieldApproximationBase  base,
void *  ptr 
)
static

◆ switchCacheBaseInterior() [3/5]

template<>
bool MoFEM::TetPolynomialBase::switchCacheBaseInterior ( FieldApproximationBase  base,
void *  ptr 
)
static

◆ switchCacheBaseInterior() [4/5]

template<>
bool MoFEM::TetPolynomialBase::switchCacheBaseInterior ( FieldApproximationBase  base,
void *  ptr 
)
static

Definition at line 2297 of file TetPolynomialBase.cpp.

2298  {
2299  return tetCacheSwitch(ptr, TetBaseCache::hdivBaseInterior[base],
2300  std::string("hdivBaseInterior") +
2301  ApproximationBaseNames[base]);
2302 }

◆ switchCacheBaseInterior() [5/5]

template<>
bool MoFEM::TetPolynomialBase::switchCacheBaseInterior ( FieldApproximationBase  base,
void *  ptr 
)
static

Definition at line 2359 of file TetPolynomialBase.cpp.

2360  {
2361  return tetCacheSwitch(ptr, TetBaseCache::l2BaseInterior[base],
2362  std::string("hdivBaseInterior") +
2363  ApproximationBaseNames[base]);
2364 }

◆ switchCacheBaseOff() [1/10]

template<int SPACE>
static void MoFEM::TetPolynomialBase::switchCacheBaseOff ( FieldApproximationBase  base,
std::vector< void * >  v 
)
static

◆ switchCacheBaseOff() [2/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOff ( FieldApproximationBase  base,
std::vector< void * >  v 
)
static

◆ switchCacheBaseOff() [3/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOff ( FieldApproximationBase  base,
std::vector< void * >  v 
)
static

◆ switchCacheBaseOff() [4/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOff ( FieldApproximationBase  base,
std::vector< void * >  v 
)
static

Definition at line 2329 of file TetPolynomialBase.cpp.

2330  {
2331  for (auto fe_ptr : v) {
2332  if (TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr)) {
2333  TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr);
2334  }
2335  if (TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr)) {
2336  TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr);
2337  }
2338  if (TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr)) {
2339  TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr);
2340  }
2341  }
2342 }

◆ switchCacheBaseOff() [5/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOff ( FieldApproximationBase  base,
std::vector< void * >  v 
)
static

Definition at line 2377 of file TetPolynomialBase.cpp.

2378  {
2379  for (auto fe_ptr : v) {
2380  if (TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr)) {
2381  TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr);
2382  }
2383  }
2384 }

◆ switchCacheBaseOff() [6/10]

template<int SPACE>
static void MoFEM::TetPolynomialBase::switchCacheBaseOff ( std::vector< void * >  v)
static

◆ switchCacheBaseOff() [7/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOff ( std::vector< void * >  v)
static

◆ switchCacheBaseOff() [8/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOff ( std::vector< void * >  v)
static

◆ switchCacheBaseOff() [9/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOff ( std::vector< void * >  v)
static

Definition at line 2352 of file TetPolynomialBase.cpp.

2352  {
2353  for (auto b = 0; b != LASTBASE; ++b) {
2354  switchCacheBaseOff<HDIV>(static_cast<FieldApproximationBase>(b), v);
2355  }
2356 }

◆ switchCacheBaseOff() [10/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOff ( std::vector< void * >  v)
static

Definition at line 2394 of file TetPolynomialBase.cpp.

2394  {
2395  for (auto b = 0; b != LASTBASE; ++b) {
2396  switchCacheBaseOff<L2>(static_cast<FieldApproximationBase>(b), v);
2397  }
2398 }

◆ switchCacheBaseOn() [1/10]

template<int SPACE>
static void MoFEM::TetPolynomialBase::switchCacheBaseOn ( FieldApproximationBase  base,
std::vector< void * >  v 
)
static

◆ switchCacheBaseOn() [2/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOn ( FieldApproximationBase  base,
std::vector< void * >  v 
)
static

◆ switchCacheBaseOn() [3/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOn ( FieldApproximationBase  base,
std::vector< void * >  v 
)
static

◆ switchCacheBaseOn() [4/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOn ( FieldApproximationBase  base,
std::vector< void * >  v 
)
static

Definition at line 2313 of file TetPolynomialBase.cpp.

2314  {
2315  for (auto fe_ptr : v) {
2316  if (!TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr)) {
2317  TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr);
2318  }
2319  if (!TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr)) {
2320  TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr);
2321  }
2322  if (!TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr)) {
2323  TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr);
2324  }
2325  }
2326 }

◆ switchCacheBaseOn() [5/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOn ( FieldApproximationBase  base,
std::vector< void * >  v 
)
static

Definition at line 2367 of file TetPolynomialBase.cpp.

2368  {
2369  for (auto fe_ptr : v) {
2370  if (!TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr)) {
2371  TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr);
2372  }
2373  }
2374 }

◆ switchCacheBaseOn() [6/10]

template<int SPACE>
static void MoFEM::TetPolynomialBase::switchCacheBaseOn ( std::vector< void * >  v)
static

◆ switchCacheBaseOn() [7/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOn ( std::vector< void * >  v)
static

◆ switchCacheBaseOn() [8/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOn ( std::vector< void * >  v)
static

◆ switchCacheBaseOn() [9/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOn ( std::vector< void * >  v)
static

Definition at line 2345 of file TetPolynomialBase.cpp.

2345  {
2346  for (auto b = 0; b != LASTBASE; ++b) {
2347  switchCacheBaseOn<HDIV>(static_cast<FieldApproximationBase>(b), v);
2348  }
2349 }

◆ switchCacheBaseOn() [10/10]

template<>
void MoFEM::TetPolynomialBase::switchCacheBaseOn ( std::vector< void * >  v)
static

Definition at line 2387 of file TetPolynomialBase.cpp.

2387  {
2388  for (auto b = 0; b != LASTBASE; ++b) {
2389  switchCacheBaseOn<L2>(static_cast<FieldApproximationBase>(b), v);
2390  }
2391 }

◆ switchCacheBrokenBaseInterior() [1/3]

template<int SPACE>
static bool MoFEM::TetPolynomialBase::switchCacheBrokenBaseInterior ( FieldApproximationBase  base,
void *  ptr 
)
static

◆ switchCacheBrokenBaseInterior() [2/3]

template<>
bool MoFEM::TetPolynomialBase::switchCacheBrokenBaseInterior ( FieldApproximationBase  base,
void *  ptr 
)
static

◆ switchCacheBrokenBaseInterior() [3/3]

template<>
bool MoFEM::TetPolynomialBase::switchCacheBrokenBaseInterior ( FieldApproximationBase  base,
void *  ptr 
)
static

Definition at line 2305 of file TetPolynomialBase.cpp.

2306  {
2307  return tetCacheSwitch(ptr, TetBaseCache::hdivBrokenBaseInterior[base],
2308  std::string("hdivBrokenBaseInterior") +
2309  ApproximationBaseNames[base]);
2310 }

Member Data Documentation

◆ cTx

EntPolynomialBaseCtx* MoFEM::TetPolynomialBase::cTx
protected

Definition at line 71 of file TetPolynomialBase.hpp.

◆ diffN_face_bubble

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

Definition at line 171 of file TetPolynomialBase.hpp.

◆ diffN_face_edge

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

Definition at line 170 of file TetPolynomialBase.hpp.

◆ diffN_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::diffN_volume_bubble
protected

Definition at line 174 of file TetPolynomialBase.hpp.

◆ diffN_volume_edge

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

Definition at line 172 of file TetPolynomialBase.hpp.

◆ diffN_volume_face

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

Definition at line 173 of file TetPolynomialBase.hpp.

◆ N_face_bubble

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

Definition at line 165 of file TetPolynomialBase.hpp.

◆ N_face_edge

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

Definition at line 164 of file TetPolynomialBase.hpp.

◆ N_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::N_volume_bubble
protected

Definition at line 168 of file TetPolynomialBase.hpp.

◆ N_volume_edge

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

Definition at line 166 of file TetPolynomialBase.hpp.

◆ N_volume_face

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

Definition at line 167 of file TetPolynomialBase.hpp.

◆ senseFaceAlpha

MatrixInt MoFEM::TetPolynomialBase::senseFaceAlpha
protected

Definition at line 162 of file TetPolynomialBase.hpp.

◆ vPtr

const void* MoFEM::TetPolynomialBase::vPtr
protected

Definition at line 70 of file TetPolynomialBase.hpp.

The documentation for this struct was generated from the following files:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::Hcurl_Ainsworth_FaceFunctions_MBTET
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
TetBaseCache::hdivBrokenBaseInterior
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > hdivBrokenBaseInterior
Definition: TetPolynomialBase.cpp:70
g
constexpr double g
Definition: shallow_wave.cpp:63
MoFEM::TetPolynomialBase::getValueH1
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Get base functions for H1 space.
Definition: TetPolynomialBase.cpp:113
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET
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:401
NBEDGE_H1
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE
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:634
LASTBASE
@ LASTBASE
Definition: definitions.h:69
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
L2_Ainsworth_ShapeFunctions_MBTET
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on tetrahedron for L2 space.
Definition: l2.c:74
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv
static boost::function< int(int)> broken_nbvolumetet_volume_hdiv
Definition: Hdiv.hpp:29
MoFEM::TetPolynomialBase::getValueL2
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
Definition: TetPolynomialBase.cpp:525
NBVOLUMETET_DEMKOWICZ_HCURL
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:110
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
MoFEM::TetPolynomialBase::N_face_bubble
ublas::vector< MatrixDouble > N_face_bubble
Definition: TetPolynomialBase.hpp:165
MoFEM::TetPolynomialBase::getValueHcurlAinsworthBase
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1779
NOBASE
@ NOBASE
Definition: definitions.h:59
TetBaseCache::BaseCacheMI
boost::multi_index_container< BaseCacheItem, boost::multi_index::indexed_by< boost::multi_index::hashed_unique< composite_key< BaseCacheItem, member< BaseCacheItem, int, &BaseCacheItem::order >, member< BaseCacheItem, int, &BaseCacheItem::nb_gauss_pts > >> > > BaseCacheMI
Definition: TetPolynomialBase.cpp:47
MoFEM::Field::maxBrokenDofsOrder
static constexpr int maxBrokenDofsOrder
max number of broken dofs
Definition: FieldMultiIndices.hpp:282
ApproximationBaseNames
const static char *const ApproximationBaseNames[]
Definition: definitions.h:72
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
NBFACETRI_AINSWORTH_EDGE_HDIV
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:130
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3], double *diff_phi_f_e[3], int gdim, int nb, 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:47
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_face_hdiv
static boost::function< int(int)> broken_nbvolumetet_face_hdiv
Definition: Hdiv.hpp:28
MoFEM::Hcurl_Demkowicz_EdgeBaseFunctions_MBTET
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:2068
MoFEM::TetPolynomialBase::setDofsSideMapHdiv
static MoFEMErrorCode setDofsSideMapHdiv(const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &dofs_side_map)
Set the Dofs Side Map Hdiv object.
Definition: TetPolynomialBase.cpp:2159
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::TetPolynomialBase::senseFaceAlpha
MatrixInt senseFaceAlpha
Definition: TetPolynomialBase.hpp:162
MoFEM::EntPolynomialBaseCtx::dAta
EntitiesFieldData & dAta
Definition: EntPolynomialBaseCtx.hpp:36
NBVOLUMETET_AINSWORTH_FACE_HDIV
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:134
TetBaseCache::hdivBaseInterior
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > hdivBaseInterior
Definition: TetPolynomialBase.cpp:68
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
MoFEM::TetPolynomialBase::N_face_edge
ublas::matrix< MatrixDouble > N_face_edge
Definition: TetPolynomialBase.hpp:164
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
MoFEM::Hdiv_Demkowicz_Interior_MBTET
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:780
MoFEM::TetPolynomialBase::getValueH1AinsworthBase
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:131
MoFEM::TetPolynomialBase::diffN_volume_edge
ublas::vector< MatrixDouble > diffN_volume_edge
Definition: TetPolynomialBase.hpp:172
NBVOLUMETET_L2
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:27
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition: Templates.hpp:2013
TetBaseCache::l2BaseInterior
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > l2BaseInterior
Definition: TetPolynomialBase.cpp:74
MoFEM::TetPolynomialBase::getValueL2AinsworthBase
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:544
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:49
MoFEM::EntitiesFieldData::facesNodes
MatrixInt facesNodes
nodes on finite element faces
Definition: EntitiesFieldData.hpp:45
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
H1_EdgeShapeFunctions_MBTET
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:274
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::AinsworthOrderHooks::broken_nbfacetri_edge_hdiv
static boost::function< int(int)> broken_nbfacetri_edge_hdiv
Definition: Hdiv.hpp:25
MoFEM::EntPolynomialBaseCtx::copyNodeBase
const FieldApproximationBase copyNodeBase
Definition: EntPolynomialBaseCtx.hpp:41
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv
static boost::function< int(int)> broken_nbvolumetet_edge_hdiv
Definition: Hdiv.hpp:27
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::TetPolynomialBase::getValueHcurl
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
Definition: TetPolynomialBase.cpp:2029
MoFEM::BernsteinBezier::generateIndicesEdgeTet
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:138
TetBaseCache::HDivBaseCacheItem
Definition: TetPolynomialBase.cpp:19
MoFEM::BernsteinBezier::generateIndicesTriTet
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:174
NBVOLUMETET_AINSWORTH_VOLUME_HDIV
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:135
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::TetPolynomialBase::diffN_volume_bubble
MatrixDouble diffN_volume_bubble
Definition: TetPolynomialBase.hpp:174
MoFEM::getFTensor2HVecFromPtr< 3, 3 >
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:968
MoFEM::TetPolynomialBase::N_volume_edge
ublas::vector< MatrixDouble > N_volume_edge
Definition: TetPolynomialBase.hpp:166
NBFACETRI_DEMKOWICZ_HDIV
#define NBFACETRI_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:139
MoFEM::Types::MatrixInt
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:76
MoFEM::TetPolynomialBase::vPtr
const void * vPtr
Definition: TetPolynomialBase.hpp:70
MoFEM::TetPolynomialBase::getValueHdivAinsworthBaseImpl
MoFEMErrorCode getValueHdivAinsworthBaseImpl(MatrixDouble &pts, MatrixDouble &shape_functions, MatrixDouble &diff_shape_functions, int volume_order, std::array< int, 4 > &faces_order, std::array< int, 3 *4 > &faces_nodes, boost::function< int(int)> broken_nbfacetri_edge_hdiv, boost::function< int(int)> broken_nbfacetri_face_hdiv, boost::function< int(int)> broken_nbvolumetet_edge_hdiv, boost::function< int(int)> broken_nbvolumetet_face_hdiv, boost::function< int(int)> broken_nbvolumetet_volume_hdiv)
Definition: TetPolynomialBase.cpp:781
TetBaseCache::HDivBaseFaceCacheMI
boost::multi_index_container< HDivBaseCacheItem, boost::multi_index::indexed_by< boost::multi_index::hashed_unique< composite_key< HDivBaseCacheItem, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::order >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::nb_gauss_pts >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n0 >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n1 >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n2 > >> > > HDivBaseFaceCacheMI
Definition: TetPolynomialBase.cpp:65
MoFEM::EntPolynomialBaseCtx::sPace
const FieldSpace sPace
Definition: EntPolynomialBaseCtx.hpp:37
MoFEM::BernsteinBezier::generateIndicesVertexTet
static MoFEMErrorCode generateIndicesVertexTet(const int N, int *alpha)
Definition: BernsteinBezier.cpp:128
NBVOLUMETET_DEMKOWICZ_HDIV
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:140
MoFEM::Hcurl_Demkowicz_VolumeBaseFunctions_MBTET
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:2468
MoFEM::TetPolynomialBase::diffN_face_edge
ublas::matrix< MatrixDouble > diffN_face_edge
Definition: TetPolynomialBase.hpp:170
NBEDGE_DEMKOWICZ_HCURL
#define NBEDGE_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:108
TetBaseCache::BaseCacheItem
Definition: TetPolynomialBase.cpp:12
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::TetPolynomialBase::N_volume_face
ublas::vector< MatrixDouble > N_volume_face
Definition: TetPolynomialBase.hpp:167
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
TetBaseCache::hDivBaseFace
static std::array< std::map< const void *, HDivBaseFaceCacheMI >, LASTBASE > hDivBaseFace
Definition: TetPolynomialBase.cpp:72
MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb, 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:174
MoFEM::AinsworthOrderHooks::broken_nbfacetri_face_hdiv
static boost::function< int(int)> broken_nbfacetri_face_hdiv
Definition: Hdiv.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::TetPolynomialBase::getValueHdiv
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
Definition: TetPolynomialBase.cpp:1744
MoFEM::TetPolynomialBase::getValueHdivAinsworthBase
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:909
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:39
MoFEM::TetPolynomialBase::diffN_volume_face
ublas::vector< MatrixDouble > diffN_volume_face
Definition: TetPolynomialBase.hpp:173
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::TetPolynomialBase::getValueHcurlDemkowiczBase
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1897
NBFACETRI_DEMKOWICZ_HCURL
#define NBFACETRI_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:109
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEM::TetPolynomialBase::getValueL2BernsteinBezierBase
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:605
convert.n
n
Definition: convert.py:82
MoFEM::EntPolynomialBaseCtx::spaceContinuity
const FieldContinuity spaceContinuity
Definition: EntPolynomialBaseCtx.hpp:38
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::Tools::shapeFunMBTET
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:747
MoFEM::TetPolynomialBase::getValueHdivDemkowiczBrokenBase
MoFEMErrorCode getValueHdivDemkowiczBrokenBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1593
NBEDGE_AINSWORTH_HCURL
#define NBEDGE_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:97
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
MoFEM::BernsteinBezier::generateIndicesTetTet
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
Definition: BernsteinBezier.cpp:203
FieldSpaceNames
const static char *const FieldSpaceNames[]
Definition: definitions.h:92
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::TetPolynomialBase::cTx
EntPolynomialBaseCtx * cTx
Definition: TetPolynomialBase.hpp:71
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::Hcurl_Ainsworth_EdgeBaseFunctions_MBTET
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
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET
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:507
NBVOLUMETET_AINSWORTH_HCURL
#define NBVOLUMETET_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:105
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
DISCONTINUOUS
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
Definition: definitions.h:101
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MoFEM::TetPolynomialBase::diffN_face_bubble
ublas::vector< MatrixDouble > diffN_face_bubble
Definition: TetPolynomialBase.hpp:171
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MoFEM::TetPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1456
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::EntPolynomialBaseCtx::fieldName
const std::string fieldName
Definition: EntPolynomialBaseCtx.hpp:40
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET
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:307
MoFEM::Hcurl_Demkowicz_FaceBaseFunctions_MBTET
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:2395
MoFEM::TetPolynomialBase::N_volume_bubble
MatrixDouble N_volume_bubble
Definition: TetPolynomialBase.hpp:168
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:57
NBFACETRI_AINSWORTH_FACE_HDIV
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:131
tetCacheSwitch
auto tetCacheSwitch(const void *ptr, T &cache, std::string cache_name)
Definition: TetPolynomialBase.cpp:2271
H1_FaceShapeFunctions_MBTET
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:373
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
NBFACETRI_AINSWORTH_HCURL
#define NBFACETRI_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:100
MoFEM::BernsteinBezier::baseFunctionsTet
static MoFEMErrorCode baseFunctionsTet(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
Definition: BernsteinBezier.cpp:443
H1_VolumeShapeFunctions_MBTET
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:475
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
MoFEM::TetPolynomialBase::getValueH1BernsteinBezierBase
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:243
CONTINUOUS
@ CONTINUOUS
Regular field.
Definition: definitions.h:100
NBVOLUMETET_AINSWORTH_EDGE_HDIV
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:133
MoFEM::Hcurl_Ainsworth_VolumeFunctions_MBTET
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
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::TetPolynomialBase::getValueHdivAinsworthBrokenBase
MoFEMErrorCode getValueHdivAinsworthBrokenBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1170
MoFEM::TetPolynomialBase
Calculate base functions on tetrahedral.
Definition: TetPolynomialBase.hpp:17
MoFEM::EntPolynomialBaseCtx::basePolynomialsType0
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Definition: EntPolynomialBaseCtx.hpp:27
Generated by Doxygen 1.8.17 and hosted at University of Glasgow