v0.14.0
Loading...
Searching...
No Matches
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 1348 of file TetPolynomialBase.cpp.

1349 {
1351
1353
1354 int nb_gauss_pts = pts.size2();
1355 if (!nb_gauss_pts)
1357
1358 if (pts.size1() < 3)
1359 SETERRQ(
1360 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1361 "Wrong dimension of pts, should be at least 3 rows with coordinates");
1362
1364 const FieldApproximationBase base = cTx->bAse;
1365 EntitiesFieldData &data = cTx->dAta;
1366 if (cTx->copyNodeBase == LASTBASE) {
1367 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
1368 false);
1370 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1371 &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
1372 } else {
1373 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
1374 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
1375 }
1376 if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
1377 (unsigned int)nb_gauss_pts) {
1378 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1379 "Base functions or nodes has wrong number of integration points "
1380 "for base %s",
1382 }
1383 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
1384 std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
1385 data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1386 }
1387
1388 switch (cTx->sPace) {
1389 case H1:
1390 CHKERR getValueH1(pts);
1391 break;
1392 case HDIV:
1393 CHKERR getValueHdiv(pts);
1394 break;
1395 case HCURL:
1396 CHKERR getValueHcurl(pts);
1397 break;
1398 case L2:
1399 CHKERR getValueL2(pts);
1400 break;
1401 default:
1402 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space");
1403 }
1404
1406}
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:680
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 19 of file TetPolynomialBase.cpp.

19 {
21
22 switch (cTx->bAse) {
26 break;
29 break;
30 default:
31 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
32 }
33
35}
@ 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 37 of file TetPolynomialBase.cpp.

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

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

1329 {
1331
1332 switch (cTx->bAse) {
1336 break;
1339 break;
1340 default:
1341 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1342 }
1343
1345}
@ 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 1079 of file TetPolynomialBase.cpp.

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

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

1062 {
1064
1065 switch (cTx->bAse) {
1068 return getValueHdivAinsworthBase(pts);
1070 return getValueHdivDemkowiczBase(pts);
1071 default:
1072 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1073 }
1074
1076}
#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 655 of file TetPolynomialBase.cpp.

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

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

431 {
433
434 switch (cTx->bAse) {
438 break;
441 break;
442 default:
443 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
444 }
445
447}
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)

◆ getValueL2AinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueL2AinsworthBase ( MatrixDouble pts)
private

Definition at line 449 of file TetPolynomialBase.cpp.

449 {
451
452 EntitiesFieldData &data = cTx->dAta;
453 const FieldApproximationBase base = cTx->bAse;
454 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
455 double *diffL, const int dim) =
457
458 int nb_gauss_pts = pts.size2();
459
460 data.dataOnEntities[MBTET][0].getN(base).resize(
461 nb_gauss_pts,
462 NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()), false);
463 data.dataOnEntities[MBTET][0].getDiffN(base).resize(
464 nb_gauss_pts,
465 3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()), false);
466
468 data.dataOnEntities[MBTET][0].getOrder(),
469 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
470 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
471 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
472 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
473 nb_gauss_pts, base_polynomials);
474
476}
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 479 of file TetPolynomialBase.cpp.

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

12 {
13
15 *iface = const_cast<TetPolynomialBase *>(this);
17}
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: