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

Calculate base functions on tetrahedral. More...

#include <src/approximation/TetPolynomialBase.hpp>

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

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 
 TetPolynomialBase ()=default
 
 ~TetPolynomialBase ()=default
 
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 refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

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 getValueHdivAinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurlAinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivDemkowiczBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurlDemkowiczBase (MatrixDouble &pts)
 

Private Attributes

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

Additional Inherited Members

- 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...
 

Detailed Description

Calculate base functions on tetrahedral.

Definition at line 31 of file TetPolynomialBase.hpp.

Constructor & Destructor Documentation

◆ TetPolynomialBase()

MoFEM::TetPolynomialBase::TetPolynomialBase ( )
default

◆ ~TetPolynomialBase()

MoFEM::TetPolynomialBase::~TetPolynomialBase ( )
default

Member Function Documentation

◆ getValue()

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

Reimplemented from MoFEM::BaseFunction.

Definition at line 1362 of file TetPolynomialBase.cpp.

1363  {
1365 
1366  cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
1367 
1368  int nb_gauss_pts = pts.size2();
1369  if (!nb_gauss_pts)
1371 
1372  if (pts.size1() < 3)
1373  SETERRQ(
1374  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1375  "Wrong dimension of pts, should be at least 3 rows with coordinates");
1376 
1378  const FieldApproximationBase base = cTx->bAse;
1379  EntitiesFieldData &data = cTx->dAta;
1380  if (cTx->copyNodeBase == LASTBASE) {
1381  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
1382  false);
1384  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1385  &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
1386  } else {
1387  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
1388  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
1389  }
1390  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
1391  (unsigned int)nb_gauss_pts) {
1392  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1393  "Base functions or nodes has wrong number of integration points "
1394  "for base %s",
1395  ApproximationBaseNames[base]);
1396  }
1397  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
1398  std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
1399  data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1400  }
1401 
1402  switch (cTx->sPace) {
1403  case H1:
1404  CHKERR getValueH1(pts);
1405  break;
1406  case HDIV:
1407  CHKERR getValueHdiv(pts);
1408  break;
1409  case HCURL:
1410  CHKERR getValueHcurl(pts);
1411  break;
1412  case L2:
1413  CHKERR getValueL2(pts);
1414  break;
1415  default:
1416  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space");
1417  }
1418 
1420 }
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ LASTBASE
Definition: definitions.h:82
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
@ HCURL
field with continuous tangents
Definition: definitions.h:99
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
static const char *const ApproximationBaseNames[]
Definition: definitions.h:85
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase copyNodeBase
const FieldApproximationBase bAse
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
EntPolynomialBaseCtx * cTx
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Get base functions for H1 space.
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
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:607
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:250
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ getValueH1()

MoFEMErrorCode TetPolynomialBase::getValueH1 ( MatrixDouble pts)
private

Get base functions for H1 space.

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

Definition at line 33 of file TetPolynomialBase.cpp.

33  {
35 
36  switch (cTx->bAse) {
40  break;
43  break;
44  default:
45  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
46  }
47 
49 }
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:75
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)

◆ getValueH1AinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueH1AinsworthBase ( MatrixDouble pts)
private

Definition at line 51 of file TetPolynomialBase.cpp.

51  {
53 
54  EntitiesFieldData &data = cTx->dAta;
55  const FieldApproximationBase base = cTx->bAse;
56  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
57  double *diffL, const int dim) =
59 
60  int nb_gauss_pts = pts.size2();
61 
62  int sense[6], order[6];
63  if (data.spacesOnEntities[MBEDGE].test(H1)) {
64  // edges
65  if (data.dataOnEntities[MBEDGE].size() != 6) {
66  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
67  }
68  double *h1_edge_n[6], *diff_h1_egde_n[6];
69  for (int ee = 0; ee != 6; ++ee) {
70  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
71  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
72  "data inconsistency");
73  }
74  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
75  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
76  int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
77  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
78  false);
79  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
80  3 * nb_dofs, false);
81  h1_edge_n[ee] =
82  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
83  diff_h1_egde_n[ee] =
84  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
85  }
87  sense, order,
88  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
89  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
90  h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
91  } else {
92  for (int ee = 0; ee != 6; ++ee) {
93  data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
94  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
95  }
96  }
97 
98  if (data.spacesOnEntities[MBTRI].test(H1)) {
99  // faces
100  if (data.dataOnEntities[MBTRI].size() != 4) {
101  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
102  }
103  double *h1_face_n[4], *diff_h1_face_n[4];
104  for (int ff = 0; ff != 4; ++ff) {
105  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
106  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
107  "data inconsistency");
108  }
109  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getOrder());
110  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
111  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
112  false);
113  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
114  3 * nb_dofs, false);
115  h1_face_n[ff] =
116  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
117  diff_h1_face_n[ff] =
118  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
119  }
120  if (data.facesNodes.size1() != 4) {
121  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
122  }
123  if (data.facesNodes.size2() != 3) {
124  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
125  }
127  &*data.facesNodes.data().begin(), order,
128  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
129  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
130  h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
131 
132  } else {
133  for (int ff = 0; ff != 4; ++ff) {
134  data.dataOnEntities[MBTRI][ff].getN(base).resize(0, false);
135  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, false);
136  }
137  }
138 
139  if (data.spacesOnEntities[MBTET].test(H1)) {
140  // volume
141  int order = data.dataOnEntities[MBTET][0].getOrder();
142  int nb_vol_dofs = NBVOLUMETET_H1(order);
143  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
144  false);
145  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
146  3 * nb_vol_dofs, false);
148  data.dataOnEntities[MBTET][0].getOrder(),
149  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
150  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
151  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
152  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
153  nb_gauss_pts, base_polynomials);
154  } else {
155  data.dataOnEntities[MBTET][0].getN(base).resize(0, 0, false);
156  data.dataOnEntities[MBTET][0].getDiffN(base).resize(0, 0, false);
157  }
158 
160 }
static Index< 'p', 3 > p
const int dim
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
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:284
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
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:485
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:383
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types

◆ getValueH1BernsteinBezierBase()

