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

Calculate base functions on tetrahedral. More...

#include <src/approximation/TetPolynomialBase.hpp>

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

Public Member Functions

MoFEMErrorCode query_interface (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...
 

Private 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 Private Member Functions

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

Private 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 2021 of file TetPolynomialBase.cpp.

2022  {
2023  MoFEMFunctionBegin;
2024 
2025  cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
2026 
2027  int nb_gauss_pts = pts.size2();
2028  if (!nb_gauss_pts)
2029  MoFEMFunctionReturnHot(0);
2030 
2031  if (pts.size1() < 3)
2032  SETERRQ(
2033  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2034  "Wrong dimension of pts, should be at least 3 rows with coordinates");
2035 
2036  if (cTx->bAse != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
2037  const FieldApproximationBase base = cTx->bAse;
2038  EntitiesFieldData &data = cTx->dAta;
2039  if (cTx->copyNodeBase == LASTBASE) {
2040  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
2041  false);
2042  CHKERR Tools::shapeFunMBTET(
2043  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
2044  &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
2045  } else {
2046  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
2047  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
2048  }
2049  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
2050  (unsigned int)nb_gauss_pts) {
2051  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2052  "Base functions or nodes has wrong number of integration points "
2053  "for base %s",
2054  ApproximationBaseNames[base]);
2055  }
2056  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
2057  std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
2058  data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
2059  }
2060 
2061  switch (cTx->spaceContinuity) {
2062  case CONTINUOUS:
2063 
2064  switch (cTx->sPace) {
2065  case H1:
2066  CHKERR getValueH1(pts);
2067  break;
2068  case HDIV:
2069  CHKERR getValueHdiv(pts);
2070  break;
2071  case HCURL:
2072  CHKERR getValueHcurl(pts);
2073  break;
2074  case L2:
2075  CHKERR getValueL2(pts);
2076  break;
2077  default:
2078  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2079  FieldSpaceNames[cTx->sPace]);
2080  }
2081  break;
2082 
2083  case DISCONTINUOUS:
2084 
2085  switch (cTx->sPace) {
2086  case HDIV:
2087  CHKERR getValueHdiv(pts);
2088  break;
2089  default:
2090  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2091  FieldSpaceNames[cTx->sPace]);
2092  }
2093  break;
2094 
2095  default:
2096  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
2097  }
2098 
2099  MoFEMFunctionReturn(0);
2100 }

◆ getValueH1()

MoFEMErrorCode TetPolynomialBase::getValueH1 ( MatrixDouble &  pts)
private

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)
private

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)
private

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)
private

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 2002 of file TetPolynomialBase.cpp.

2002  {
2003  MoFEMFunctionBegin;
2004 
2005  switch (cTx->bAse) {
2006  case AINSWORTH_LEGENDRE_BASE:
2007  case AINSWORTH_LOBATTO_BASE:
2008  CHKERR getValueHcurlAinsworthBase(pts);
2009  break;
2010  case DEMKOWICZ_JACOBI_BASE:
2011  CHKERR getValueHcurlDemkowiczBase(pts);
2012  break;
2013  default:
2014  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
2015  }
2016 
2017  MoFEMFunctionReturn(0);
2018 }

◆ getValueHcurlAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlAinsworthBase ( MatrixDouble &  pts)
private

Definition at line 1752 of file TetPolynomialBase.cpp.

