v0.13.1
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
MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, MoFEM::UnknownInterface **iface) const
 
virtual MoFEMErrorCode getValue (MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
virtual MoFEMErrorCode getValue (MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
- Public Member Functions inherited from MoFEM::BaseFunctionUnknownInterface
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
virtual ~BaseFunctionUnknownInterface ()=default
 
- Public Member Functions inherited from MoFEM::UnknownInterface
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
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 19 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 1350 of file TetPolynomialBase.cpp.

1351 {
1353
1355
1356 int nb_gauss_pts = pts.size2();
1357 if (!nb_gauss_pts)
1359
1360 if (pts.size1() < 3)
1361 SETERRQ(
1362 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1363 "Wrong dimension of pts, should be at least 3 rows with coordinates");
1364
1366 const FieldApproximationBase base = cTx->bAse;
1367 EntitiesFieldData &data = cTx->dAta;
1368 if (cTx->copyNodeBase == LASTBASE) {
1369 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
1370 false);
1372 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1373 &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
1374 } else {
1375 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
1376 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
1377 }
1378 if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
1379 (unsigned int)nb_gauss_pts) {
1380 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1381 "Base functions or nodes has wrong number of integration points "
1382 "for base %s",
1384 }
1385 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
1386 std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
1387 data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1388 }
1389
1390 switch (cTx->sPace) {
1391 case H1:
1392 CHKERR getValueH1(pts);
1393 break;
1394 case HDIV:
1395 CHKERR getValueHdiv(pts);
1396 break;
1397 case HCURL:
1398 CHKERR getValueHcurl(pts);
1399 break;
1400 case L2:
1401 CHKERR getValueL2(pts);
1402 break;
1403 default:
1404 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space");
1405 }
1406
1408}
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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:673
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
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 21 of file TetPolynomialBase.cpp.

21 {
23
24 switch (cTx->bAse) {
28 break;
31 break;
32 default:
33 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
34 }
35
37}
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)

◆ getValueH1AinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueH1AinsworthBase ( MatrixDouble pts)
private

Definition at line 39 of file TetPolynomialBase.cpp.

39 {
41
42 EntitiesFieldData &data = cTx->dAta;
43 const FieldApproximationBase base = cTx->bAse;
44 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
45 double *diffL, const int dim) =
47
48 int nb_gauss_pts = pts.size2();
49
50 int sense[6], order[6];
51 if (data.spacesOnEntities[MBEDGE].test(H1)) {
52 // edges
53 if (data.dataOnEntities[MBEDGE].size() != 6) {
54 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
55 }
56 double *h1_edge_n[6], *diff_h1_egde_n[6];
57 for (int ee = 0; ee != 6; ++ee) {
58 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
59 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
60 "data inconsistency");
61 }
62 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
63 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
64 int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
65 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
66 false);
67 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
68 3 * nb_dofs, false);
69 h1_edge_n[ee] =
70 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
71 diff_h1_egde_n[ee] =
72 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
73 }
75 sense, order,
76 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
77 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
78 h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
79 } else {
80 for (int ee = 0; ee != 6; ++ee) {
81 data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
82 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
83 }
84 }
85
86 if (data.spacesOnEntities[MBTRI].test(H1)) {
87 // faces
88 if (data.dataOnEntities[MBTRI].size() != 4) {
89 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
90 }
91 double *h1_face_n[4], *diff_h1_face_n[4];
92 for (int ff = 0; ff != 4; ++ff) {
93 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
94 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
95 "data inconsistency");
96 }
97 int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getOrder());
98 order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
99 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
100 false);
101 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
102 3 * nb_dofs, false);
103 h1_face_n[ff] =
104 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
105 diff_h1_face_n[ff] =
106 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
107 }
108 if (data.facesNodes.size1() != 4) {
109 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
110 }
111 if (data.facesNodes.size2() != 3) {
112 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
113 }
115 &*data.facesNodes.data().begin(), order,
116 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
117 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
118 h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
119
120 } else {
121 for (int ff = 0; ff != 4; ++ff) {
122 data.dataOnEntities[MBTRI][ff].getN(base).resize(0, false);
123 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, false);
124 }
125 }
126
127 if (data.spacesOnEntities[MBTET].test(H1)) {
128 // volume
129 int order = data.dataOnEntities[MBTET][0].getOrder();
130 int nb_vol_dofs = NBVOLUMETET_H1(order);
131 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
132 false);
133 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
134 3 * nb_vol_dofs, false);
136 data.dataOnEntities[MBTET][0].getOrder(),
137 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
138 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
139 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
140 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
141 nb_gauss_pts, base_polynomials);
142 } else {
143 data.dataOnEntities[MBTET][0].getN(base).resize(0, 0, false);
144 data.dataOnEntities[MBTET][0].getDiffN(base).resize(0, 0, false);
145 }
146
148}
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:274
#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:475
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
#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 151 of file TetPolynomialBase.cpp.

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