MoFEMErrorCode TetPolynomialBase::getValueH1BernsteinBezierBase ( MatrixDouble pts)
private

Definition at line 163 of file TetPolynomialBase.cpp.

163  {
165 
166  EntitiesFieldData &data = cTx->dAta;
167  const std::string field_name = cTx->fieldName;
168  const int nb_gauss_pts = pts.size2();
169 
170  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
171  (unsigned int)nb_gauss_pts)
172  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
173  "Base functions or nodes has wrong number of integration points "
174  "for base %s",
176  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
177 
178  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
179  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
180  if (!ptr)
181  ptr.reset(new MatrixInt());
182  return *ptr;
183  };
184 
185  auto get_base = [field_name](auto &data) -> MatrixDouble & {
186  auto &ptr = data.getBBNSharedPtr(field_name);
187  if (!ptr)
188  ptr.reset(new MatrixDouble());
189  return *ptr;
190  };
191 
192  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
193  auto &ptr = data.getBBDiffNSharedPtr(field_name);
194  if (!ptr)
195  ptr.reset(new MatrixDouble());
196  return *ptr;
197  };
198 
199  auto get_alpha_by_name_ptr =
200  [](auto &data,
201  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
202  return data.getBBAlphaIndicesSharedPtr(field_name);
203  };
204 
205  auto get_base_by_name_ptr =
206  [](auto &data,
207  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
208  return data.getBBNSharedPtr(field_name);
209  };
210 
211  auto get_diff_base_by_name_ptr =
212  [](auto &data,
213  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
214  return data.getBBDiffNSharedPtr(field_name);
215  };
216 
217  auto get_alpha_by_order_ptr =
218  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
219  return data.getBBAlphaIndicesByOrderSharedPtr(o);
220  };
221 
222  auto get_base_by_order_ptr =
223  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
224  return data.getBBNByOrderSharedPtr(o);
225  };
226 
227  auto get_diff_base_by_order_ptr =
228  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
229  return data.getBBDiffNByOrderSharedPtr(o);
230  };
231 
232  auto &vert_ent_data = data.dataOnEntities[MBVERTEX][0];
233  auto &vertex_alpha = get_alpha(vert_ent_data);
234  vertex_alpha.resize(4, 4, false);
235  vertex_alpha.clear();
236  for (int n = 0; n != 4; ++n)
237  vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
238 
239  auto &vert_get_n = get_base(vert_ent_data);
240  auto &vert_get_diff_n = get_diff_base(vert_ent_data);
241  vert_get_n.resize(nb_gauss_pts, 4, false);
242  vert_get_diff_n.resize(nb_gauss_pts, 12, false);
244  1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
245  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &vert_get_n(0, 0),
246  &vert_get_diff_n(0, 0));
247  for (int n = 0; n != 4; ++n) {
248  const double f = boost::math::factorial<double>(
249  data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
250  for (int g = 0; g != nb_gauss_pts; ++g) {
251  vert_get_n(g, n) *= f;
252  for (int d = 0; d != 3; ++d)
253  vert_get_diff_n(g, 3 * n + d) *= f;
254  }
255  }
256 
257  // edges
258  if (data.spacesOnEntities[MBEDGE].test(H1)) {
259  if (data.dataOnEntities[MBEDGE].size() != 6)
260  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
261  "Wrong size of ent data");
262 
263  constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
264  {0, 3}, {1, 3}, {2, 3}};
265  for (int ee = 0; ee != 6; ++ee) {
266  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
267  const int sense = ent_data.getSense();
268  if (sense == 0)
269  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
270  "Sense of the edge unknown");
271  const int order = ent_data.getOrder();
272  const int nb_dofs = NBEDGE_H1(order);
273 
274  if (nb_dofs) {
275  if (get_alpha_by_order_ptr(ent_data, order)) {
276  get_alpha_by_name_ptr(ent_data, field_name) =
277  get_alpha_by_order_ptr(ent_data, order);
278  get_base_by_name_ptr(ent_data, field_name) =
279  get_base_by_order_ptr(ent_data, order);
280  get_diff_base_by_name_ptr(ent_data, field_name) =
281  get_diff_base_by_order_ptr(ent_data, order);
282  } else {
283  auto &get_n = get_base(ent_data);
284  auto &get_diff_n = get_diff_base(ent_data);
285  get_n.resize(nb_gauss_pts, nb_dofs, false);
286  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
287 
288  auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
289  edge_alpha.resize(nb_dofs, 4, false);
291  &edge_alpha(0, 0));
292  if (sense == -1) {
293  for (int i = 0; i != edge_alpha.size1(); ++i) {
294  int a = edge_alpha(i, edges_nodes[ee][0]);
295  edge_alpha(i, edges_nodes[ee][0]) =
296  edge_alpha(i, edges_nodes[ee][1]);
297  edge_alpha(i, edges_nodes[ee][1]) = a;
298  }
299  }
301  order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
302  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
303  &get_diff_n(0, 0));
304 
305  get_alpha_by_order_ptr(ent_data, order) =
306  get_alpha_by_name_ptr(ent_data, field_name);
307  get_base_by_order_ptr(ent_data, order) =
308  get_base_by_name_ptr(ent_data, field_name);
309  get_diff_base_by_order_ptr(ent_data, order) =
310  get_diff_base_by_name_ptr(ent_data, field_name);
311  }
312  }
313  }
314  } else {
315  for (int ee = 0; ee != 6; ++ee) {
316  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
317  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
318  auto &get_n = get_base(ent_data);
319  auto &get_diff_n = get_diff_base(ent_data);
320  get_n.resize(nb_gauss_pts, 0, false);
321  get_diff_n.resize(nb_gauss_pts, 0, false);
322  }
323  }
324 
325  // face
326  if (data.spacesOnEntities[MBTRI].test(H1)) {
327  if (data.dataOnEntities[MBTRI].size() != 4)
328  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
329  "Wrong size of ent data");
330  if (data.facesNodes.size1() != 4)
331  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
332  if (data.facesNodes.size2() != 3)
333  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
334 
335  for (int ff = 0; ff != 4; ++ff) {
336  auto &ent_data = data.dataOnEntities[MBTRI][ff];
337  const int order = ent_data.getOrder();
338  const int nb_dofs = NBFACETRI_H1(order);
339 
340  if (nb_dofs) {
341  if (get_alpha_by_order_ptr(ent_data, order)) {
342  get_alpha_by_name_ptr(ent_data, field_name) =
343  get_alpha_by_order_ptr(ent_data, order);
344  get_base_by_name_ptr(ent_data, field_name) =
345  get_base_by_order_ptr(ent_data, order);
346  get_diff_base_by_name_ptr(ent_data, field_name) =
347  get_diff_base_by_order_ptr(ent_data, order);
348  } else {
349 
350  auto &get_n = get_base(ent_data);
351  auto &get_diff_n = get_diff_base(ent_data);
352  get_n.resize(nb_gauss_pts, nb_dofs, false);
353  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
354 
355  auto &face_alpha = get_alpha(ent_data);
356  face_alpha.resize(nb_dofs, 4, false);
357 
359  &face_alpha(0, 0));
360  senseFaceAlpha.resize(face_alpha.size1(), face_alpha.size2(), false);
361  senseFaceAlpha.clear();
362  constexpr int tri_nodes[4][3] = {
363  {0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
364  for (int d = 0; d != nb_dofs; ++d)
365  for (int n = 0; n != 3; ++n)
366  senseFaceAlpha(d, data.facesNodes(ff, n)) =
367  face_alpha(d, tri_nodes[ff][n]);
368  face_alpha.swap(senseFaceAlpha);
370  order, lambda.size1(), face_alpha.size1(), &face_alpha(0, 0),
371  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
372  &get_diff_n(0, 0));
373 
374  get_alpha_by_order_ptr(ent_data, order) =
375  get_alpha_by_name_ptr(ent_data, field_name);
376  get_base_by_order_ptr(ent_data, order) =
377  get_base_by_name_ptr(ent_data, field_name);
378  get_diff_base_by_order_ptr(ent_data, order) =
379  get_diff_base_by_name_ptr(ent_data, field_name);
380  }
381  }
382  }
383  } else {
384  for (int ff = 0; ff != 4; ++ff) {
385  auto &ent_data = data.dataOnEntities[MBTRI][ff];
386  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
387  auto &get_n = get_base(ent_data);
388  auto &get_diff_n = get_diff_base(ent_data);
389  get_n.resize(nb_gauss_pts, 0, false);
390  get_diff_n.resize(nb_gauss_pts, 0, false);
391  }
392  }
393 
394  if (data.spacesOnEntities[MBTET].test(H1)) {
395  if (data.dataOnEntities[MBTET].size() != 1)
396  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
397  "Wrong size ent of ent data");
398 
399  auto &ent_data = data.dataOnEntities[MBTET][0];
400  const int order = ent_data.getOrder();
401  const int nb_dofs = NBVOLUMETET_H1(order);
402  if (get_alpha_by_order_ptr(ent_data, order)) {
403  get_alpha_by_name_ptr(ent_data, field_name) =
404  get_alpha_by_order_ptr(ent_data, order);
405  get_base_by_name_ptr(ent_data, field_name) =
406  get_base_by_order_ptr(ent_data, order);
407  get_diff_base_by_name_ptr(ent_data, field_name) =
408  get_diff_base_by_order_ptr(ent_data, order);
409  } else {
410 
411  auto &get_n = get_base(ent_data);
412  auto &get_diff_n = get_diff_base(ent_data);
413  get_n.resize(nb_gauss_pts, nb_dofs, false);
414  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
415  if (nb_dofs) {
416  auto &tet_alpha = get_alpha(ent_data);
417  tet_alpha.resize(nb_dofs, 4, false);
418 
421  order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
422  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
423  &get_diff_n(0, 0));
424 
425  get_alpha_by_order_ptr(ent_data, order) =
426  get_alpha_by_name_ptr(ent_data, field_name);
427  get_base_by_order_ptr(ent_data, order) =
428  get_base_by_name_ptr(ent_data, field_name);
429  get_diff_base_by_order_ptr(ent_data, order) =
430  get_diff_base_by_name_ptr(ent_data, field_name);
431  }
432  }
433  } else {
434  auto &ent_data = data.dataOnEntities[MBTET][0];
435  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
436  auto &get_n = get_base(ent_data);
437  auto &get_diff_n = get_diff_base(ent_data);
438  get_n.resize(nb_gauss_pts, 0, false);
439  get_diff_n.resize(nb_gauss_pts, 0, false);
440  }
441 
443 }
static Index< 'n', 3 > n
static Index< 'o', 3 > o
constexpr double a
@ NOBASE
Definition: definitions.h:72
constexpr double lambda
FTensor::Index< 'i', SPACE_DIM > i
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:87
constexpr double g
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
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)
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)