1752  {
1753  MoFEMFunctionBegin;
1754 
1755  EntitiesFieldData &data = cTx->dAta;
1756  const FieldApproximationBase base = cTx->bAse;
1757  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
1758  double *diffL, const int dim) =
1759  cTx->basePolynomialsType0;
1760 
1761  int nb_gauss_pts = pts.size2();
1762 
1763  // edges
1764  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1765  int sense[6], order[6];
1766  if (data.dataOnEntities[MBEDGE].size() != 6) {
1767  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1768  }
1769  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1770  for (int ee = 0; ee != 6; ee++) {
1771  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1772  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1773  "data inconsistency");
1774  }
1775  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1776  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1777  int nb_dofs =
1778  NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1779  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1780  3 * nb_dofs, false);
1781  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1782  9 * nb_dofs, false);
1783  hcurl_edge_n[ee] =
1784  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1785  diff_hcurl_edge_n[ee] =
1786  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1787  }
1788  CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(
1789  sense, order,
1790  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1791  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1792  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
1793  } else {
1794  for (int ee = 0; ee != 6; ee++) {
1795  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1796  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1797  false);
1798  }
1799  }
1800 
1801  // triangles
1802  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1803  int order[4];
1804  // faces
1805  if (data.dataOnEntities[MBTRI].size() != 4) {
1806  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1807  }
1808  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1809  for (int ff = 0; ff != 4; ff++) {
1810  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1811  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1812  "data inconsistency");
1813  }
1814  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1815  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
1816  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1817  3 * nb_dofs, false);
1818  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1819  9 * nb_dofs, false);
1820  hcurl_base_n[ff] =
1821  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1822  diff_hcurl_base_n[ff] =
1823  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1824  }
1825  if (data.facesNodes.size1() != 4) {
1826  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1827  }
1828  if (data.facesNodes.size2() != 3) {
1829  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1830  }
1831  CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET(
1832  &*data.facesNodes.data().begin(), order,
1833  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1834  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1835  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
1836  } else {
1837  for (int ff = 0; ff != 4; ff++) {
1838  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1839  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1840  false);
1841  }
1842  }
1843 
1844  if (data.spacesOnEntities[MBTET].test(HCURL)) {
1845 
1846  // volume
1847  int order = data.dataOnEntities[MBTET][0].getOrder();
1848  int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
1849  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1850  3 * nb_vol_dofs, false);
1851  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1852  9 * nb_vol_dofs, false);
1853  CHKERR Hcurl_Ainsworth_VolumeFunctions_MBTET(
1854  data.dataOnEntities[MBTET][0].getOrder(),
1855  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1856  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1857  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1858  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1859  nb_gauss_pts, base_polynomials);
1860 
1861  } else {
1862  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1863  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1864  }
1865 
1866  MoFEMFunctionReturn(0);
1867 }

◆ getValueHcurlDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlDemkowiczBase ( MatrixDouble &  pts)
private

Definition at line 1870 of file TetPolynomialBase.cpp.

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

◆ getValueHdiv()

MoFEMErrorCode TetPolynomialBase::getValueHdiv ( MatrixDouble &  pts)
private

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 1717 of file TetPolynomialBase.cpp.

1717  {
1718  MoFEMFunctionBeginHot;
1719 
1720  switch (cTx->spaceContinuity) {
1721  case CONTINUOUS:
1722  switch (cTx->bAse) {
1723  case AINSWORTH_LEGENDRE_BASE:
1724  case AINSWORTH_LOBATTO_BASE:
1725  return getValueHdivAinsworthBase(pts);
1726  case DEMKOWICZ_JACOBI_BASE:
1727  return getValueHdivDemkowiczBase(pts);
1728  default:
1729  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1730  }
1731  break;
1732  case DISCONTINUOUS:
1733  switch (cTx->bAse) {
1734  case AINSWORTH_LEGENDRE_BASE:
1735  case AINSWORTH_LOBATTO_BASE:
1736  return getValueHdivAinsworthBrokenBase(pts);
1737  case DEMKOWICZ_JACOBI_BASE:
1738  return getValueHdivDemkowiczBrokenBase(pts);
1739  default:
1740  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1741  }
1742  break;
1743 
1744  default:
1745  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1746  }
1747 
1748  MoFEMFunctionReturnHot(0);
1749 }

