v0.13.1
VolumeElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1/** \file VolumeElementForcesAndSourcesCore.cpp
2
3\brief Implementation of volume element
4
5*/
6
7
8
9namespace MoFEM {
10
12 Interface &m_field, const EntityType type)
13 : ForcesAndSourcesCore(m_field), vOlume(elementMeasure),
14 meshPositionsFieldName("MESH_NODE_POSITIONS"), coords(24), jAc(3, 3),
15 invJac(3, 3), opSetInvJacH1(invJac),
16 opContravariantPiolaTransform(elementMeasure, jAc),
17 opCovariantPiolaTransform(invJac), opSetInvJacHdivAndHcurl(invJac),
18 tJac(&jAc(0, 0), &jAc(0, 1), &jAc(0, 2), &jAc(1, 0), &jAc(1, 1),
19 &jAc(1, 2), &jAc(2, 0), &jAc(2, 1), &jAc(2, 2)),
20 tInvJac(&invJac(0, 0), &invJac(0, 1), &invJac(0, 2), &invJac(1, 0),
21 &invJac(1, 1), &invJac(1, 2), &invJac(2, 0), &invJac(2, 1),
22 &invJac(2, 2)) {}
23
26
27 int order_data = getMaxDataOrder();
28 int order_row = getMaxRowOrder();
29 int order_col = getMaxColOrder();
30 int rule = getRule(order_row, order_col, order_data);
31 const auto type = numeredEntFiniteElementPtr->getEntType();
32
33 auto calc_base_for_hex = [&]() {
35 const size_t nb_gauss_pts = gaussPts.size2();
36 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
37 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
38 base.resize(nb_gauss_pts, 8, false);
39 diff_base.resize(nb_gauss_pts, 24, false);
40 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
41 const double ksi = gaussPts(0, gg);
42 const double zeta = gaussPts(1, gg);
43 const double eta = gaussPts(2, gg);
44 base(gg, 0) = N_MBHEX0(ksi, zeta, eta);
45 base(gg, 1) = N_MBHEX1(ksi, zeta, eta);
46 base(gg, 2) = N_MBHEX2(ksi, zeta, eta);
47 base(gg, 3) = N_MBHEX3(ksi, zeta, eta);
48 base(gg, 4) = N_MBHEX4(ksi, zeta, eta);
49 base(gg, 5) = N_MBHEX5(ksi, zeta, eta);
50 base(gg, 6) = N_MBHEX6(ksi, zeta, eta);
51 base(gg, 7) = N_MBHEX7(ksi, zeta, eta);
52 diff_base(gg, 0 * 3 + 0) = diffN_MBHEX0x(zeta, eta);
53 diff_base(gg, 1 * 3 + 0) = diffN_MBHEX1x(zeta, eta);
54 diff_base(gg, 2 * 3 + 0) = diffN_MBHEX2x(zeta, eta);
55 diff_base(gg, 3 * 3 + 0) = diffN_MBHEX3x(zeta, eta);
56 diff_base(gg, 4 * 3 + 0) = diffN_MBHEX4x(zeta, eta);
57 diff_base(gg, 5 * 3 + 0) = diffN_MBHEX5x(zeta, eta);
58 diff_base(gg, 6 * 3 + 0) = diffN_MBHEX6x(zeta, eta);
59 diff_base(gg, 7 * 3 + 0) = diffN_MBHEX7x(zeta, eta);
60 diff_base(gg, 0 * 3 + 1) = diffN_MBHEX0y(ksi, eta);
61 diff_base(gg, 1 * 3 + 1) = diffN_MBHEX1y(ksi, eta);
62 diff_base(gg, 2 * 3 + 1) = diffN_MBHEX2y(ksi, eta);
63 diff_base(gg, 3 * 3 + 1) = diffN_MBHEX3y(ksi, eta);
64 diff_base(gg, 4 * 3 + 1) = diffN_MBHEX4y(ksi, eta);
65 diff_base(gg, 5 * 3 + 1) = diffN_MBHEX5y(ksi, eta);
66 diff_base(gg, 6 * 3 + 1) = diffN_MBHEX6y(ksi, eta);
67 diff_base(gg, 7 * 3 + 1) = diffN_MBHEX7y(ksi, eta);
68 diff_base(gg, 0 * 3 + 2) = diffN_MBHEX0z(ksi, zeta);
69 diff_base(gg, 1 * 3 + 2) = diffN_MBHEX1z(ksi, zeta);
70 diff_base(gg, 2 * 3 + 2) = diffN_MBHEX2z(ksi, zeta);
71 diff_base(gg, 3 * 3 + 2) = diffN_MBHEX3z(ksi, zeta);
72 diff_base(gg, 4 * 3 + 2) = diffN_MBHEX4z(ksi, zeta);
73 diff_base(gg, 5 * 3 + 2) = diffN_MBHEX5z(ksi, zeta);
74 diff_base(gg, 6 * 3 + 2) = diffN_MBHEX6z(ksi, zeta);
75 diff_base(gg, 7 * 3 + 2) = diffN_MBHEX7z(ksi, zeta);
76 }
78 };
79
80 auto set_integration_pts_for_hex = [&]() {
83 rule);
85 };
86
87 auto set_integration_pts_for_tet = [&]() {
89 if (rule < QUAD_3D_TABLE_SIZE) {
90 if (QUAD_3D_TABLE[rule]->dim != 3) {
91 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "wrong dimension");
92 }
93 if (QUAD_3D_TABLE[rule]->order < rule) {
95 "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
96 }
97 size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
98 gaussPts.resize(4, nb_gauss_pts, false);
99 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
100 &gaussPts(0, 0), 1);
101 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
102 &gaussPts(1, 0), 1);
103 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
104 &gaussPts(2, 0), 1);
105 cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
106 &gaussPts(3, 0), 1);
107
108 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
109 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
110 base.resize(nb_gauss_pts, 4, false);
111 diff_base.resize(nb_gauss_pts, 12, false);
112 double *shape_ptr = &*base.data().begin();
113 cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
114 1);
115 double *diff_shape_ptr = &*diff_base.data().begin();
116 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
117 for (int nn = 0; nn != 4; ++nn) {
118 for (int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
119 *diff_shape_ptr = Tools::diffShapeFunMBTET[3 * nn + dd];
120 }
121 }
122 }
123
124 } else {
126 "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
127 }
129 };
130
131 if (rule >= 0) {
132 switch (type) {
133 case MBTET:
134 CHKERR set_integration_pts_for_tet();
135 break;
136 case MBHEX:
137 CHKERR set_integration_pts_for_hex();
138 CHKERR calc_base_for_hex();
139 break;
140 default:
141 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
142 "Element type not implemented: %d", type);
143 }
144 } else {
145 // If rule is negative, set user defined integration points
146 CHKERR setGaussPts(order_row, order_col, order_data);
147 const size_t nb_gauss_pts = gaussPts.size2();
148 if (nb_gauss_pts) {
149 switch (type) {
150 case MBTET: {
151
152 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
153 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
154 base.resize(nb_gauss_pts, 4, false);
155 diff_base.resize(nb_gauss_pts, 12, false);
156 CHKERR Tools::shapeFunMBTET(&*base.data().begin(), &gaussPts(0, 0),
157 &gaussPts(1, 0), &gaussPts(2, 0),
158 nb_gauss_pts);
159 double *diff_shape_ptr = &*diff_base.data().begin();
160 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
161 for (int nn = 0; nn != 4; ++nn) {
162 for (int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
163 *diff_shape_ptr = Tools::diffShapeFunMBTET[3 * nn + dd];
164 }
165 }
166 }
167
168 } break;
169 case MBHEX:
170 CHKERR calc_base_for_hex();
171 break;
172 default:
173 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
174 "Element type not implemented: %d", type);
175 }
176 }
177 }
178
180}
181
185 const auto ent = numeredEntFiniteElementPtr->getEnt();
186 const auto type = numeredEntFiniteElementPtr->getEntType();
187 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
188 coords.resize(3 * num_nodes, false);
189 CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
190
191 auto get_tet_t_diff_n = [&]() {
195 };
196
197 auto get_hex_t_diff_n = [&]() {
202 };
203
204 auto get_t_diff_n = [&]() {
205 if (type == MBTET)
206 return get_tet_t_diff_n();
207 return get_hex_t_diff_n();
208 };
209
210 auto t_diff_n = get_t_diff_n();
211
213 &coords[0], &coords[1], &coords[2]);
216 jAc.clear();
217
218 for (size_t n = 0; n != num_nodes; ++n) {
219 tJac(i, j) += t_coords(i) * t_diff_n(j);
220 ++t_coords;
221 ++t_diff_n;
222 }
225
226 if (type == MBTET)
227 vOlume *= G_TET_W1[0] / 6.;
228
230}
231
235 // Get coords at Gauss points
237
238 auto &shape_functions = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
239 double *shape_functions_ptr = &*shape_functions.data().begin();
240 const size_t nb_base_functions = shape_functions.size2();
241 const size_t nb_gauss_pts = gaussPts.size2();
242 coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
243 coordsAtGaussPts.clear();
244 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_coords_at_gauss_ptr(
245 &coordsAtGaussPts(0, 0), &coordsAtGaussPts(0, 1),
246 &coordsAtGaussPts(0, 2));
248 shape_functions_ptr);
249 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
251 &coords[0], &coords[1], &coords[2]);
252 for (int bb = 0; bb != nb_base_functions; ++bb) {
253 t_coords_at_gauss_ptr(i) += t_coords(i) * t_shape_functions;
254 ++t_coords;
255 ++t_shape_functions;
256 };
257 ++t_coords_at_gauss_ptr;
258 }
260}
261
265
268 // H1
269 if ((dataH1.spacesOnEntities[MBEDGE]).test(H1)) {
270 CHKERR getEntitySense<MBEDGE>(dataH1);
271 CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
272 }
273 if ((dataH1.spacesOnEntities[MBTRI]).test(H1)) {
274 CHKERR getEntitySense<MBTRI>(dataH1);
275 CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
276 }
277 if ((dataH1.spacesOnEntities[MBTET]).test(H1)) {
278 CHKERR getEntityDataOrder<MBTET>(dataH1, H1);
279 }
280 if ((dataH1.spacesOnEntities[MBQUAD]).test(H1)) {
281 CHKERR getEntitySense<MBQUAD>(dataH1);
282 CHKERR getEntityDataOrder<MBQUAD>(dataH1, H1);
283 }
284 if ((dataH1.spacesOnEntities[MBHEX]).test(H1)) {
285 CHKERR getEntityDataOrder<MBHEX>(dataH1, H1);
286 }
287 // Hcurl
288 if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL)) {
289 CHKERR getEntitySense<MBEDGE>(dataHcurl);
290 CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
291 dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
292 }
293 if ((dataH1.spacesOnEntities[MBTRI]).test(HCURL)) {
295 CHKERR getEntitySense<MBTRI>(dataHcurl);
296 CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
297 dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
298 }
299 if ((dataH1.spacesOnEntities[MBTET]).test(HCURL)) {
300 CHKERR getEntityDataOrder<MBTET>(dataHcurl, HCURL);
301 dataHcurl.spacesOnEntities[MBTET].set(HCURL);
302 }
303 if ((dataH1.spacesOnEntities[MBQUAD]).test(HCURL)) {
306 CHKERR getEntitySense<MBQUAD>(dataHcurl);
307 CHKERR getEntityDataOrder<MBQUAD>(dataHcurl, HCURL);
308 dataHcurl.spacesOnEntities[MBQUAD].set(HCURL);
309 }
310 if ((dataH1.spacesOnEntities[MBHEX]).test(HCURL)) {
311 CHKERR getEntityDataOrder<MBHEX>(dataHcurl, HCURL);
312 dataHcurl.spacesOnEntities[MBHEX].set(HCURL);
313 }
314 // Hdiv
315 if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV)) {
317 CHKERR getEntitySense<MBTRI>(dataHdiv);
318 CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
319 dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
320 }
321 if ((dataH1.spacesOnEntities[MBTET]).test(HDIV)) {
322 CHKERR getEntityDataOrder<MBTET>(dataHdiv, HDIV);
323 dataHdiv.spacesOnEntities[MBTET].set(HDIV);
324 }
325 if ((dataH1.spacesOnEntities[MBQUAD]).test(HDIV)) {
328 CHKERR getEntitySense<MBQUAD>(dataHdiv);
329 CHKERR getEntityDataOrder<MBQUAD>(dataHdiv, HDIV);
330 dataHdiv.spacesOnEntities[MBQUAD].set(HDIV);
331 }
332 if ((dataH1.spacesOnEntities[MBHEX]).test(HDIV)) {
333 CHKERR getEntityDataOrder<MBHEX>(dataHdiv, HDIV);
334 dataHdiv.spacesOnEntities[MBHEX].set(HDIV);
335 }
336 // L2
337 if ((dataH1.spacesOnEntities[MBTET]).test(L2)) {
338 CHKERR getEntityDataOrder<MBTET>(dataL2, L2);
339 dataL2.spacesOnEntities[MBTET].set(L2);
340 }
341 if ((dataH1.spacesOnEntities[MBHEX]).test(L2)) {
342 CHKERR getEntityDataOrder<MBHEX>(dataL2, L2);
343 dataL2.spacesOnEntities[MBHEX].set(L2);
344 }
346}
347
350
351 if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
352 // OK that code is not nice.
353 MatrixDouble new_diff_n;
354 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
358 dataH1.dataOnEntities[MBVERTEX][0];
359 if ((data.getDiffN(base).size1() == 4) &&
360 (data.getDiffN(base).size2() == 3)) {
361 const size_t nb_gauss_pts = gaussPts.size2();
362 const size_t nb_base_functions = 4;
363 new_diff_n.resize(nb_gauss_pts, 3 * nb_base_functions, false);
364 double *new_diff_n_ptr = &*new_diff_n.data().begin();
366 new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
367 double *t_diff_n_ptr = &*data.getDiffN(base).data().begin();
368 for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
370 t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
371 for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
372 t_new_diff_n(i) = t_diff_n(i);
373 ++t_new_diff_n;
374 ++t_diff_n;
375 }
376 }
377 data.getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(),
378 false);
379 data.getDiffN(base).swap(new_diff_n);
380 }
381 }
382 }
383
384 if (dataH1.spacesOnEntities[MBVERTEX].test(H1))
386 if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
389 }
390 for (EntityType t = CN::TypeDimensionMap[2].first;
391 t <= CN::TypeDimensionMap[2].second; ++t) {
392 if (dataH1.spacesOnEntities[t].test(HDIV)) {
394 // Fix for tetrahedrons
395 if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
396 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
398 for (auto t : {MBTRI, MBTET}) {
399 for (auto &d : dataHdiv.dataOnEntities[t]) {
400 d.getN(base) /= 6;
401 d.getDiffN(base) /= 6;
402 }
403 }
404 }
405 }
407 break;
408 }
409 }
410 for (EntityType t = CN::TypeDimensionMap[3].first;
411 t <= CN::TypeDimensionMap[3].second; ++t) {
412 if (dataH1.spacesOnEntities[t].test(L2)) {
414 break;
415 }
416 }
417
419}
420
425 if (!(ptrFE = dynamic_cast<VolumeElementForcesAndSourcesCore *>(ptr)))
426 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
427 "User operator and finite element do not work together");
429}
430
433
434 const auto type = numeredEntFiniteElementPtr->getEntType();
436 switch (type) {
437 case MBTET:
439 boost::shared_ptr<BaseFunction>(new TetPolynomialBase());
440 break;
441 case MBHEX:
443 boost::shared_ptr<BaseFunction>(new HexPolynomialBase());
444 break;
445 default:
447 }
450 }
451
455 if (gaussPts.size2() == 0)
460
462
463 // Iterate over operators
465
467}
468
469} // namespace MoFEM
MatrixDouble invJac
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ NOBASE
Definition: definitions.h:59
#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
#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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const int dim
#define diffN_MBHEX3z(x, y)
Definition: fem_tools.h:98
#define diffN_MBHEX7z(x, y)
Definition: fem_tools.h:102
#define diffN_MBHEX7y(x, z)
Definition: fem_tools.h:94
#define diffN_MBHEX6x(y, z)
Definition: fem_tools.h:85
#define N_MBHEX7(x, y, z)
Definition: fem_tools.h:78
#define N_MBHEX3(x, y, z)
Definition: fem_tools.h:74
#define N_MBHEX5(x, y, z)
Definition: fem_tools.h:76
#define diffN_MBHEX1y(x, z)
Definition: fem_tools.h:88
#define diffN_MBHEX3y(x, z)
Definition: fem_tools.h:90
static const double G_TET_W1[]
Definition: fem_tools.h:1115
#define diffN_MBHEX4y(x, z)
Definition: fem_tools.h:91
#define diffN_MBHEX5x(y, z)
Definition: fem_tools.h:84
#define diffN_MBHEX1z(x, y)
Definition: fem_tools.h:96
#define diffN_MBHEX4z(x, y)
Definition: fem_tools.h:99
#define diffN_MBHEX6y(x, z)
Definition: fem_tools.h:93
#define diffN_MBHEX6z(x, y)
Definition: fem_tools.h:101
#define N_MBHEX4(x, y, z)
Definition: fem_tools.h:75
#define N_MBHEX0(x, y, z)
Definition: fem_tools.h:71
#define N_MBHEX6(x, y, z)
Definition: fem_tools.h:77
#define diffN_MBHEX5z(x, y)
Definition: fem_tools.h:100
#define N_MBHEX2(x, y, z)
Definition: fem_tools.h:73
#define diffN_MBHEX0x(y, z)
Definition: fem_tools.h:79
#define diffN_MBHEX1x(y, z)
Definition: fem_tools.h:80
#define diffN_MBHEX2z(x, y)
Definition: fem_tools.h:97
#define diffN_MBHEX5y(x, z)
Definition: fem_tools.h:92
#define diffN_MBHEX4x(y, z)
Definition: fem_tools.h:83
#define N_MBHEX1(x, y, z)
Definition: fem_tools.h:72
#define diffN_MBHEX3x(y, z)
Definition: fem_tools.h:82
#define diffN_MBHEX2y(x, z)
Definition: fem_tools.h:89
#define diffN_MBHEX0y(x, z)
Definition: fem_tools.h:87
#define diffN_MBHEX0z(x, y)
Definition: fem_tools.h:95
#define diffN_MBHEX2x(y, z)
Definition: fem_tools.h:81
#define diffN_MBHEX7x(y, z)
Definition: fem_tools.h:86
FTensor::Index< 'n', SPACE_DIM > n
constexpr double eta
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
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
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1204
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1193
constexpr double t
plate stiffness
Definition: plate.cpp:59
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
MatrixInt facesNodesOrder
order of face nodes on element
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
int getMaxRowOrder() const
Get max order of approximation for field in rows.
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
MoFEMErrorCode getFaceNodes(EntitiesFieldData &data) const
Get nodes on faces.
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
MatrixDouble coordsAtGaussPts
coordinated at gauss points
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
int getMaxColOrder() const
Get max order of approximation for field in columns.
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
Calculate base functions on tetrahedral.
Calculate base functions on tetrahedral.
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForHex(MatrixDouble &pts, const int edge0, const int edge1, const int edge2)
Definition: Tools.cpp:629
static constexpr std::array< double, 24 > diffShapeFunMBHEXAtCenter
Definition: Tools.hpp:218
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
virtual MoFEMErrorCode transformBaseFunctions()
Transform base functions based on geometric element Jacobian.
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
MoFEMErrorCode operator()()
function is run for every finite element
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
virtual MoFEMErrorCode calculateVolumeAndJacobian()
Calculate element volume and Jacobian.
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
OpSetContravariantPiolaTransform opContravariantPiolaTransform
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
int order
Definition: quad.h:28
int npoints
Definition: quad.h:29