◆ getValueHcurl()

MoFEMErrorCode TetPolynomialBase::getValueHcurl ( MatrixDouble pts)
private

Get base functions for Hcurl space.

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

Definition at line 1343 of file TetPolynomialBase.cpp.

1343  {
1345 
1346  switch (cTx->bAse) {
1350  break;
1351  case DEMKOWICZ_JACOBI_BASE:
1353  break;
1354  default:
1355  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1356  }
1357 
1359 }
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)

◆ getValueHcurlAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlAinsworthBase ( MatrixDouble pts)
private

Definition at line 1093 of file TetPolynomialBase.cpp.

1093  {
1095 
1096  EntitiesFieldData &data = cTx->dAta;
1097  const FieldApproximationBase base = cTx->bAse;
1098  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
1099  double *diffL, const int dim) =
1101 
1102  int nb_gauss_pts = pts.size2();
1103 
1104  // edges
1105  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1106  int sense[6], order[6];
1107  if (data.dataOnEntities[MBEDGE].size() != 6) {
1108  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1109  }
1110  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1111  for (int ee = 0; ee != 6; ee++) {
1112  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1113  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1114  "data inconsistency");
1115  }
1116  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1117  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1118  int nb_dofs = NBEDGE_AINSWORTH_HCURL(
1119  data.dataOnEntities[MBEDGE][ee].getOrder());
1120  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1121  3 * nb_dofs, false);
1122  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1123  9 * nb_dofs, false);
1124  hcurl_edge_n[ee] =
1125  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1126  diff_hcurl_edge_n[ee] =
1127  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1128  }
1130  sense, order,
1131  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1132  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1133  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
1134  } else {
1135  for (int ee = 0; ee != 6; ee++) {
1136  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1137  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1138  false);
1139  }
1140  }
1141 
1142  // triangles
1143  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1144  int order[4];
1145  // faces
1146  if (data.dataOnEntities[MBTRI].size() != 4) {
1147  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1148  }
1149  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1150  for (int ff = 0; ff != 4; ff++) {
1151  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1152  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1153  "data inconsistency");
1154  }
1155  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1156  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
1157  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1158  3 * nb_dofs, false);
1159  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1160  9 * nb_dofs, false);
1161  hcurl_base_n[ff] =
1162  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1163  diff_hcurl_base_n[ff] =
1164  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1165  }
1166  if (data.facesNodes.size1() != 4) {
1167  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1168  }
1169  if (data.facesNodes.size2() != 3) {
1170  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1171  }
1173  &*data.facesNodes.data().begin(), order,
1174  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1175  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1176  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
1177  } else {
1178  for (int ff = 0; ff != 4; ff++) {
1179  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1180  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1181  false);
1182  }
1183  }
1184 
1185  if (data.spacesOnEntities[MBTET].test(HCURL)) {
1186 
1187  // volume
1188  int order = data.dataOnEntities[MBTET][0].getOrder();
1189  int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
1190  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1191  3 * nb_vol_dofs, false);
1192  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1193  9 * nb_vol_dofs, false);
1195  data.dataOnEntities[MBTET][0].getOrder(),
1196  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1197  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1198  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1199  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1200  nb_gauss_pts, base_polynomials);
1201 
1202  } else {
1203  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1204  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1205  }
1206 
1208 }
#define NBVOLUMETET_AINSWORTH_HCURL(P)
#define NBFACETRI_AINSWORTH_HCURL(P)
#define NBEDGE_AINSWORTH_HCURL(P)
MoFEMErrorCode Hcurl_Ainsworth_VolumeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H-curl volume base functions.
Definition: Hcurl.cpp:1403
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Edge based H-curl base functions on tetrahedral.
Definition: Hcurl.cpp:16
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

◆ getValueHcurlDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlDemkowiczBase ( MatrixDouble pts)
private

Definition at line 1211 of file TetPolynomialBase.cpp.

1211  {
1213 
1214  EntitiesFieldData &data = cTx->dAta;
1215  const FieldApproximationBase base = cTx->bAse;
1216  if (base != DEMKOWICZ_JACOBI_BASE) {
1217  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1218  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1219  "but base is %s",
1220  ApproximationBaseNames[base]);
1221  }
1222 
1223  int nb_gauss_pts = pts.size2();
1224 
1225  // edges
1226  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1227  int sense[6], order[6];
1228  if (data.dataOnEntities[MBEDGE].size() != 6) {
1229  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1230  "wrong size of data structure, expected space for six edges "
1231  "but is %d",
1232  data.dataOnEntities[MBEDGE].size());
1233  }
1234  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1235  for (int ee = 0; ee != 6; ee++) {
1236  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1237  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1238  "orintation of edges is not set");
1239  }
1240  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1241  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1242  int nb_dofs = NBEDGE_DEMKOWICZ_HCURL(
1243  data.dataOnEntities[MBEDGE][ee].getOrder());
1244  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1245  3 * nb_dofs, false);
1246  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1247  9 * nb_dofs, false);
1248  hcurl_edge_n[ee] =
1249  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1250  diff_hcurl_edge_n[ee] =
1251  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1252  }
1254  sense, order,
1255  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1256  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1257  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
1258  } else {
1259  // No DOFs on edges, resize base function matrices, indicating that no
1260  // dofs on them.
1261  for (int ee = 0; ee != 6; ee++) {
1262  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1263  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1264  false);
1265  }
1266  }
1267 
1268  // triangles
1269  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1270  int order[4];
1271  // faces
1272  if (data.dataOnEntities[MBTRI].size() != 4) {
1273  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1274  "data structure for storing face h-curl base have wrong size "
1275  "should be four but is %d",
1276  data.dataOnEntities[MBTRI].size());
1277  }
1278  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1279  for (int ff = 0; ff != 4; ff++) {
1280  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1281  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1282  "orintation of face is not set");
1283  }
1284  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1285  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order[ff]);
1286  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1287  3 * nb_dofs, false);
1288  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1289  9 * nb_dofs, false);
1290  hcurl_base_n[ff] =
1291  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1292  diff_hcurl_base_n[ff] =
1293  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1294  }
1295  if (data.facesNodes.size1() != 4) {
1296  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1297  "data inconsistency, should be four faces");
1298  }
1299  if (data.facesNodes.size2() != 3) {
1300  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1301  "data inconsistency, should be three nodes on face");
1302  }
1304  &*data.facesNodes.data().begin(), order,
1305  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1306  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1307  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
1308  } else {
1309  // No DOFs on faces, resize base function matrices, indicating that no
1310  // dofs on them.
1311  for (int ff = 0; ff != 4; ff++) {
1312  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1313  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1314  false);
1315  }
1316  }
1317 
1318  if (data.spacesOnEntities[MBTET].test(HCURL)) {
1319  // volume
1320  int order = data.dataOnEntities[MBTET][0].getOrder();
1321  int nb_vol_dofs = NBVOLUMETET_DEMKOWICZ_HCURL(order);
1322  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1323  3 * nb_vol_dofs, false);
1324  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1325  9 * nb_vol_dofs, false);
1327  data.dataOnEntities[MBTET][0].getOrder(),
1328  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1329  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1330  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1331  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1332  nb_gauss_pts);
1333  } else {
1334  // No DOFs on faces, resize base function matrices, indicating that no
1335  // dofs on them.
1336  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1337  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1338  }
1339 
1341 }
#define NBEDGE_DEMKOWICZ_HCURL(P)
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
#define NBFACETRI_DEMKOWICZ_HCURL(P)
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
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
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