1331 {
1333
1334 switch (cTx->bAse) {
1338 break;
1341 break;
1342 default:
1343 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1344 }
1345
1347}
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)

◆ getValueHcurlAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlAinsworthBase ( MatrixDouble pts)
private

Definition at line 1081 of file TetPolynomialBase.cpp.

1081 {
1083
1084 EntitiesFieldData &data = cTx->dAta;
1085 const FieldApproximationBase base = cTx->bAse;
1086 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
1087 double *diffL, const int dim) =
1089
1090 int nb_gauss_pts = pts.size2();
1091
1092 // edges
1093 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1094 int sense[6], order[6];
1095 if (data.dataOnEntities[MBEDGE].size() != 6) {
1096 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1097 }
1098 double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1099 for (int ee = 0; ee != 6; ee++) {
1100 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1101 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1102 "data inconsistency");
1103 }
1104 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1105 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1106 int nb_dofs = NBEDGE_AINSWORTH_HCURL(
1107 data.dataOnEntities[MBEDGE][ee].getOrder());
1108 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1109 3 * nb_dofs, false);
1110 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1111 9 * nb_dofs, false);
1112 hcurl_edge_n[ee] =
1113 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1114 diff_hcurl_edge_n[ee] =
1115 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1116 }
1118 sense, order,
1119 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1120 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1121 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
1122 } else {
1123 for (int ee = 0; ee != 6; ee++) {
1124 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1125 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1126 false);
1127 }
1128 }
1129
1130 // triangles
1131 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1132 int order[4];
1133 // faces
1134 if (data.dataOnEntities[MBTRI].size() != 4) {
1135 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1136 }
1137 double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1138 for (int ff = 0; ff != 4; ff++) {
1139 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1140 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1141 "data inconsistency");
1142 }
1143 order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1144 int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
1145 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1146 3 * nb_dofs, false);
1147 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1148 9 * nb_dofs, false);
1149 hcurl_base_n[ff] =
1150 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1151 diff_hcurl_base_n[ff] =
1152 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1153 }
1154 if (data.facesNodes.size1() != 4) {
1155 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1156 }
1157 if (data.facesNodes.size2() != 3) {
1158 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1159 }
1161 &*data.facesNodes.data().begin(), order,
1162 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1163 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1164 hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
1165 } else {
1166 for (int ff = 0; ff != 4; ff++) {
1167 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1168 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1169 false);
1170 }
1171 }
1172
1173 if (data.spacesOnEntities[MBTET].test(HCURL)) {
1174
1175 // volume
1176 int order = data.dataOnEntities[MBTET][0].getOrder();
1177 int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
1178 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1179 3 * nb_vol_dofs, false);
1180 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1181 9 * nb_vol_dofs, false);
1183 data.dataOnEntities[MBTET][0].getOrder(),
1184 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1185 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1186 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1187 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1188 nb_gauss_pts, base_polynomials);
1189
1190 } else {
1191 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1192 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1193 }
1194
1196}
#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 1199 of file TetPolynomialBase.cpp.

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

1064 {
1066
1067 switch (cTx->bAse) {
1070 return getValueHdivAinsworthBase(pts);
1072 return getValueHdivDemkowiczBase(pts);
1073 default:
1074 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1075 }
1076
1078}
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)

◆ getValueHdivAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBase ( MatrixDouble pts)
private

Definition at line 657 of file TetPolynomialBase.cpp.

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

