v0.14.0
Loading...
Searching...
No Matches
FatPrismElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1/** \file FatPrismElementForcesAndSourcesCore.cpp
2
3\brief Implementation of fat prism element
4
5*/
6
7
8
9namespace MoFEM {
10
12 Interface &m_field)
13 : VolumeElementForcesAndSourcesCore(m_field, MBPRISM),
14 dataH1TrianglesOnly(MBPRISM), dataH1TroughThickness(MBPRISM),
15 opHOCoordsAndNormals(hoCoordsAtGaussPtsF3, nOrmals_at_GaussPtF3,
16 tAngent1_at_GaussPtF3, tAngent2_at_GaussPtF3,
17 hoCoordsAtGaussPtsF4, nOrmals_at_GaussPtF4,
18 tAngent1_at_GaussPtF4, tAngent2_at_GaussPtF4) {
20 "Problem with creation data on element");
21}
22
25
26 if (numeredEntFiniteElementPtr->getEntType() != MBPRISM)
28
29 auto get_fe_coordinates = [&]() {
32 int num_nodes;
33 const EntityHandle *conn;
34 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
35 coords.resize(num_nodes * 3, false);
36 CHKERR mField.get_moab().get_coords(conn, num_nodes,
37 &*coords.data().begin());
39 };
40
41 auto calculate_area_of_triangles = [&] {
43 normal.resize(6, false);
46 aRea[0] = cblas_dnrm2(3, &normal[0], 1) * 0.5;
47 aRea[1] = cblas_dnrm2(3, &normal[3], 1) * 0.5;
49 };
50
51 CHKERR get_fe_coordinates();
52 CHKERR calculate_area_of_triangles();
53
57
58 auto get_h1_base_data = [&](auto &dataH1) {
60 CHKERR getEntitySense<MBEDGE>(dataH1);
61 CHKERR getEntitySense<MBTRI>(dataH1);
62 CHKERR getEntitySense<MBQUAD>(dataH1);
63 CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
64 CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
65 CHKERR getEntityDataOrder<MBQUAD>(dataH1, H1);
66 CHKERR getEntityDataOrder<MBPRISM>(dataH1, H1);
67 // Triangles only
68 CHKERR getEntitySense<MBEDGE>(dataH1TrianglesOnly);
69 CHKERR getEntitySense<MBTRI>(dataH1TrianglesOnly);
70 CHKERR getEntityDataOrder<MBEDGE>(dataH1TrianglesOnly, H1);
71 CHKERR getEntityDataOrder<MBTRI>(dataH1TrianglesOnly, H1);
72 // Through thickness
73 CHKERR getEntitySense<MBEDGE>(dataH1TroughThickness);
74 CHKERR getEntityDataOrder<MBEDGE>(dataH1TroughThickness, H1);
76 };
77
78 // H1
79 if ((dataH1.spacesOnEntities[MBEDGE]).test(H1))
80 CHKERR get_h1_base_data(dataH1);
81
82 // Hdiv
83 if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV))
84 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented yet");
85
86 // Hcurl
87 if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL))
88 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented yet");
89
90 // L2
91 if ((dataH1.spacesOnEntities[MBPRISM]).test(L2))
92 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented yet");
93
94 // get approx. on triangles, i.e. faces 3 and 4
95 auto set_gauss_points_on_triangle = [&](int &nb_gauss_pts_on_faces) {
97 int order_triangles_only = 1;
98 int valid_edges[] = {1, 1, 1, 0, 0, 0, 1, 1, 1};
99 for (unsigned int ee = 0; ee < 9; ee++) {
100 if (!valid_edges[ee])
101 continue;
102 order_triangles_only = std::max(
103 order_triangles_only,
104 dataH1TrianglesOnly.dataOnEntities[MBEDGE][ee].getOrder());
105 }
106 for (unsigned int ff = 3; ff <= 4; ff++) {
107 order_triangles_only = std::max(
108 order_triangles_only,
109 dataH1TrianglesOnly.dataOnEntities[MBTRI][ff].getOrder());
110 }
111 for (unsigned int qq = 0; qq < 3; qq++) {
112 order_triangles_only = std::max(
113 order_triangles_only,
114 dataH1TroughThickness.dataOnEntities[MBQUAD][qq].getOrder());
115 }
116 order_triangles_only = std::max(
117 order_triangles_only,
118 dataH1TroughThickness.dataOnEntities[MBPRISM][0].getOrder());
119
120 // integration pts on the triangles surfaces
121 nb_gauss_pts_on_faces = 0;
122 int rule = getRuleTrianglesOnly(order_triangles_only);
123 if (rule >= 0) {
124 if (rule < QUAD_2D_TABLE_SIZE) {
125 if (QUAD_2D_TABLE[rule]->dim != 2) {
126 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
127 }
128 if (QUAD_2D_TABLE[rule]->order < rule) {
129 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
130 "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
131 }
132 nb_gauss_pts_on_faces = QUAD_2D_TABLE[rule]->npoints;
133 gaussPtsTrianglesOnly.resize(3, nb_gauss_pts_on_faces, false);
134 cblas_dcopy(nb_gauss_pts_on_faces, &QUAD_2D_TABLE[rule]->points[1], 3,
135 &gaussPtsTrianglesOnly(0, 0), 1);
136 cblas_dcopy(nb_gauss_pts_on_faces, &QUAD_2D_TABLE[rule]->points[2], 3,
137 &gaussPtsTrianglesOnly(1, 0), 1);
138 cblas_dcopy(nb_gauss_pts_on_faces, QUAD_2D_TABLE[rule]->weights, 1,
139 &gaussPtsTrianglesOnly(2, 0), 1);
140 dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(
141 nb_gauss_pts_on_faces, 3, false);
142 double *shape_ptr = &*dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0]
143 .getN(NOBASE)
144 .data()
145 .begin();
146 cblas_dcopy(3 * nb_gauss_pts_on_faces, QUAD_2D_TABLE[rule]->points, 1,
147 shape_ptr, 1);
148 dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(
149 1, 6, false);
150 std::copy(Tools::diffShapeFunMBTRI.begin(),
153 .getDiffN(NOBASE)
154 .data()
155 .begin());
156 } else
157 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
158 "rule > quadrature order %d < %d", rule, QUAD_2D_TABLE_SIZE);
159
160 } else {
161 CHKERR setGaussPtsTrianglesOnly(order_triangles_only);
162 nb_gauss_pts_on_faces = gaussPtsTrianglesOnly.size2();
163 if (nb_gauss_pts_on_faces == 0)
165 dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(
166 nb_gauss_pts_on_faces, 3, false);
167 if (nb_gauss_pts_on_faces) {
169 .getN(NOBASE)
170 .data()
171 .begin(),
173 &gaussPtsTrianglesOnly(1, 0), nb_gauss_pts_on_faces);
174 dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(
175 1, 6, false);
176 std::copy(Tools::diffShapeFunMBTRI.begin(),
179 .getDiffN(NOBASE)
180 .data()
181 .begin());
182 }
183 }
185 };
186
187 // approx. trough prism thickness
188 auto set_gauss_points_through_thickness =
189 [&](int &nb_gauss_pts_through_thickness) {
191 nb_gauss_pts_through_thickness = 0;
192 int order_thickness = 1;
193 for (unsigned int ee = 3; ee <= 5; ee++) {
194 order_thickness = std::max(
195 order_thickness,
196 dataH1TroughThickness.dataOnEntities[MBEDGE][ee].getOrder());
197 }
198 for (unsigned int qq = 0; qq < 3; qq++) {
199 order_thickness = std::max(
200 order_thickness,
201 dataH1TroughThickness.dataOnEntities[MBQUAD][qq].getOrder());
202 }
203 order_thickness = std::max(
204 order_thickness,
205 dataH1TroughThickness.dataOnEntities[MBPRISM][0].getOrder());
206 // integration points
207 int rule = getRuleThroughThickness(order_thickness);
208 if (rule >= 0) {
209 if (rule < QUAD_1D_TABLE_SIZE) {
210 if (QUAD_1D_TABLE[rule]->dim != 1) {
211 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
212 "wrong dimension");
213 }
214 if (QUAD_1D_TABLE[rule]->order < rule) {
215 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
216 "wrong order %d != %d", QUAD_1D_TABLE[rule]->order,
217 rule);
218 }
219 nb_gauss_pts_through_thickness = QUAD_1D_TABLE[rule]->npoints;
220 gaussPtsThroughThickness.resize(2, nb_gauss_pts_through_thickness,
221 false);
222 cblas_dcopy(nb_gauss_pts_through_thickness,
223 &QUAD_1D_TABLE[rule]->points[1], 2,
224 &gaussPtsThroughThickness(0, 0), 1);
225 cblas_dcopy(nb_gauss_pts_through_thickness,
226 QUAD_1D_TABLE[rule]->weights, 1,
227 &gaussPtsThroughThickness(1, 0), 1);
228 } else {
229 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
230 "rule > quadrature order %d < %d", rule,
232 nb_gauss_pts_through_thickness = 0;
233 }
234 } else {
235 CHKERR setGaussPtsThroughThickness(order_thickness);
236 nb_gauss_pts_through_thickness = gaussPtsThroughThickness.size2();
237 }
239 };
240
241 // Generate integration pts.
242 auto set_gauss_points_in_volume = [&](int nb_gauss_pts_on_faces,
243 int nb_gauss_pts_through_thickness,
244 int &nb_gauss_pts) {
246 nb_gauss_pts = nb_gauss_pts_on_faces * nb_gauss_pts_through_thickness;
247 gaussPts.resize(4, nb_gauss_pts, false);
248 int gg = 0;
249 for (int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
250 for (int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
251 gaussPts(0, gg) = gaussPtsTrianglesOnly(0, ggf);
252 gaussPts(1, gg) = gaussPtsTrianglesOnly(1, ggf);
253 gaussPts(2, gg) = gaussPtsThroughThickness(0, ggt);
254 gaussPts(3, gg) =
256 }
257 }
259 };
260
261 int nb_gauss_pts, nb_gauss_pts_through_thickness, nb_gauss_pts_on_faces;
262 CHKERR set_gauss_points_on_triangle(nb_gauss_pts_on_faces);
263 if (!nb_gauss_pts_on_faces)
265 CHKERR set_gauss_points_through_thickness(nb_gauss_pts_through_thickness);
266 if (!nb_gauss_pts_through_thickness)
268 CHKERR set_gauss_points_in_volume(
269 nb_gauss_pts_on_faces, nb_gauss_pts_through_thickness, nb_gauss_pts);
270
271 auto calc_coordinates_at_triangles = [&]() {
273 coordsAtGaussPtsTrianglesOnly.resize(nb_gauss_pts_on_faces, 6, false);
274 for (int gg = 0; gg < nb_gauss_pts_on_faces; gg++) {
275 for (int dd = 0; dd < 3; dd++) {
277 cblas_ddot(3,
278 &dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(
279 NOBASE)(gg, 0),
280 1, &coords[dd], 3);
282 cblas_ddot(3,
283 &dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(
284 NOBASE)(gg, 0),
285 1, &coords[9 + dd], 3);
286 }
287 }
289 };
290
291 auto calc_vertex_base_on_prism = [&]() {
293 // Calculate "nobase" base functions on prism, this is cartesian product
294 // of base functions on triangles with base functions through thickness
295 // FIXME: This could be effectively implemented with tensors
296 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 6,
297 false);
298 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(nb_gauss_pts, 18,
299 false);
300 for (int dd = 0; dd != 6; dd++) {
301 int gg = 0;
302 for (int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
303 int ddd = dd > 2 ? dd - 3 : dd;
304 double tri_n = dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(
305 NOBASE)(ggf, ddd);
306 double dksi_tri_n =
307 dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(
308 0, 2 * ddd + 0);
309 double deta_tri_n =
310 dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(
311 0, 2 * ddd + 1);
312 for (int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
313 double zeta = gaussPtsThroughThickness(0, ggt);
314 double dzeta, edge_shape;
315 if (dd < 3) {
316 dzeta = diffN_MBEDGE0;
317 edge_shape = N_MBEDGE0(zeta);
318 } else {
319 dzeta = diffN_MBEDGE1;
320 edge_shape = N_MBEDGE1(zeta);
321 }
322 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, dd) =
323 tri_n * edge_shape;
324 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(gg, 3 * dd + 0) =
325 dksi_tri_n * edge_shape;
326 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(gg, 3 * dd + 1) =
327 deta_tri_n * edge_shape;
328 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(gg, 3 * dd + 2) =
329 tri_n * dzeta;
330 }
331 }
332 }
334 };
335
336 auto calc_base_on_prism = [&]() {
338 // Calculate base functions on prism
339 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
340 if (dataH1.bAse.test(b)) {
341 switch (static_cast<FieldApproximationBase>(b)) {
344 if (dataH1.spacesOnEntities[MBVERTEX].test(H1)) {
346 gaussPts,
347 boost::shared_ptr<BaseFunctionCtx>(
352 static_cast<FieldApproximationBase>(b), NOBASE)));
353 }
354 if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
355 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
356 "Not yet implemented");
357 }
358 if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
359 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
360 "Not yet implemented");
361 }
362 if (dataH1.spacesOnEntities[MBTET].test(L2)) {
363 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
364 "Not yet implemented");
365 }
366 break;
367 default:
368 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
369 "Not yet implemented");
370 }
371 }
372 }
374 };
375
376 auto calc_coordinate_on_prism = [&]() {
378 /// Calculate coordinates at integration points
379 coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
380 for (int gg = 0; gg < nb_gauss_pts; gg++) {
381 for (int dd = 0; dd < 3; dd++) {
382 coordsAtGaussPts(gg, dd) = cblas_ddot(
383 6, &dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0), 1,
384 &coords[dd], 3);
385 }
386 }
388 };
389
390 auto calculate_volume = [&]() {
391 auto get_t_w = [&] {
393 };
394
395 auto get_t_coords = [&]() {
397 &coords[0], &coords[1], &coords[2]);
398 };
399
400 FTensor::Index<'i', 3> i;
401 FTensor::Index<'j', 3> j;
403
404 const size_t nb_gauss_pts = gaussPts.size2();
405 auto t_diff_n =
406 dataH1.dataOnEntities[MBVERTEX][0].getFTensor1DiffN<3>(NOBASE);
407
408 double vol = 0;
409 auto t_w = get_t_w();
410 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
411
412 auto t_coords = get_t_coords();
413 t_jac(i, j) = 0;
414 for (size_t n = 0; n != 6; ++n) {
415 t_jac(i, j) += t_coords(i) * t_diff_n(j);
416 ++t_diff_n;
417 ++t_coords;
418 }
419
420 double det;
421 CHKERR determinantTensor3by3(t_jac, det);
422 vol += det * t_w / 2;
423
424 ++t_w;
425 }
426
427 return vol;
428 };
429
430 auto calc_ho_triangle_face_normals = [&]() {
432
433 auto check_field = [&]() {
434 auto field_it =
436 if (field_it != fieldsPtr->get<FieldName_mi_tag>().end())
437 if ((numeredEntFiniteElementPtr->getBitFieldIdData() &
438 (*field_it)->getId())
439 .any())
440 return true;
441 return false;
442 };
443
444 // Check if field meshPositionsFieldName exist
445 if (check_field()) {
446 hoCoordsAtGaussPtsF3.resize(nb_gauss_pts_on_faces, 3, false);
447 nOrmals_at_GaussPtF3.resize(nb_gauss_pts_on_faces, 3, false);
448 tAngent1_at_GaussPtF3.resize(nb_gauss_pts_on_faces, 3, false);
449 tAngent2_at_GaussPtF3.resize(nb_gauss_pts_on_faces, 3, false);
450 hoCoordsAtGaussPtsF4.resize(nb_gauss_pts_on_faces, 3, false);
451 nOrmals_at_GaussPtF4.resize(nb_gauss_pts_on_faces, 3, false);
452 tAngent1_at_GaussPtF4.resize(nb_gauss_pts_on_faces, 3, false);
453 tAngent2_at_GaussPtF4.resize(nb_gauss_pts_on_faces, 3, false);
454 const auto bit_number =
457 CHKERR getEntityFieldData(dataH1TrianglesOnly, bit_number, MBEDGE);
458 CHKERR getEntityFieldData(dataH1TrianglesOnly, bit_number, MBEDGE);
461 } else {
462 hoCoordsAtGaussPtsF3.resize(0, 0, false);
463 nOrmals_at_GaussPtF3.resize(0, 0, false);
464 tAngent1_at_GaussPtF3.resize(0, 0, false);
465 tAngent2_at_GaussPtF3.resize(0, 0, false);
466 hoCoordsAtGaussPtsF4.resize(0, 0, false);
467 nOrmals_at_GaussPtF4.resize(0, 0, false);
468 tAngent1_at_GaussPtF4.resize(0, 0, false);
469 tAngent2_at_GaussPtF4.resize(0, 0, false);
470 }
472 };
473
474 CHKERR calc_coordinates_at_triangles();
475 CHKERR calc_vertex_base_on_prism();
476 CHKERR calc_base_on_prism();
477 CHKERR calc_coordinate_on_prism();
478 vOlume = calculate_volume();
479 CHKERR calc_ho_triangle_face_normals();
480
481 // Iterate over operators
483
485}
486
487MoFEMErrorCode FatPrismElementForcesAndSourcesCore::UserDataOperator::setPtrFE(
490 if (!(ptrFE = dynamic_cast<FatPrismElementForcesAndSourcesCore *>(ptr)))
491 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
492 "User operator and finite element do not work together");
494}
495
496} // namespace MoFEM
FieldApproximationBase
approximation base
Definition definitions.h:58
@ LASTBASE
Definition definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
@ NOBASE
Definition definitions.h:59
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ 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 ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
const int dim
#define diffN_MBEDGE0
derivative of edge shape function
Definition fem_tools.h:107
#define diffN_MBEDGE1
derivative of edge shape function
Definition fem_tools.h:108
#define N_MBEDGE0(x)
edge shape function
Definition fem_tools.h:105
#define N_MBEDGE1(x)
edge shape function
Definition fem_tools.h:106
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition fem_tools.c:182
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
double zeta
Viscous hardening.
Definition plastic.cpp:177
#define QUAD_2D_TABLE_SIZE
Definition quad.h:174
static QUAD *const QUAD_2D_TABLE[]
Definition quad.h:175
static QUAD *const QUAD_1D_TABLE[]
Definition quad.h:164
#define QUAD_1D_TABLE_SIZE
Definition quad.h:163
const Field_multiIndex * fieldsPtr
raw pointer to fields container
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
Deprecated interface functions.
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::bitset< LASTBASE > bAse
bases on element
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
virtual MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
virtual MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
MoFEMErrorCode operator()()
function is run for every finite element
Class used to pass element data to calculate base functions on fat prism.
Calculate base functions on tetrahedralFIXME: Need moab and mofem finite element structure to work (t...
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
MultiIndex Tag for field name.
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
MoFEMErrorCode getEntityFieldData(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MatrixDouble coordsAtGaussPts
coordinated at gauss points
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode getNodesFieldData(EntitiesFieldData &data, const int bit_number) const
Get data on nodes.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition Tools.hpp:104
static MoFEMErrorCode getTriNormal(const double *coords, double *normal, double *d_normal=nullptr)
Get the Tri Normal objectGet triangle normal.
Definition Tools.cpp:355
int npoints
Definition quad.h:29