◆ getValueHdiv()

MoFEMErrorCode TetPolynomialBase::getValueHdiv ( MatrixDouble pts)
private

Get base functions for Hdiv space.

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

Definition at line 1076 of file TetPolynomialBase.cpp.

1076  {
1078 
1079  switch (cTx->bAse) {
1082  return getValueHdivAinsworthBase(pts);
1083  case DEMKOWICZ_JACOBI_BASE:
1084  return getValueHdivDemkowiczBase(pts);
1085  default:
1086  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1087  }
1088 
1090 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)

◆ getValueHdivAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBase ( MatrixDouble pts)
private

Definition at line 669 of file TetPolynomialBase.cpp.

669  {
671 
672  EntitiesFieldData &data = cTx->dAta;
673  const FieldApproximationBase base = cTx->bAse;
674  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
675  double *diffL, const int dim) =
677 
678  int nb_gauss_pts = pts.size2();
679 
680  // face shape functions
681 
682  double *phi_f_e[4][3];
683  double *phi_f[4];
684  double *diff_phi_f_e[4][3];
685  double *diff_phi_f[4];
686 
687  N_face_edge.resize(4, 3, false);
688  N_face_bubble.resize(4, false);
689  diffN_face_edge.resize(4, 3, false);
690  diffN_face_bubble.resize(4, false);
691 
692  int faces_order[4];
693  for (int ff = 0; ff != 4; ++ff) {
694  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
695  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
696  }
697  faces_order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
698  // three edges on face
699  for (int ee = 0; ee < 3; ee++) {
700  N_face_edge(ff, ee).resize(
701  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
702  false);
703  diffN_face_edge(ff, ee).resize(
704  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
705  false);
706  phi_f_e[ff][ee] = &*N_face_edge(ff, ee).data().begin();
707  diff_phi_f_e[ff][ee] = &*diffN_face_edge(ff, ee).data().begin();
708  }
709  N_face_bubble[ff].resize(nb_gauss_pts,
710  3 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
711  false);
712  diffN_face_bubble[ff].resize(
713  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
714  false);
715  phi_f[ff] = &*(N_face_bubble[ff].data().begin());
716  diff_phi_f[ff] = &*(diffN_face_bubble[ff].data().begin());
717  }
718 
720  &data.facesNodes(0, 0), faces_order,
721  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
722  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f_e,
723  diff_phi_f_e, nb_gauss_pts, base_polynomials);
724 
726  &data.facesNodes(0, 0), faces_order,
727  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
728  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f,
729  diff_phi_f, nb_gauss_pts, base_polynomials);
730 
731  // volume shape functions
732 
733  double *phi_v_e[6];
734  double *phi_v_f[4];
735  double *phi_v;
736  double *diff_phi_v_e[6];
737  double *diff_phi_v_f[4];
738  double *diff_phi_v;
739 
740  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
741 
742  N_volume_edge.resize(6, false);
743  diffN_volume_edge.resize(6, false);
744  for (int ee = 0; ee != 6; ++ee) {
745  N_volume_edge[ee].resize(
746  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
747  diffN_volume_edge[ee].resize(
748  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
749  phi_v_e[ee] = &*(N_volume_edge[ee].data().begin());
750  diff_phi_v_e[ee] = &*(diffN_volume_edge[ee].data().begin());
751  }
753  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
754  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_e,
755  diff_phi_v_e, nb_gauss_pts, base_polynomials);
756 
757  N_volume_face.resize(4, false);
758  diffN_volume_face.resize(4, false);
759  for (int ff = 0; ff != 4; ++ff) {
760  N_volume_face[ff].resize(
761  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
762  diffN_volume_face[ff].resize(
763  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
764  phi_v_f[ff] = &*(N_volume_face[ff].data().begin());
765  diff_phi_v_f[ff] = &*(diffN_volume_face[ff].data().begin());
766  }
768  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
769  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_f,
770  diff_phi_v_f, nb_gauss_pts, base_polynomials);
771 
772  N_volume_bubble.resize(
773  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
774  diffN_volume_bubble.resize(
775  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
776  phi_v = &*(N_volume_bubble.data().begin());
777  diff_phi_v = &*(diffN_volume_bubble.data().begin());
779  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
780  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v, diff_phi_v,
781  nb_gauss_pts, base_polynomials);
782 
783  // Set shape functions into data structure Shape functions hast to be put
784  // in arrays in order which guarantee hierarchical series of degrees of
785  // freedom, i.e. in other words dofs form sub-entities has to be group
786  // by order.
787 
790 
791  // faces
792  if (data.dataOnEntities[MBTRI].size() != 4) {
793  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
794  }
795  for (int ff = 0; ff != 4; ff++) {
796  data.dataOnEntities[MBTRI][ff].getN(base).resize(
797  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
798  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
799  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
800  if (NBFACETRI_AINSWORTH_HDIV(faces_order[ff]) == 0)
801  continue;
802  // face
803  double *base_ptr =
804  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
805  FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
806  &base_ptr[HVEC2], 3);
807  double *diff_base_ptr =
808  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
810  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
811  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
812  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
813  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
814  &diff_base_ptr[HVEC2_2], 9);
815  // face-face
816  boost::shared_ptr<FTensor::Tensor1<double *, 3>> t_base_f;
817  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_base_f;
818  if (NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]) > 0) {
819  base_ptr = phi_f[ff];
820  t_base_f = boost::shared_ptr<FTensor::Tensor1<double *, 3>>(
821  new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
822  &base_ptr[HVEC2], 3));
823  diff_base_ptr = diff_phi_f[ff];
824  t_diff_base_f = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
826  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
827  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
828  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
829  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
830  &diff_base_ptr[HVEC2_2], 9));
831  }
832  // edge-face
833  base_ptr = phi_f_e[ff][0];
834  FTensor::Tensor1<double *, 3> t_base_f_e0(base_ptr, &base_ptr[HVEC1],
835  &base_ptr[HVEC2], 3);
836  diff_base_ptr = diff_phi_f_e[ff][0];
837  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e0(
838  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
839  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
840  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
841  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
842  &diff_base_ptr[HVEC2_2], 9);
843  base_ptr = phi_f_e[ff][1];
844  FTensor::Tensor1<double *, 3> t_base_f_e1(base_ptr, &base_ptr[HVEC1],
845  &base_ptr[HVEC2], 3);
846  diff_base_ptr = diff_phi_f_e[ff][1];
847  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e1(
848  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
849  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
850  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
851  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
852  &diff_base_ptr[HVEC2_2], 9);
853  base_ptr = phi_f_e[ff][2];
854  FTensor::Tensor1<double *, 3> t_base_f_e2(base_ptr, &base_ptr[HVEC1],
855  &base_ptr[HVEC2], 3);
856  diff_base_ptr = diff_phi_f_e[ff][2];
857  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e2(
858  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
859  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
860  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
861  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
862  &diff_base_ptr[HVEC2_2], 9);
863  for (int gg = 0; gg != nb_gauss_pts; gg++) {
864  for (int oo = 0; oo != faces_order[ff]; oo++) {
865  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
866  dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
867  t_base(i) = t_base_f_e0(i);
868  ++t_base;
869  ++t_base_f_e0;
870  t_diff_base(i, j) = t_diff_base_f_e0(i, j);
871  ++t_diff_base;
872  ++t_diff_base_f_e0;
873  t_base(i) = t_base_f_e1(i);
874  ++t_base;
875  ++t_base_f_e1;
876  t_diff_base(i, j) = t_diff_base_f_e1(i, j);
877  ++t_diff_base;
878  ++t_diff_base_f_e1;
879  t_base(i) = t_base_f_e2(i);
880  ++t_base;
881  ++t_base_f_e2;
882  t_diff_base(i, j) = t_diff_base_f_e2(i, j);
883  ++t_diff_base;
884  ++t_diff_base_f_e2;
885  }
886  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
887  dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
888  t_base(i) = (*t_base_f)(i);
889  ++t_base;
890  ++(*t_base_f);
891  t_diff_base(i, j) = (*t_diff_base_f)(i, j);
892  ++t_diff_base;
893  ++(*t_diff_base_f);
894  }
895  }
896  }
897  }
898 
899  // volume
900  data.dataOnEntities[MBTET][0].getN(base).resize(
901  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
902  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
903  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
904  if (NBVOLUMETET_AINSWORTH_HDIV(volume_order) > 0) {
905  double *base_ptr =
906  &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
907  FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
908  &base_ptr[HVEC2], 3);
909  double *diff_base_ptr =
910  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
912  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
913  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
914  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
915  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
916  &diff_base_ptr[HVEC2_2], 9);
917  // edges
918  std::vector<FTensor::Tensor1<double *, 3>> t_base_v_e;
919  t_base_v_e.reserve(6);
920  std::vector<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v_e;
921  t_diff_base_v_e.reserve(6);
922  for (int ee = 0; ee != 6; ee++) {
923  base_ptr = phi_v_e[ee];
924  diff_base_ptr = diff_phi_v_e[ee];
925  t_base_v_e.push_back(FTensor::Tensor1<double *, 3>(
926  base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
927  t_diff_base_v_e.push_back(FTensor::Tensor2<double *, 3, 3>(
928  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
929  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
930  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
931  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
932  &diff_base_ptr[HVEC2_2], 9));
933  }
934  // faces
935  std::vector<FTensor::Tensor1<double *, 3>> t_base_v_f;
936  t_base_v_f.reserve(4);
937  std::vector<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v_f;
938  t_diff_base_v_f.reserve(4);
939  if (NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order) > 0) {
940  for (int ff = 0; ff != 4; ff++) {
941  base_ptr = phi_v_f[ff];
942  diff_base_ptr = diff_phi_v_f[ff];
943  t_base_v_f.push_back(FTensor::Tensor1<double *, 3>(
944  base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
945  t_diff_base_v_f.push_back(FTensor::Tensor2<double *, 3, 3>(
946  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
947  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
948  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
949  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
950  &diff_base_ptr[HVEC2_2], 9));
951  }
952  }
953  boost::shared_ptr<FTensor::Tensor1<double *, 3>> t_base_v;
954  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v;
955  if (NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order) > 0) {
956  base_ptr = phi_v;
957  t_base_v = boost::shared_ptr<FTensor::Tensor1<double *, 3>>(
958  new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
959  &base_ptr[HVEC2], 3));
960  diff_base_ptr = diff_phi_v;
961  t_diff_base_v = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
963  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
964  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
965  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
966  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
967  &diff_base_ptr[HVEC2_2], 9));
968  }
969  for (int gg = 0; gg != nb_gauss_pts; gg++) {
970  for (int oo = 0; oo < volume_order; oo++) {
971  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
972  dd < NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
973  for (int ee = 0; ee < 6; ee++) {
974  t_base(i) = t_base_v_e[ee](i);
975  ++t_base;
976  ++t_base_v_e[ee];
977  t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
978  ++t_diff_base;
979  ++t_diff_base_v_e[ee];
980  }
981  }
982  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
983  dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
984  for (int ff = 0; ff < 4; ff++) {
985  t_base(i) = t_base_v_f[ff](i);
986  ++t_base;
987  ++t_base_v_f[ff];
988  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
989  ++t_diff_base;
990  ++t_diff_base_v_f[ff];
991  }
992  }
993  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
994  dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
995  t_base(i) = (*t_base_v)(i);
996  ++t_base;
997  ++(*t_base_v);
998  t_diff_base(i, j) = (*t_diff_base_v)(i, j);
999  ++t_diff_base;
1000  ++(*t_diff_base_v);
1001  }
1002  }
1003  }
1004  }
1005 
1007 }
@ HVEC1
Definition: definitions.h:199
@ HVEC2
Definition: definitions.h:199
@ HVEC1_1
Definition: definitions.h:209
@ HVEC0_1
Definition: definitions.h:208
@ HVEC1_0
Definition: definitions.h:206
@ HVEC2_1
Definition: definitions.h:210
@ HVEC1_2
Definition: definitions.h:212
@ HVEC2_2
Definition: definitions.h:213
@ HVEC2_0
Definition: definitions.h:207
@ HVEC0_2
Definition: definitions.h:211
@ HVEC0_0
Definition: definitions.h:205
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)
#define NBVOLUMETET_AINSWORTH_HDIV(P)
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
#define NBFACETRI_AINSWORTH_HDIV(P)
FTensor::Index< 'j', 3 > j
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
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[], double *diff_phi_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
Definition: Hdiv.cpp:137
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3], double *diff_phi_f_e[4][3], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
Definition: Hdiv.cpp:13
MoFEMErrorCode Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: Hdiv.cpp:368
MoFEMErrorCode Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_e[6], double *diff_phi_v_e[6], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base function, Edge-based interior (volume) functions by Ainsworth .
Definition: Hdiv.cpp:280
MoFEMErrorCode Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Interior bubble functions by Ainsworth .
Definition: Hdiv.cpp:484
ublas::matrix< MatrixDouble > diffN_face_edge
ublas::vector< MatrixDouble > diffN_volume_face
ublas::vector< MatrixDouble > diffN_face_bubble
ublas::vector< MatrixDouble > N_volume_edge
ublas::vector< MatrixDouble > N_volume_face
ublas::vector< MatrixDouble > N_face_bubble
ublas::matrix< MatrixDouble > N_face_edge
ublas::vector< MatrixDouble > diffN_volume_edge