997 {
999
1000 EntitiesFieldData &data = cTx->dAta;
1001 const FieldApproximationBase base = cTx->bAse;
1002 if (base != DEMKOWICZ_JACOBI_BASE) {
1003 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1004 "This should be used only with DEMKOWICZ_JACOBI_BASE "
1005 "but base is %s",
1007 }
1008 int nb_gauss_pts = pts.size2();
1009
1010 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1011
1012 int p_f[4];
1013 double *phi_f[4];
1014 double *diff_phi_f[4];
1015
1016 // Calculate base function on tet faces
1017 for (int ff = 0; ff != 4; ff++) {
1018 int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1019 int order = volume_order > face_order ? volume_order : face_order;
1020 data.dataOnEntities[MBTRI][ff].getN(base).resize(
1021 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1022 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1023 nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1024 p_f[ff] = order;
1025 phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1026 diff_phi_f[ff] =
1027 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1029 continue;
1031 &data.facesNodes(ff, 0), order,
1032 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1033 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1034 phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1035 }
1036
1037 // Calculate base functions in tet interior
1038 if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
1039 data.dataOnEntities[MBTET][0].getN(base).resize(
1040 nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1041 data.dataOnEntities[MBTET][0].getDiffN(base).resize(
1042 nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1043 double *phi_v = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1044 double *diff_phi_v =
1045 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1047 volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1048 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
1049 diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
1050 }
1051
1052 // Set size of face base correctly
1053 for (int ff = 0; ff != 4; ff++) {
1054 int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1055 data.dataOnEntities[MBTRI][ff].getN(base).resize(
1056 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1057 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1058 nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1059 }
1060
1062}
#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 433 of file TetPolynomialBase.cpp.

433 {
435
436 switch (cTx->bAse) {
440 break;
443 break;
444 default:
445 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
446 }
447
449}
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)

◆ getValueL2AinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueL2AinsworthBase ( MatrixDouble pts)
private

Definition at line 451 of file TetPolynomialBase.cpp.

451 {
453
454 EntitiesFieldData &data = cTx->dAta;
455 const FieldApproximationBase base = cTx->bAse;
456 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
457 double *diffL, const int dim) =
459
460 int nb_gauss_pts = pts.size2();
461
462 data.dataOnEntities[MBTET][0].getN(base).resize(
463 nb_gauss_pts,
464 NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()), false);
465 data.dataOnEntities[MBTET][0].getDiffN(base).resize(
466 nb_gauss_pts,
467 3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()), false);
468
470 data.dataOnEntities[MBTET][0].getOrder(),
471 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
472 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
473 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
474 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
475 nb_gauss_pts, base_polynomials);
476
478}
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
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.

◆ getValueL2BernsteinBezierBase()

MoFEMErrorCode TetPolynomialBase::getValueL2BernsteinBezierBase ( MatrixDouble pts)
private

Definition at line 481 of file TetPolynomialBase.cpp.

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

14 {
15
17 *iface = const_cast<TetPolynomialBase *>(this);
19}
Calculate base functions on tetrahedral.

Member Data Documentation

◆ cTx

EntPolynomialBaseCtx* MoFEM::TetPolynomialBase::cTx
private

Definition at line 31 of file TetPolynomialBase.hpp.

◆ diffN_face_bubble

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

Definition at line 84 of file TetPolynomialBase.hpp.

◆ diffN_face_edge

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

Definition at line 83 of file TetPolynomialBase.hpp.

◆ diffN_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::diffN_volume_bubble
private

Definition at line 87 of file TetPolynomialBase.hpp.

◆ diffN_volume_edge

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

Definition at line 85 of file TetPolynomialBase.hpp.

◆ diffN_volume_face

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

Definition at line 86 of file TetPolynomialBase.hpp.

◆ N_face_bubble

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

Definition at line 78 of file TetPolynomialBase.hpp.

◆ N_face_edge

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

Definition at line 77 of file TetPolynomialBase.hpp.

◆ N_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::N_volume_bubble
private

Definition at line 81 of file TetPolynomialBase.hpp.

◆ N_volume_edge

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

Definition at line 79 of file TetPolynomialBase.hpp.

◆ N_volume_face

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

Definition at line 80 of file TetPolynomialBase.hpp.

◆ senseFaceAlpha

MatrixInt MoFEM::TetPolynomialBase::senseFaceAlpha
private

Definition at line 102 of file TetPolynomialBase.hpp.


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