◆ getValueHdivAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBase ( MatrixDouble &  pts)
private

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  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1012  for (int oo = 0; oo != face_order; oo++) {
1013 
1014  // face-edge
1015  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(
1016  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(oo));
1017  dd !=
1018  NBFACETRI_AINSWORTH_EDGE_HDIV(
1019  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(oo + 1));
1020  dd++) {
1021  for (int ee = 0; ee != 3; ++ee) {
1022  t_base(i) = t_base_f_e[ff * 3 + ee](i);
1023  ++t_base;
1024  ++t_base_f_e[ff * 3 + ee];
1025  }
1026  for (int ee = 0; ee != 3; ++ee) {
1027  t_diff_base(i, j) = t_diff_base_f_e[ff * 3 + ee](i, j);
1028  ++t_diff_base;
1029  ++t_diff_base_f_e[ff * 3 + ee];
1030  }
1031  }
1032 
1033  // face-face
1034  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(
1035  AinsworthOrderHooks::broken_nbfacetri_face_hdiv(oo));
1036  dd !=
1037  NBFACETRI_AINSWORTH_FACE_HDIV(
1038  AinsworthOrderHooks::broken_nbfacetri_face_hdiv(oo + 1));
1039  dd++) {
1040  t_base(i) = t_base_f_f[ff](i);
1041  ++t_base;
1042  ++t_base_f_f[ff];
1043  t_diff_base(i, j) = t_diff_base_f_f[ff](i, j);
1044  ++t_diff_base;
1045  ++t_diff_base_f_f[ff];
1046  }
1047  }
1048  }
1049  }
1050  }
1051 
1052  // volume
1053  int volume_dofs =
1054  6 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(
1055  AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order)) +
1056  4 * NBVOLUMETET_AINSWORTH_FACE_HDIV(
1057  AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order)) +
1058  NBVOLUMETET_AINSWORTH_VOLUME_HDIV(
1059  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order));
1060  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * volume_dofs,
1061  false);
1062  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1063  9 * volume_dofs, false);
1064  if (volume_dofs) {
1065  double *base_ptr =
1066  &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1067  double *diff_base_ptr =
1068  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1069  auto t_base = getFTensor1FromPtr<3>(base_ptr);
1070  auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1071 
1072  // volume-edge
1073  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_e[] = {
1074  getFTensor1FromPtr<3>(&*N_volume_edge[0].data().begin()),
1075  getFTensor1FromPtr<3>(&*N_volume_edge[1].data().begin()),
1076  getFTensor1FromPtr<3>(&*N_volume_edge[2].data().begin()),
1077  getFTensor1FromPtr<3>(&*N_volume_edge[3].data().begin()),
1078  getFTensor1FromPtr<3>(&*N_volume_edge[4].data().begin()),
1079  getFTensor1FromPtr<3>(&*N_volume_edge[5].data().begin())};
1080  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_e[] = {
1081  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[0].data().begin()),
1082  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[1].data().begin()),
1083  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[2].data().begin()),
1084  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[3].data().begin()),
1085  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[4].data().begin()),
1086  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[5].data().begin())};
1087 
1088  // volume-faces
1089  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_f[] = {
1090  getFTensor1FromPtr<3>(&*N_volume_face[0].data().begin()),
1091  getFTensor1FromPtr<3>(&*N_volume_face[1].data().begin()),
1092  getFTensor1FromPtr<3>(&*N_volume_face[2].data().begin()),
1093  getFTensor1FromPtr<3>(&*N_volume_face[3].data().begin())};
1094  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1095  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[0].data().begin()),
1096  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[1].data().begin()),
1097  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[2].data().begin()),
1098  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[3].data().begin())};
1099 
1100  // volume-bubble
1101  base_ptr = &*(N_volume_bubble.data().begin());
1102  diff_base_ptr = &*(diffN_volume_bubble.data().begin());
1103  auto t_base_v = getFTensor1FromPtr<3>(base_ptr);
1104  auto t_diff_base_v = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1105 
1106  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1107  for (int oo = 0; oo < volume_order; oo++) {
1108 
1109  // volume-edge
1110  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(
1111  AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(oo));
1112  dd !=
1113  NBVOLUMETET_AINSWORTH_EDGE_HDIV(
1114  AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(oo + 1));
1115  dd++) {
1116  for (int ee = 0; ee < 6; ee++) {
1117  t_base(i) = t_base_v_e[ee](i);
1118  ++t_base;
1119  ++t_base_v_e[ee];
1120  t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
1121  ++t_diff_base;
1122  ++t_diff_base_v_e[ee];
1123  }
1124  }
1125 
1126  // volume-face
1127  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(
1128  AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(oo));
1129  dd <
1130  NBVOLUMETET_AINSWORTH_FACE_HDIV(
1131  AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(oo + 1));
1132  dd++) {
1133  for (int ff = 0; ff < 4; ff++) {
1134  t_base(i) = t_base_v_f[ff](i);
1135  ++t_base;
1136  ++t_base_v_f[ff];
1137  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1138  ++t_diff_base;
1139  ++t_diff_base_v_f[ff];
1140  }
1141  }
1142 
1143  // volume-bubble
1144  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(
1145  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(oo));
1146  dd <
1147  NBVOLUMETET_AINSWORTH_VOLUME_HDIV(
1148  AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(oo + 1));
1149  dd++) {
1150  t_base(i) = t_base_v(i);
1151  ++t_base;
1152  ++t_base_v;
1153  t_diff_base(i, j) = t_diff_base_v(i, j);
1154  ++t_diff_base;
1155  ++t_diff_base_v;
1156  }
1157  }
1158  }
1159  }
1160 
1161  MoFEMFunctionReturn(0);
1162 }

◆ 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 
)
private

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)
private

Definition at line 1165 of file TetPolynomialBase.cpp.

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

◆ getValueHdivDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBase ( MatrixDouble &  pts)
private

Definition at line 1429 of file TetPolynomialBase.cpp.

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

◆ getValueHdivDemkowiczBrokenBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBrokenBase ( MatrixDouble &  pts)
private

Definition at line 1566 of file TetPolynomialBase.cpp.

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

◆ getValueL2()

MoFEMErrorCode TetPolynomialBase::getValueL2 ( MatrixDouble &  pts)
private

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)
private

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  MoFEMFunctionBeginHot(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)
private

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 2102 of file TetPolynomialBase.cpp.

2107  {
2108  MoFEMFunctionBegin;
2109 
2110  switch (continuity) {
2111  case DISCONTINUOUS:
2112 
2113  switch (space) {
2114  case HDIV:
2115  CHKERR setDofsSideMapHdiv(continuity, base, dofs_side_map);
2116  break;
2117  default:
2118  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2119  FieldSpaceNames[space]);
2120  }
2121  break;
2122 
2123  default:
2124  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2125  "Unknown (or not implemented) continuity");
2126  }
2127 
2128  MoFEMFunctionReturn(0);
2129 }

◆ setDofsSideMapHdiv()

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

Set the Dofs Side Map Hdiv object.

Parameters
space
continuity
base
dofs_side_map
Returns
MoFEMErrorCode

Definition at line 2132 of file TetPolynomialBase.cpp.

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

◆ 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 2277 of file TetPolynomialBase.cpp.

2278  {
2279  return tetCacheSwitch(ptr, TetBaseCache::hDivBaseFace[base],
2280  std::string("hDivBaseFace") +
2281  ApproximationBaseNames[base]);
2282 }

◆ 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 2285 of file TetPolynomialBase.cpp.

2286  {
2287  return tetCacheSwitch(ptr, TetBaseCache::hdivBaseInterior[base],
2288  std::string("hdivBaseInterior") +
2289  ApproximationBaseNames[base]);
2290 }

◆ switchCacheBaseInterior() [5/5]

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

Definition at line 2347 of file TetPolynomialBase.cpp.

2348  {
2349  return tetCacheSwitch(ptr, TetBaseCache::l2BaseInterior[base],
2350  std::string("hdivBaseInterior") +
2351  ApproximationBaseNames[base]);
2352 }

◆ 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 2317 of file TetPolynomialBase.cpp.

2318  {
2319  for (auto fe_ptr : v) {
2320  if (TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr)) {
2321  TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr);
2322  }
2323  if (TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr)) {
2324  TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr);
2325  }
2326  if (TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr)) {
2327  TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr);
2328  }
2329  }
2330 }