◆ getValueHdivDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBase ( MatrixDouble pts)
private

Definition at line 1009 of file TetPolynomialBase.cpp.

1009  {
1011 
1012  EntitiesFieldData &data = cTx->dAta;
1013  const FieldApproximationBase base = cTx->bAse;
1014  if (base != DEMKOWICZ_JACOBI_BASE) {
1015  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1016  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1017  "but base is %s",
1018  ApproximationBaseNames[base]);
1019  }
1020  int nb_gauss_pts = pts.size2();
1021 
1022  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1023 
1024  int p_f[4];
1025  double *phi_f[4];
1026  double *diff_phi_f[4];
1027 
1028  // Calculate base function on tet faces
1029  for (int ff = 0; ff != 4; ff++) {
1030  int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1031  int order = volume_order > face_order ? volume_order : face_order;
1032  data.dataOnEntities[MBTRI][ff].getN(base).resize(
1033  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1034  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1035  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1036  p_f[ff] = order;
1037  phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1038  diff_phi_f[ff] =
1039  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1040  if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
1041  continue;
1043  &data.facesNodes(ff, 0), order,
1044  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1045  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1046  phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1047  }
1048 
1049  // Calculate base functions in tet interior
1050  if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
1051  data.dataOnEntities[MBTET][0].getN(base).resize(
1052  nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1053  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
1054  nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1055  double *phi_v = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1056  double *diff_phi_v =
1057  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1059  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1060  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
1061  diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
1062  }
1063 
1064  // Set size of face base correctly
1065  for (int ff = 0; ff != 4; ff++) {
1066  int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1067  data.dataOnEntities[MBTRI][ff].getN(base).resize(
1068  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1069  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1070  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1071  }
1072 
1074 }
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
#define NBFACETRI_DEMKOWICZ_HDIV(P)
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition: Hdiv.cpp:617
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
Definition: Hdiv.cpp:762

◆ getValueL2()

MoFEMErrorCode TetPolynomialBase::getValueL2 ( MatrixDouble pts)
private

Get base functions for L2 space.

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

Definition at line 445 of file TetPolynomialBase.cpp.

445  {
447 
448  switch (cTx->bAse) {
452  break;
455  break;
456  default:
457  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
458  }
459 
461 }
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)

