v0.14.0
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 const auto type = numeredEntFiniteElementPtr->getEntType();
29
30 auto get_rule_by_type = [&]() {
31 switch (type) {
32 case MBHEX:
33 return getRule(order_row + 1, order_col + 1, order_data + 1);
34 default:
35 return getRule(order_row, order_col, order_data);
36 }
37 };
38
39 const int rule = get_rule_by_type();
40
41 auto calc_base_for_tet = [&]() {
43 const size_t nb_gauss_pts = gaussPts.size2();
44 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
45 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
46 base.resize(nb_gauss_pts, 4, false);
47 diff_base.resize(nb_gauss_pts, 12, false);
48 CHKERR Tools::shapeFunMBTET(&*base.data().begin(), &gaussPts(0, 0),
49 &gaussPts(1, 0), &gaussPts(2, 0), nb_gauss_pts);
50 double *diff_shape_ptr = &*diff_base.data().begin();
51 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
52 for (int nn = 0; nn != 4; ++nn) {
53 for (int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
54 *diff_shape_ptr = Tools::diffShapeFunMBTET[3 * nn + dd];
55 }
56 }
57 }
59 };
60
61 auto calc_base_for_hex = [&]() {
63 const size_t nb_gauss_pts = gaussPts.size2();
64 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
65 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
66 base.resize(nb_gauss_pts, 8, false);
67 diff_base.resize(nb_gauss_pts, 24, false);
68 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
69 const double ksi = gaussPts(0, gg);
70 const double zeta = gaussPts(1, gg);
71 const double eta = gaussPts(2, gg);
72 base(gg, 0) = N_MBHEX0(ksi, zeta, eta);
73 base(gg, 1) = N_MBHEX1(ksi, zeta, eta);
74 base(gg, 2) = N_MBHEX2(ksi, zeta, eta);
75 base(gg, 3) = N_MBHEX3(ksi, zeta, eta);
76 base(gg, 4) = N_MBHEX4(ksi, zeta, eta);
77 base(gg, 5) = N_MBHEX5(ksi, zeta, eta);
78 base(gg, 6) = N_MBHEX6(ksi, zeta, eta);
79 base(gg, 7) = N_MBHEX7(ksi, zeta, eta);
80 diff_base(gg, 0 * 3 + 0) = diffN_MBHEX0x(zeta, eta);
81 diff_base(gg, 1 * 3 + 0) = diffN_MBHEX1x(zeta, eta);
82 diff_base(gg, 2 * 3 + 0) = diffN_MBHEX2x(zeta, eta);
83 diff_base(gg, 3 * 3 + 0) = diffN_MBHEX3x(zeta, eta);
84 diff_base(gg, 4 * 3 + 0) = diffN_MBHEX4x(zeta, eta);
85 diff_base(gg, 5 * 3 + 0) = diffN_MBHEX5x(zeta, eta);
86 diff_base(gg, 6 * 3 + 0) = diffN_MBHEX6x(zeta, eta);
87 diff_base(gg, 7 * 3 + 0) = diffN_MBHEX7x(zeta, eta);
88 diff_base(gg, 0 * 3 + 1) = diffN_MBHEX0y(ksi, eta);
89 diff_base(gg, 1 * 3 + 1) = diffN_MBHEX1y(ksi, eta);
90 diff_base(gg, 2 * 3 + 1) = diffN_MBHEX2y(ksi, eta);
91 diff_base(gg, 3 * 3 + 1) = diffN_MBHEX3y(ksi, eta);
92 diff_base(gg, 4 * 3 + 1) = diffN_MBHEX4y(ksi, eta);
93 diff_base(gg, 5 * 3 + 1) = diffN_MBHEX5y(ksi, eta);
94 diff_base(gg, 6 * 3 + 1) = diffN_MBHEX6y(ksi, eta);
95 diff_base(gg, 7 * 3 + 1) = diffN_MBHEX7y(ksi, eta);
96 diff_base(gg, 0 * 3 + 2) = diffN_MBHEX0z(ksi, zeta);
97 diff_base(gg, 1 * 3 + 2) = diffN_MBHEX1z(ksi, zeta);
98 diff_base(gg, 2 * 3 + 2) = diffN_MBHEX2z(ksi, zeta);
99 diff_base(gg, 3 * 3 + 2) = diffN_MBHEX3z(ksi, zeta);
100 diff_base(gg, 4 * 3 + 2) = diffN_MBHEX4z(ksi, zeta);
101 diff_base(gg, 5 * 3 + 2) = diffN_MBHEX5z(ksi, zeta);
102 diff_base(gg, 6 * 3 + 2) = diffN_MBHEX6z(ksi, zeta);
103 diff_base(gg, 7 * 3 + 2) = diffN_MBHEX7z(ksi, zeta);
104 }
106 };
107
108 auto set_integration_pts_for_hex = [&]() {
111 rule);
113 };
114
115 auto set_integration_pts_for_tet = [&]() {
117 if (rule < QUAD_3D_TABLE_SIZE) {
118 if (QUAD_3D_TABLE[rule]->dim != 3) {
119 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
120 }
121 if (QUAD_3D_TABLE[rule]->order < rule) {
123 "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
124 }
125 size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
126 gaussPts.resize(4, nb_gauss_pts, false);
127 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
128 &gaussPts(0, 0), 1);
129 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
130 &gaussPts(1, 0), 1);
131 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
132 &gaussPts(2, 0), 1);
133 cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
134 &gaussPts(3, 0), 1);
135
136 // CHKERR calc_base_for_tet();
137
138 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
139 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
140 base.resize(nb_gauss_pts, 4, false);
141 diff_base.resize(nb_gauss_pts, 12, false);
142 double *shape_ptr = &*base.data().begin();
143 cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
144 1);
145 double *diff_shape_ptr = &*diff_base.data().begin();
146 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
147 for (int nn = 0; nn != 4; ++nn) {
148 for (int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
149 *diff_shape_ptr = Tools::diffShapeFunMBTET[3 * nn + dd];
150 }
151 }
152 }
153
154 } else {
155 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
156 "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
157 }
159 };
160
161 if (rule >= 0) {
162 switch (type) {
163 case MBTET:
164 CHKERR set_integration_pts_for_tet();
165 break;
166 case MBHEX:
167 CHKERR set_integration_pts_for_hex();
168 CHKERR calc_base_for_hex();
169 break;
170 default:
171 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
172 "Element type not implemented: %d", type);
173 }
174 } else {
175 // If rule is negative, set user defined integration points
176 CHKERR setGaussPts(order_row, order_col, order_data);
177 const size_t nb_gauss_pts = gaussPts.size2();
178 if (nb_gauss_pts) {
179 switch (type) {
180 case MBTET: {
181 CHKERR calc_base_for_tet();
182 } break;
183 case MBHEX:
184 CHKERR calc_base_for_hex();
185 break;
186 default:
187 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
188 "Element type not implemented: %d", type);
189 }
190 }
191 }
192
194}
195
198 const auto ent = numeredEntFiniteElementPtr->getEnt();
199 const auto type = numeredEntFiniteElementPtr->getEntType();
200 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
201 coords.resize(3 * num_nodes, false);
202 CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
203
204 auto get_tet_t_diff_n = [&]() {
208 };
209
210 auto get_hex_t_diff_n = [&]() {
215 };
216
217 auto get_t_diff_n = [&]() {
218 if (type == MBTET)
219 return get_tet_t_diff_n();
220 return get_hex_t_diff_n();
221 };
222
223 auto t_diff_n = get_t_diff_n();
224
226 &coords[0], &coords[1], &coords[2]);
229 jAc.clear();
230
231 for (size_t n = 0; n != num_nodes; ++n) {
232 tJac(i, j) += t_coords(i) * t_diff_n(j);
233 ++t_coords;
234 ++t_diff_n;
235 }
238
239 if (type == MBTET)
240 vOlume *= G_TET_W1[0] / 6.;
241
242#ifndef NDEBUG
243 if (vOlume <= std::numeric_limits<double>::epsilon() || vOlume != vOlume) {
244 MOFEM_LOG("SELF", Sev::warning) << "Volume is zero " << vOlume;
245 MOFEM_LOG("SELF", Sev::warning) << "Coords " << coords << endl;
246 }
247#endif
248
250}
251
255 // Get coords at Gauss points
257
258 auto &shape_functions = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
259 double *shape_functions_ptr = &*shape_functions.data().begin();
260 const size_t nb_base_functions = shape_functions.size2();
261 const size_t nb_gauss_pts = gaussPts.size2();
262 coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
263 coordsAtGaussPts.clear();
264 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_coords_at_gauss_ptr(
265 &coordsAtGaussPts(0, 0), &coordsAtGaussPts(0, 1),
266 &coordsAtGaussPts(0, 2));
268 shape_functions_ptr);
269 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
271 &coords[0], &coords[1], &coords[2]);
272 for (int bb = 0; bb != nb_base_functions; ++bb) {
273 t_coords_at_gauss_ptr(i) += t_coords(i) * t_shape_functions;
274 ++t_coords;
275 ++t_shape_functions;
276 };
277 ++t_coords_at_gauss_ptr;
278 }
280}
281
285
288 // H1
289 if ((dataH1.spacesOnEntities[MBEDGE]).test(H1)) {
290 CHKERR getEntitySense<MBEDGE>(dataH1);
291 CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
292 }
293 if ((dataH1.spacesOnEntities[MBTRI]).test(H1)) {
294 CHKERR getEntitySense<MBTRI>(dataH1);
295 CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
296 }
297 if ((dataH1.spacesOnEntities[MBTET]).test(H1)) {
298 CHKERR getEntityDataOrder<MBTET>(dataH1, H1);
299 }
300 if ((dataH1.spacesOnEntities[MBQUAD]).test(H1)) {
301 CHKERR getEntitySense<MBQUAD>(dataH1);
302 CHKERR getEntityDataOrder<MBQUAD>(dataH1, H1);
303 }
304 if ((dataH1.spacesOnEntities[MBHEX]).test(H1)) {
305 CHKERR getEntityDataOrder<MBHEX>(dataH1, H1);
306 }
307 // Hcurl
308 if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL)) {
309 CHKERR getEntitySense<MBEDGE>(dataHcurl);
310 CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
311 dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
312 }
313 if ((dataH1.spacesOnEntities[MBTRI]).test(HCURL)) {
315 CHKERR getEntitySense<MBTRI>(dataHcurl);
316 CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
317 dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
318 }
319 if ((dataH1.spacesOnEntities[MBTET]).test(HCURL)) {
320 CHKERR getEntityDataOrder<MBTET>(dataHcurl, HCURL);
321 dataHcurl.spacesOnEntities[MBTET].set(HCURL);
322 }
323 if ((dataH1.spacesOnEntities[MBQUAD]).test(HCURL)) {
326 CHKERR getEntitySense<MBQUAD>(dataHcurl);
327 CHKERR getEntityDataOrder<MBQUAD>(dataHcurl, HCURL);
328 dataHcurl.spacesOnEntities[MBQUAD].set(HCURL);
329 }
330 if ((dataH1.spacesOnEntities[MBHEX]).test(HCURL)) {
331 CHKERR getEntityDataOrder<MBHEX>(dataHcurl, HCURL);
332 dataHcurl.spacesOnEntities[MBHEX].set(HCURL);
333 }
334 // Hdiv
335 if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV)) {
337 CHKERR getEntitySense<MBTRI>(dataHdiv);
338 CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
339 dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
340 }
341 if ((dataH1.spacesOnEntities[MBTET]).test(HDIV)) {
342 CHKERR getEntityDataOrder<MBTET>(dataHdiv, HDIV);
343 dataHdiv.spacesOnEntities[MBTET].set(HDIV);
344 }
345 if ((dataH1.spacesOnEntities[MBQUAD]).test(HDIV)) {
348 CHKERR getEntitySense<MBQUAD>(dataHdiv);
349 CHKERR getEntityDataOrder<MBQUAD>(dataHdiv, HDIV);
350 dataHdiv.spacesOnEntities[MBQUAD].set(HDIV);
351 }
352 if ((dataH1.spacesOnEntities[MBHEX]).test(HDIV)) {
353 CHKERR getEntityDataOrder<MBHEX>(dataHdiv, HDIV);
354 dataHdiv.spacesOnEntities[MBHEX].set(HDIV);
355 }
356 // L2
357 if ((dataH1.spacesOnEntities[MBTET]).test(L2)) {
358 CHKERR getEntityDataOrder<MBTET>(dataL2, L2);
359 dataL2.spacesOnEntities[MBTET].set(L2);
360 }
361 if ((dataH1.spacesOnEntities[MBHEX]).test(L2)) {
362 CHKERR getEntityDataOrder<MBHEX>(dataL2, L2);
363 dataL2.spacesOnEntities[MBHEX].set(L2);
364 }
366}
367
370
371 if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
372 // OK that code is not nice.
373 MatrixDouble new_diff_n;
374 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
376 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
378 if ((data.getDiffN(base).size1() == 4) &&
379 (data.getDiffN(base).size2() == 3)) {
380 const size_t nb_gauss_pts = gaussPts.size2();
381 const size_t nb_base_functions = 4;
382 new_diff_n.resize(nb_gauss_pts, 3 * nb_base_functions, false);
383 double *new_diff_n_ptr = &*new_diff_n.data().begin();
385 new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
386 double *t_diff_n_ptr = &*data.getDiffN(base).data().begin();
387 for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
389 t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
390 for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
391 t_new_diff_n(i) = t_diff_n(i);
392 ++t_new_diff_n;
393 ++t_diff_n;
394 }
395 }
396 data.getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(),
397 false);
398 data.getDiffN(base).swap(new_diff_n);
399 }
400 }
401 }
402
403 if (dataH1.spacesOnEntities[MBVERTEX].test(H1))
405 if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
408 }
409 for (EntityType t = CN::TypeDimensionMap[2].first;
410 t <= CN::TypeDimensionMap[2].second; ++t) {
411 if (dataH1.spacesOnEntities[t].test(HDIV)) {
413 // Fix for tetrahedrons
414 if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
415 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
416 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
417 for (auto t : {MBTRI, MBTET}) {
418 for (auto &d : dataHdiv.dataOnEntities[t]) {
419 d.getN(base) /= 6;
420 d.getDiffN(base) /= 6;
421 }
422 }
423 }
424 }
426 break;
427 }
428 }
429 for (EntityType t = CN::TypeDimensionMap[3].first;
430 t <= CN::TypeDimensionMap[3].second; ++t) {
431 if (dataH1.spacesOnEntities[t].test(L2)) {
433 break;
434 }
435 }
436
438}
439
440MoFEMErrorCode VolumeElementForcesAndSourcesCore::UserDataOperator::setPtrFE(
443 if (!(ptrFE = dynamic_cast<VolumeElementForcesAndSourcesCore *>(ptr)))
444 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
445 "User operator and finite element do not work together");
447}
448
451
452 const auto type = numeredEntFiniteElementPtr->getEntType();
453 if (type != lastEvaluatedElementEntityType) {
454 switch (type) {
455 case MBTET:
457 boost::shared_ptr<BaseFunction>(new TetPolynomialBase());
458 break;
459 case MBHEX:
461 boost::shared_ptr<BaseFunction>(new HexPolynomialBase());
462 break;
463 default:
465 }
468 }
469
473 if (gaussPts.size2() == 0)
478
480
481 // Iterate over operators
483
485}
486
487} // 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
constexpr int order
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
double eta
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
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:1559
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1511
double zeta
Viscous hardening.
Definition: plastic.cpp:177
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:661
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:680
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