◆ switchCacheBaseOff() [5/10]

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

Definition at line 2365 of file TetPolynomialBase.cpp.

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

◆ 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 2340 of file TetPolynomialBase.cpp.

2340  {
2341  for (auto b = 0; b != LASTBASE; ++b) {
2342  switchCacheBaseOff<HDIV>(static_cast<FieldApproximationBase>(b), v);
2343  }
2344 }

◆ switchCacheBaseOff() [10/10]

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

Definition at line 2382 of file TetPolynomialBase.cpp.

2382  {
2383  for (auto b = 0; b != LASTBASE; ++b) {
2384  switchCacheBaseOff<L2>(static_cast<FieldApproximationBase>(b), v);
2385  }
2386 }

◆ 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 2301 of file TetPolynomialBase.cpp.

2302  {
2303  for (auto fe_ptr : v) {
2304  if (!TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr)) {
2305  TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr);
2306  }
2307  if (!TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr)) {
2308  TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr);
2309  }
2310  if (!TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr)) {
2311  TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr);
2312  }
2313  }
2314 }

◆ switchCacheBaseOn() [5/10]

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

Definition at line 2355 of file TetPolynomialBase.cpp.

2356  {
2357  for (auto fe_ptr : v) {
2358  if (!TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr)) {
2359  TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr);
2360  }
2361  }
2362 }

◆ 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 2333 of file TetPolynomialBase.cpp.

2333  {
2334  for (auto b = 0; b != LASTBASE; ++b) {
2335  switchCacheBaseOn<HDIV>(static_cast<FieldApproximationBase>(b), v);
2336  }
2337 }

◆ switchCacheBaseOn() [10/10]

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

Definition at line 2375 of file TetPolynomialBase.cpp.

2375  {
2376  for (auto b = 0; b != LASTBASE; ++b) {
2377  switchCacheBaseOn<L2>(static_cast<FieldApproximationBase>(b), v);
2378  }
2379 }

◆ 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 2293 of file TetPolynomialBase.cpp.

2294  {
2295  return tetCacheSwitch(ptr, TetBaseCache::hdivBrokenBaseInterior[base],
2296  std::string("hdivBrokenBaseInterior") +
2297  ApproximationBaseNames[base]);
2298 }

Member Data Documentation

◆ cTx

EntPolynomialBaseCtx* MoFEM::TetPolynomialBase::cTx
private

Definition at line 71 of file TetPolynomialBase.hpp.

◆ diffN_face_bubble

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

Definition at line 171 of file TetPolynomialBase.hpp.

◆ diffN_face_edge

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

Definition at line 170 of file TetPolynomialBase.hpp.

◆ diffN_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::diffN_volume_bubble
private

Definition at line 174 of file TetPolynomialBase.hpp.

◆ diffN_volume_edge

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

Definition at line 172 of file TetPolynomialBase.hpp.

◆ diffN_volume_face

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

Definition at line 173 of file TetPolynomialBase.hpp.

◆ N_face_bubble

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

Definition at line 165 of file TetPolynomialBase.hpp.

◆ N_face_edge

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

Definition at line 164 of file TetPolynomialBase.hpp.

◆ N_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::N_volume_bubble
private

Definition at line 168 of file TetPolynomialBase.hpp.

◆ N_volume_edge

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

Definition at line 166 of file TetPolynomialBase.hpp.

◆ N_volume_face

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

Definition at line 167 of file TetPolynomialBase.hpp.

◆ senseFaceAlpha

MatrixInt MoFEM::TetPolynomialBase::senseFaceAlpha
private

Definition at line 162 of file TetPolynomialBase.hpp.

◆ vPtr

const void* MoFEM::TetPolynomialBase::vPtr
private

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
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:1752
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:2132
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:2011
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:2002
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:1717
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:1870
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:738
MoFEM::TetPolynomialBase::getValueHdivDemkowiczBrokenBase
MoFEMErrorCode getValueHdivDemkowiczBrokenBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1566
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::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:1429
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:2259
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:1165
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