◆ getValueL2AinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueL2AinsworthBase ( MatrixDouble pts)
private

Definition at line 463 of file TetPolynomialBase.cpp.

463  {
465 
466  EntitiesFieldData &data = cTx->dAta;
467  const FieldApproximationBase base = cTx->bAse;
468  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
469  double *diffL, const int dim) =
471 
472  int nb_gauss_pts = pts.size2();
473 
474  data.dataOnEntities[MBTET][0].getN(base).resize(
475  nb_gauss_pts,
476  NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()), false);
477  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
478  nb_gauss_pts,
479  3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()), false);
480 
482  data.dataOnEntities[MBTET][0].getOrder(),
483  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
484  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
485  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
486  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
487  nb_gauss_pts, base_polynomials);
488 
490 }
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on tetrahedron for L2 space.
Definition: l2.c:84
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.

◆ getValueL2BernsteinBezierBase()

MoFEMErrorCode TetPolynomialBase::getValueL2BernsteinBezierBase ( MatrixDouble pts)
private

Definition at line 493 of file TetPolynomialBase.cpp.

493  {
495 
496  EntitiesFieldData &data = cTx->dAta;
497  const std::string field_name = cTx->fieldName;
498  const int nb_gauss_pts = pts.size2();
499 
500  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
501  (unsigned int)nb_gauss_pts)
502  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
503  "Base functions or nodes has wrong number of integration points "
504  "for base %s",
506  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
507 
508  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
509  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
510  if (!ptr)
511  ptr.reset(new MatrixInt());
512  return *ptr;
513  };
514 
515  auto get_base = [field_name](auto &data) -> MatrixDouble & {
516  auto &ptr = data.getBBNSharedPtr(field_name);
517  if (!ptr)
518  ptr.reset(new MatrixDouble());
519  return *ptr;
520  };
521 
522  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
523  auto &ptr = data.getBBDiffNSharedPtr(field_name);
524  if (!ptr)
525  ptr.reset(new MatrixDouble());
526  return *ptr;
527  };
528 
529  auto get_alpha_by_name_ptr =
530  [](auto &data,
531  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
532  return data.getBBAlphaIndicesSharedPtr(field_name);
533  };
534 
535  auto get_base_by_name_ptr =
536  [](auto &data,
537  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
538  return data.getBBNSharedPtr(field_name);
539  };
540 
541  auto get_diff_base_by_name_ptr =
542  [](auto &data,
543  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
544  return data.getBBDiffNSharedPtr(field_name);
545  };
546 
547  auto get_alpha_by_order_ptr =
548  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
549  return data.getBBAlphaIndicesByOrderSharedPtr(o);
550  };
551 
552  auto get_base_by_order_ptr =
553  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
554  return data.getBBNByOrderSharedPtr(o);
555  };
556 
557  auto get_diff_base_by_order_ptr =
558  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
559  return data.getBBDiffNByOrderSharedPtr(o);
560  };
561 
562  if (data.spacesOnEntities[MBTET].test(L2)) {
563  if (data.dataOnEntities[MBTET].size() != 1)
564  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
565  "Wrong size ent of ent data");
566 
567  auto &ent_data = data.dataOnEntities[MBTET][0];
568  const int order = ent_data.getOrder();
569  const int nb_dofs = NBVOLUMETET_L2(order);
570 
571  if (get_alpha_by_order_ptr(ent_data, order)) {
572  get_alpha_by_name_ptr(ent_data, field_name) =
573  get_alpha_by_order_ptr(ent_data, order);
574  get_base_by_name_ptr(ent_data, field_name) =
575  get_base_by_order_ptr(ent_data, order);
576  get_diff_base_by_name_ptr(ent_data, field_name) =
577  get_diff_base_by_order_ptr(ent_data, order);
578  } else {
579 
580  auto &get_n = get_base(ent_data);
581  auto &get_diff_n = get_diff_base(ent_data);
582  get_n.resize(nb_gauss_pts, nb_dofs, false);
583  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
584 
585  if (nb_dofs) {
586 
587  if (order == 0) {
588 
589  if (nb_dofs != 1)
590  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
591  "Inconsistent number of DOFs");
592 
593  auto &tri_alpha = get_alpha(ent_data);
594  tri_alpha.clear();
595  get_n(0, 0) = 1;
596  get_diff_n.clear();
597 
598  } else {
599 
600  if (nb_dofs != 4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order) +
602  nb_dofs != NBVOLUMETET_L2(order))
603  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
604  "Inconsistent number of DOFs");
605 
606  auto &tet_alpha = get_alpha(ent_data);
607  tet_alpha.resize(nb_dofs, 4, false);
608 
610  &tet_alpha(0, 0));
611  if (order > 1) {
612  std::array<int, 6> edge_n{order, order, order, order, order, order};
613  std::array<int *, 6> tet_edge_ptr{
614  &tet_alpha(4, 0),
615  &tet_alpha(4 + 1 * NBEDGE_H1(order), 0),
616  &tet_alpha(4 + 2 * NBEDGE_H1(order), 0),
617  &tet_alpha(4 + 3 * NBEDGE_H1(order), 0),
618  &tet_alpha(4 + 4 * NBEDGE_H1(order), 0),
619  &tet_alpha(4 + 5 * NBEDGE_H1(order), 0)};
621  tet_edge_ptr.data());
622  if (order > 2) {
623  std::array<int, 6> face_n{order, order, order, order};
624  std::array<int *, 6> tet_face_ptr{
625  &tet_alpha(4 + 6 * NBEDGE_H1(order), 0),
626  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 1 * NBFACETRI_H1(order),
627  0),
628  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 2 * NBFACETRI_H1(order),
629  0),
630  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 3 * NBFACETRI_H1(order),
631  0),
632  };
634  face_n.data(), tet_face_ptr.data());
635  if (order > 3)
637  order,
638  &tet_alpha(
639  4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order), 0));
640  }
641  }
642 
644  order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
645  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
646  &get_diff_n(0, 0));
647 
648  get_alpha_by_order_ptr(ent_data, order) =
649  get_alpha_by_name_ptr(ent_data, field_name);
650  get_base_by_order_ptr(ent_data, order) =
651  get_base_by_name_ptr(ent_data, field_name);
652  get_diff_base_by_order_ptr(ent_data, order) =
653  get_diff_base_by_name_ptr(ent_data, field_name);
654  }
655  }
656  }
657  } else {
658  auto &ent_data = data.dataOnEntities[MBTET][0];
659  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
660  auto &get_n = get_base(ent_data);
661  auto &get_diff_n = get_diff_base(ent_data);
662  get_n.resize(nb_gauss_pts, 0, false);
663  get_diff_n.resize(nb_gauss_pts, 0, false);
664  }
665 
667 }
static MoFEMErrorCode generateIndicesVertexTet(const int N, int *alpha)

