v0.13.1
Loading...
Searching...
No Matches
VolumeElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1/** \file VolumeElementForcesAndSourcesCore.cpp
2
3\brief Implementation of volume element
4
5*/
6
7namespace MoFEM {
8
10 Interface &m_field, const EntityType type)
11 : ForcesAndSourcesCore(m_field), vOlume(elementMeasure),
12 meshPositionsFieldName("MESH_NODE_POSITIONS"), coords(24), jAc(3, 3),
13 invJac(3, 3), opSetInvJacH1(invJac),
14 opContravariantPiolaTransform(elementMeasure, jAc),
15 opCovariantPiolaTransform(invJac), opSetInvJacHdivAndHcurl(invJac),
16 tJac(&jAc(0, 0), &jAc(0, 1), &jAc(0, 2), &jAc(1, 0), &jAc(1, 1),
17 &jAc(1, 2), &jAc(2, 0), &jAc(2, 1), &jAc(2, 2)),
18 tInvJac(&invJac(0, 0), &invJac(0, 1), &invJac(0, 2), &invJac(1, 0),
19 &invJac(1, 1), &invJac(1, 2), &invJac(2, 0), &invJac(2, 1),
20 &invJac(2, 2)) {}
21
24
25 int order_data = getMaxDataOrder();
26 int order_row = getMaxRowOrder();
27 int order_col = getMaxColOrder();
28 int rule = getRule(order_row, order_col, order_data);
29 const auto type = numeredEntFiniteElementPtr->getEntType();
30
31 auto calc_base_for_tet = [&]() {
33 const size_t nb_gauss_pts = gaussPts.size2();
34 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
35 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
36 base.resize(nb_gauss_pts, 4, false);
37 diff_base.resize(nb_gauss_pts, 12, false);
38 CHKERR Tools::shapeFunMBTET(&*base.data().begin(), &gaussPts(0, 0),
39 &gaussPts(1, 0), &gaussPts(2, 0), nb_gauss_pts);
40 double *diff_shape_ptr = &*diff_base.data().begin();
41 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
42 for (int nn = 0; nn != 4; ++nn) {
43 for (int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
44 *diff_shape_ptr = Tools::diffShapeFunMBTET[3 * nn + dd];
45 }
46 }
47 }
49 };
50
51 auto calc_base_for_hex = [&]() {
53 const size_t nb_gauss_pts = gaussPts.size2();
54 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
55 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
56 base.resize(nb_gauss_pts, 8, false);
57 diff_base.resize(nb_gauss_pts, 24, false);
58 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
59 const double ksi = gaussPts(0, gg);
60 const double zeta = gaussPts(1, gg);
61 const double eta = gaussPts(2, gg);
62 base(gg, 0) = N_MBHEX0(ksi, zeta, eta);
63 base(gg, 1) = N_MBHEX1(ksi, zeta, eta);
64 base(gg, 2) = N_MBHEX2(ksi, zeta, eta);
65 base(gg, 3) = N_MBHEX3(ksi, zeta, eta);
66 base(gg, 4) = N_MBHEX4(ksi, zeta, eta);
67 base(gg, 5) = N_MBHEX5(ksi, zeta, eta);
68 base(gg, 6) = N_MBHEX6(ksi, zeta, eta);
69 base(gg, 7) = N_MBHEX7(ksi, zeta, eta);
70 diff_base(gg, 0 * 3 + 0) = diffN_MBHEX0x(zeta, eta);
71 diff_base(gg, 1 * 3 + 0) = diffN_MBHEX1x(zeta, eta);
72 diff_base(gg, 2 * 3 + 0) = diffN_MBHEX2x(zeta, eta);
73 diff_base(gg, 3 * 3 + 0) = diffN_MBHEX3x(zeta, eta);
74 diff_base(gg, 4 * 3 + 0) = diffN_MBHEX4x(zeta, eta);
75 diff_base(gg, 5 * 3 + 0) = diffN_MBHEX5x(zeta, eta);
76 diff_base(gg, 6 * 3 + 0) = diffN_MBHEX6x(zeta, eta);
77 diff_base(gg, 7 * 3 + 0) = diffN_MBHEX7x(zeta, eta);
78 diff_base(gg, 0 * 3 + 1) = diffN_MBHEX0y(ksi, eta);
79 diff_base(gg, 1 * 3 + 1) = diffN_MBHEX1y(ksi, eta);
80 diff_base(gg, 2 * 3 + 1) = diffN_MBHEX2y(ksi, eta);
81 diff_base(gg, 3 * 3 + 1) = diffN_MBHEX3y(ksi, eta);
82 diff_base(gg, 4 * 3 + 1) = diffN_MBHEX4y(ksi, eta);
83 diff_base(gg, 5 * 3 + 1) = diffN_MBHEX5y(ksi, eta);
84 diff_base(gg, 6 * 3 + 1) = diffN_MBHEX6y(ksi, eta);
85 diff_base(gg, 7 * 3 + 1) = diffN_MBHEX7y(ksi, eta);
86 diff_base(gg, 0 * 3 + 2) = diffN_MBHEX0z(ksi, zeta);
87 diff_base(gg, 1 * 3 + 2) = diffN_MBHEX1z(ksi, zeta);
88 diff_base(gg, 2 * 3 + 2) = diffN_MBHEX2z(ksi, zeta);
89 diff_base(gg, 3 * 3 + 2) = diffN_MBHEX3z(ksi, zeta);
90 diff_base(gg, 4 * 3 + 2) = diffN_MBHEX4z(ksi, zeta);
91 diff_base(gg, 5 * 3 + 2) = diffN_MBHEX5z(ksi, zeta);
92 diff_base(gg, 6 * 3 + 2) = diffN_MBHEX6z(ksi, zeta);
93 diff_base(gg, 7 * 3 + 2) = diffN_MBHEX7z(ksi, zeta);
94 }
96 };
97
98 auto set_integration_pts_for_hex = [&]() {
101 rule);
103 };
104
105 auto set_integration_pts_for_tet = [&]() {
107 if (rule < QUAD_3D_TABLE_SIZE) {
108 if (QUAD_3D_TABLE[rule]->dim != 3) {
109 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
110 }
111 if (QUAD_3D_TABLE[rule]->order < rule) {
113 "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
114 }
115 size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
116 gaussPts.resize(4, nb_gauss_pts, false);
117 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
118 &gaussPts(0, 0), 1);
119 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
120 &gaussPts(1, 0), 1);
121 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
122 &gaussPts(2, 0), 1);
123 cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
124 &gaussPts(3, 0), 1);
125
126 CHKERR calc_base_for_tet();
127
128 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
129 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
130 base.resize(nb_gauss_pts, 4, false);
131 diff_base.resize(nb_gauss_pts, 12, false);
132 double *shape_ptr = &*base.data().begin();
133 cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
134 1);
135 double *diff_shape_ptr = &*diff_base.data().begin();
136 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
137 for (int nn = 0; nn != 4; ++nn) {
138 for (int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
139 *diff_shape_ptr = Tools::diffShapeFunMBTET[3 * nn + dd];
140 }
141 }
142 }
143
144 } else {
145 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
146 "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
147 }
149 };
150
151 if (rule >= 0) {
152 switch (type) {
153 case MBTET:
154 CHKERR set_integration_pts_for_tet();
155 break;
156 case MBHEX:
157 CHKERR set_integration_pts_for_hex();
158 CHKERR calc_base_for_hex();
159 break;
160 default:
161 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
162 "Element type not implemented: %d", type);
163 }
164 } else {
165 // If rule is negative, set user defined integration points
166 CHKERR setGaussPts(order_row, order_col, order_data);
167 const size_t nb_gauss_pts = gaussPts.size2();
168 if (nb_gauss_pts) {
169 switch (type) {
170 case MBTET: {
171 CHKERR calc_base_for_tet();
172 } break;
173 case MBHEX:
174 CHKERR calc_base_for_hex();
175 break;
176 default:
177 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
178 "Element type not implemented: %d", type);
179 }
180 }
181 }
182
184}
185
188 const auto ent = numeredEntFiniteElementPtr->getEnt();
189 const auto type = numeredEntFiniteElementPtr->getEntType();
190 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
191 coords.resize(3 * num_nodes, false);
192 CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
193
194 auto get_tet_t_diff_n = [&]() {
198 };
199
200 auto get_hex_t_diff_n = [&]() {
205 };
206
207 auto get_t_diff_n = [&]() {
208 if (type == MBTET)
209 return get_tet_t_diff_n();
210 return get_hex_t_diff_n();
211 };
212
213 auto t_diff_n = get_t_diff_n();
214
216 &coords[0], &coords[1], &coords[2]);
219 jAc.clear();
220
221 for (size_t n = 0; n != num_nodes; ++n) {
222 tJac(i, j) += t_coords(i) * t_diff_n(j);
223 ++t_coords;
224 ++t_diff_n;
225 }
228
229 if (type == MBTET)
230 vOlume *= G_TET_W1[0] / 6.;
231
233}
234
238 // Get coords at Gauss points
240
241 auto &shape_functions = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
242 double *shape_functions_ptr = &*shape_functions.data().begin();
243 const size_t nb_base_functions = shape_functions.size2();
244 const size_t nb_gauss_pts = gaussPts.size2();
245 coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
246 coordsAtGaussPts.clear();
247 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_coords_at_gauss_ptr(
248 &coordsAtGaussPts(0, 0), &coordsAtGaussPts(0, 1),
249 &coordsAtGaussPts(0, 2));
251 shape_functions_ptr);
252 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
254 &coords[0], &coords[1], &coords[2]);
255 for (int bb = 0; bb != nb_base_functions; ++bb) {
256 t_coords_at_gauss_ptr(i) += t_coords(i) * t_shape_functions;
257 ++t_coords;
258 ++t_shape_functions;
259 };
260 ++t_coords_at_gauss_ptr;
261 }
263}
264
268
271 // H1
272 if ((dataH1.spacesOnEntities[MBEDGE]).test(H1)) {
273 CHKERR getEntitySense<MBEDGE>(dataH1);
274 CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
275 }
276 if ((dataH1.spacesOnEntities[MBTRI]).test(H1)) {
277 CHKERR getEntitySense<MBTRI>(dataH1);
278 CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
279 }
280 if ((dataH1.spacesOnEntities[MBTET]).test(H1)) {
281 CHKERR getEntityDataOrder<MBTET>(dataH1, H1);
282 }
283 if ((dataH1.spacesOnEntities[MBQUAD]).test(H1)) {
284 CHKERR getEntitySense<MBQUAD>(dataH1);
285 CHKERR getEntityDataOrder<MBQUAD>(dataH1, H1);
286 }
287 if ((dataH1.spacesOnEntities[MBHEX]).test(H1)) {
288 CHKERR getEntityDataOrder<MBHEX>(dataH1, H1);
289 }
290 // Hcurl
291 if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL)) {
292 CHKERR getEntitySense<MBEDGE>(dataHcurl);
293 CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
294 dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
295 }
296 if ((dataH1.spacesOnEntities[MBTRI]).test(HCURL)) {
298 CHKERR getEntitySense<MBTRI>(dataHcurl);
299 CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
300 dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
301 }
302 if ((dataH1.spacesOnEntities[MBTET]).test(HCURL)) {
303 CHKERR getEntityDataOrder<MBTET>(dataHcurl, HCURL);
304 dataHcurl.spacesOnEntities[MBTET].set(HCURL);
305 }
306 if ((dataH1.spacesOnEntities[MBQUAD]).test(HCURL)) {
309 CHKERR getEntitySense<MBQUAD>(dataHcurl);
310 CHKERR getEntityDataOrder<MBQUAD>(dataHcurl, HCURL);
311 dataHcurl.spacesOnEntities[MBQUAD].set(HCURL);
312 }
313 if ((dataH1.spacesOnEntities[MBHEX]).test(HCURL)) {
314 CHKERR getEntityDataOrder<MBHEX>(dataHcurl, HCURL);
315 dataHcurl.spacesOnEntities[MBHEX].set(HCURL);
316 }
317 // Hdiv
318 if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV)) {
320 CHKERR getEntitySense<MBTRI>(dataHdiv);
321 CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
322 dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
323 }
324 if ((dataH1.spacesOnEntities[MBTET]).test(HDIV)) {
325 CHKERR getEntityDataOrder<MBTET>(dataHdiv, HDIV);
326 dataHdiv.spacesOnEntities[MBTET].set(HDIV);
327 }
328 if ((dataH1.spacesOnEntities[MBQUAD]).test(HDIV)) {
331 CHKERR getEntitySense<MBQUAD>(dataHdiv);
332 CHKERR getEntityDataOrder<MBQUAD>(dataHdiv, HDIV);
333 dataHdiv.spacesOnEntities[MBQUAD].set(HDIV);
334 }
335 if ((dataH1.spacesOnEntities[MBHEX]).test(HDIV)) {
336 CHKERR getEntityDataOrder<MBHEX>(dataHdiv, HDIV);
337 dataHdiv.spacesOnEntities[MBHEX].set(HDIV);
338 }
339 // L2
340 if ((dataH1.spacesOnEntities[MBTET]).test(L2)) {
341 CHKERR getEntityDataOrder<MBTET>(dataL2, L2);
342 dataL2.spacesOnEntities[MBTET].set(L2);
343 }
344 if ((dataH1.spacesOnEntities[MBHEX]).test(L2)) {
345 CHKERR getEntityDataOrder<MBHEX>(dataL2, L2);
346 dataL2.spacesOnEntities[MBHEX].set(L2);
347 }
349}
350
353
354 if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
355 // OK that code is not nice.
356 MatrixDouble new_diff_n;
357 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
359 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
361 if ((data.getDiffN(base).size1() == 4) &&
362 (data.getDiffN(base).size2() == 3)) {
363 const size_t nb_gauss_pts = gaussPts.size2();
364 const size_t nb_base_functions = 4;
365 new_diff_n.resize(nb_gauss_pts, 3 * nb_base_functions, false);
366 double *new_diff_n_ptr = &*new_diff_n.data().begin();
368 new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
369 double *t_diff_n_ptr = &*data.getDiffN(base).data().begin();
370 for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
372 t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
373 for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
374 t_new_diff_n(i) = t_diff_n(i);
375 ++t_new_diff_n;
376 ++t_diff_n;
377 }
378 }
379 data.getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(),
380 false);
381 data.getDiffN(base).swap(new_diff_n);
382 }
383 }
384 }
385
386 if (dataH1.spacesOnEntities[MBVERTEX].test(H1))
388 if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
391 }
392 for (EntityType t = CN::TypeDimensionMap[2].first;
393 t <= CN::TypeDimensionMap[2].second; ++t) {
394 if (dataH1.spacesOnEntities[t].test(HDIV)) {
396 // Fix for tetrahedrons
397 if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
398 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
399 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
400 for (auto t : {MBTRI, MBTET}) {
401 for (auto &d : dataHdiv.dataOnEntities[t]) {
402 d.getN(base) /= 6;
403 d.getDiffN(base) /= 6;
404 }
405 }
406 }
407 }
409 break;
410 }
411 }
412 for (EntityType t = CN::TypeDimensionMap[3].first;
413 t <= CN::TypeDimensionMap[3].second; ++t) {
414 if (dataH1.spacesOnEntities[t].test(L2)) {
416 break;
417 }
418 }
419
421}
422
423MoFEMErrorCode VolumeElementForcesAndSourcesCore::UserDataOperator::setPtrFE(
426 if (!(ptrFE = dynamic_cast<VolumeElementForcesAndSourcesCore *>(ptr)))
427 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
428 "User operator and finite element do not work together");
430}
431
434
435 const auto type = numeredEntFiniteElementPtr->getEntType();
436 if (type != lastEvaluatedElementEntityType) {
437 switch (type) {
438 case MBTET:
440 boost::shared_ptr<BaseFunction>(new TetPolynomialBase());
441 break;
442 case MBHEX:
444 boost::shared_ptr<BaseFunction>(new HexPolynomialBase());
445 break;
446 default:
448 }
451 }
452
456 if (gaussPts.size2() == 0)
461
463
464 // Iterate over operators
466
468}
469
470} // 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 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: Common.hpp:10
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:1323
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1312
double zeta
Definition: plastic.cpp:104
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