◆ query_interface()

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

Reimplemented from MoFEM::BaseFunction.

Definition at line 25 of file TetPolynomialBase.cpp.

26  {
27 
29  *iface = const_cast<TetPolynomialBase *>(this);
31 }
Calculate base functions on tetrahedral.

Member Data Documentation

◆ cTx

EntPolynomialBaseCtx* MoFEM::TetPolynomialBase::cTx
private

Definition at line 43 of file TetPolynomialBase.hpp.

◆ diffN_face_bubble

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

Definition at line 96 of file TetPolynomialBase.hpp.

◆ diffN_face_edge

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

Definition at line 95 of file TetPolynomialBase.hpp.

◆ diffN_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::diffN_volume_bubble
private

Definition at line 99 of file TetPolynomialBase.hpp.

◆ diffN_volume_edge

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

Definition at line 97 of file TetPolynomialBase.hpp.

◆ diffN_volume_face

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

Definition at line 98 of file TetPolynomialBase.hpp.

◆ N_face_bubble

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

Definition at line 90 of file TetPolynomialBase.hpp.

◆ N_face_edge

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

Definition at line 89 of file TetPolynomialBase.hpp.

◆ N_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::N_volume_bubble
private

Definition at line 93 of file TetPolynomialBase.hpp.

◆ N_volume_edge

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

Definition at line 91 of file TetPolynomialBase.hpp.

◆ N_volume_face

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

Definition at line 92 of file TetPolynomialBase.hpp.

◆ senseFaceAlpha

MatrixInt MoFEM::TetPolynomialBase::senseFaceAlpha
private

Definition at line 114 of file TetPolynomialBase.hpp.


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