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
11FatPrismElementForcesAndSourcesCore::FatPrismElementForcesAndSourcesCore(
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) {
19 CHK_THROW_MESSAGE(createDataOnElement(MBPRISM),
20 "Problem with creation data on element");
21}
22
23MoFEMErrorCode FatPrismElementForcesAndSourcesCore::operator()() {
24 MoFEMFunctionBegin;
25
26 if (numeredEntFiniteElementPtr->getEntType() != MBPRISM)
27 MoFEMFunctionReturnHot(0);
28
29 auto get_fe_coordinates = [&]() {
30 MoFEMFunctionBegin;
31 EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
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());
38 MoFEMFunctionReturn(0);
39 };
40
41 auto calculate_area_of_triangles = [&] {
42 MoFEMFunctionBegin;
43 normal.resize(6, false);
44 CHKERR Tools::getTriNormal(&coords[0], &normal[0]);
45 CHKERR Tools::getTriNormal(&coords[9], &normal[3]);
46 aRea[0] = cblas_dnrm2(3, &normal[0], 1) * 0.5;
47 aRea[1] = cblas_dnrm2(3, &normal[3], 1) * 0.5;
48 MoFEMFunctionReturn(0);
49 };
50
51 CHKERR get_fe_coordinates();
52 CHKERR calculate_area_of_triangles();
53
54 CHKERR getSpacesAndBaseOnEntities(dataH1);
55 CHKERR getSpacesAndBaseOnEntities(dataH1TrianglesOnly);
56 CHKERR getSpacesAndBaseOnEntities(dataH1TroughThickness);
57
58 auto get_h1_base_data = [&](auto &dataH1) {
59 MoFEMFunctionBegin;
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);
75 MoFEMFunctionReturn(0);
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) {
96 MoFEMFunctionBegin;
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(),
151 Tools::diffShapeFunMBTRI.end(),
152 &*dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0]
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)
164 MoFEMFunctionReturnHot(0);
165 dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(
166 nb_gauss_pts_on_faces, 3, false);
167 if (nb_gauss_pts_on_faces) {
168 CHKERR ShapeMBTRI(&*dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0]
169 .getN(NOBASE)
170 .data()
171 .begin(),
172 &gaussPtsTrianglesOnly(0, 0),
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(),
177 Tools::diffShapeFunMBTRI.end(),
178 &*dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0]
179 .getDiffN(NOBASE)
180 .data()
181 .begin());
182 }
183 }
184 MoFEMFunctionReturn(0);
185 };
186
187 // approx. trough prism thickness
188 auto set_gauss_points_through_thickness =
189 [&](int &nb_gauss_pts_through_thickness) {
190 MoFEMFunctionBegin;
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,
231 QUAD_1D_TABLE_SIZE);
232 nb_gauss_pts_through_thickness = 0;
233 }
234 } else {
235 CHKERR setGaussPtsThroughThickness(order_thickness);
236 nb_gauss_pts_through_thickness = gaussPtsThroughThickness.size2();
237 }
238 MoFEMFunctionReturn(0);
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) {
245 MoFEMFunctionBegin;
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) =
255 gaussPtsTrianglesOnly(2, ggf) * gaussPtsThroughThickness(1, ggt);
256 }
257 }
258 MoFEMFunctionReturn(0);
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)
264 MoFEMFunctionReturnHot(0);
265 CHKERR set_gauss_points_through_thickness(nb_gauss_pts_through_thickness);
266 if (!nb_gauss_pts_through_thickness)
267 MoFEMFunctionReturnHot(0);
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 = [&]() {
272 MoFEMFunctionBegin;
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++) {
276 coordsAtGaussPtsTrianglesOnly(gg, dd) =
277 cblas_ddot(3,
278 &dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(
279 NOBASE)(gg, 0),
280 1, &coords[dd], 3);
281 coordsAtGaussPtsTrianglesOnly(gg, 3 + dd) =
282 cblas_ddot(3,
283 &dataH1TrianglesOnly.dataOnEntities[MBVERTEX][0].getN(
284 NOBASE)(gg, 0),
285 1, &coords[9 + dd], 3);
286 }
287 }
288 MoFEMFunctionReturn(0);
289 };
290
291 auto calc_vertex_base_on_prism = [&]() {
292 MoFEMFunctionBegin;
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 }
333 MoFEMFunctionReturn(0);
334 };
335
336 auto calc_base_on_prism = [&]() {
337 MoFEMFunctionBegin;
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)) {
342 case AINSWORTH_LEGENDRE_BASE:
343 case AINSWORTH_LOBATTO_BASE:
344 if (dataH1.spacesOnEntities[MBVERTEX].test(H1)) {
345 CHKERR FatPrismPolynomialBase().getValue(
346 gaussPts,
347 boost::shared_ptr<BaseFunctionCtx>(
348 new FatPrismPolynomialBaseCtx(
349 dataH1, dataH1TrianglesOnly, dataH1TroughThickness,
350 gaussPtsTrianglesOnly, gaussPtsThroughThickness,
351 mField.get_moab(), numeredEntFiniteElementPtr.get(), H1,
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 }
373 MoFEMFunctionReturn(0);
374 };
375
376 auto calc_coordinate_on_prism = [&]() {
377 MoFEMFunctionBegin;
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 }
387 MoFEMFunctionReturn(0);
388 };
389
390 auto calculate_volume = [&]() {
391 auto get_t_w = [&] {
392 return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(&gaussPts(2, 0));
393 };
394
395 auto get_t_coords = [&]() {
396 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
397 &coords[0], &coords[1], &coords[2]);
398 };
399
400 FTensor::Index<'i', 3> i;
401 FTensor::Index<'j', 3> j;
402 FTensor::Tensor2<double, 3, 3> t_jac;
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 = [&]() {
431 MoFEMFunctionBegin;
432
433 auto check_field = [&]() {
434 auto field_it =
435 fieldsPtr->get<FieldName_mi_tag>().find(meshPositionsFieldName);
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 =
455 mField.get_field_bit_number(meshPositionsFieldName);
456 CHKERR getNodesFieldData(dataH1TrianglesOnly, bit_number);
457 CHKERR getEntityFieldData(dataH1TrianglesOnly, bit_number, MBEDGE);
458 CHKERR getEntityFieldData(dataH1TrianglesOnly, bit_number, MBEDGE);
459 CHKERR opHOCoordsAndNormals.opRhs(dataH1TrianglesOnly);
460 CHKERR opHOCoordsAndNormals.calculateNormals();
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 }
471 MoFEMFunctionReturn(0);
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
482 CHKERR loopOverOperators();
483
484 MoFEMFunctionReturn(0);
485}
486
487MoFEMErrorCode FatPrismElementForcesAndSourcesCore::UserDataOperator::setPtrFE(
488 ForcesAndSourcesCore *ptr) {
489 MoFEMFunctionBeginHot;
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");
493 MoFEMFunctionReturnHot(0);
494}
495
496} // namespace MoFEM
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
LASTBASE
@ LASTBASE
Definition: definitions.h:69
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
NOBASE
@ NOBASE
Definition: definitions.h:59
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
H1
@ H1
continuous field
Definition: definitions.h:85
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEMFunctionBegin
#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
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
order
constexpr int order
Definition: dg_projection.cpp:18
dim
const int dim
Definition: elec_phys_2D.cpp:12
diffN_MBEDGE0
#define diffN_MBEDGE0
derivative of edge shape function
Definition: fem_tools.h:107
diffN_MBEDGE1
#define diffN_MBEDGE1
derivative of edge shape function
Definition: fem_tools.h:108
N_MBEDGE0
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:105
N_MBEDGE1
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:106
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:182
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::field_it
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
Definition: Projection10NodeCoordsOnField.cpp:123
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1557
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:177
QUAD_2D_TABLE_SIZE
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
QUAD_1D_TABLE
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
QUAD_1D_TABLE_SIZE
#define QUAD_1D_TABLE_SIZE
Definition: quad.h:163
MoFEM::BasicMethod::fieldsPtr
const Field_multiIndex * fieldsPtr
raw pointer to fields container
Definition: LoopMethods.hpp:252
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DataOperator::opRhs
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
Definition: DataOperators.cpp:140
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:50
MoFEM::EntitiesFieldData::bAse
std::bitset< LASTBASE > bAse
bases on element
Definition: EntitiesFieldData.hpp:45
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:56
MoFEM::FEMethod::numeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Definition: LoopMethods.hpp:383
MoFEM::FatPrismElementForcesAndSourcesCore::setGaussPtsTrianglesOnly
virtual MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
Definition: FatPrismElementForcesAndSourcesCore.hpp:39
MoFEM::FatPrismElementForcesAndSourcesCore::setGaussPtsThroughThickness
virtual MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
Definition: FatPrismElementForcesAndSourcesCore.hpp:45
MoFEM::FatPrismElementForcesAndSourcesCore::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: FatPrismElementForcesAndSourcesCore.cpp:23
MoFEM::FatPrismPolynomialBaseCtx
Class used to pass element data to calculate base functions on fat prism.
Definition: FatPrismPolynomialBase.hpp:23
MoFEM::FatPrismPolynomialBase
Calculate base functions on tetrahedralFIXME: Need moab and mofem finite element structure to work (t...
Definition: FatPrismPolynomialBase.hpp:56
MoFEM::FatPrismPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: FatPrismPolynomialBase.cpp:51
MoFEM::FieldName_mi_tag
MultiIndex Tag for field name.
Definition: TagMultiIndices.hpp:67
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
MoFEM::ForcesAndSourcesCore::loopOverOperators
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
Definition: ForcesAndSourcesCore.cpp:1347
MoFEM::ForcesAndSourcesCore::getSpacesAndBaseOnEntities
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
Definition: ForcesAndSourcesCore.cpp:990
MoFEM::ForcesAndSourcesCore::getEntityFieldData
MoFEMErrorCode getEntityFieldData(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
Definition: ForcesAndSourcesCore.cpp:770
MoFEM::ForcesAndSourcesCore::coordsAtGaussPts
MatrixDouble coordsAtGaussPts
coordinated at gauss points
Definition: ForcesAndSourcesCore.hpp:541
MoFEM::ForcesAndSourcesCore::gaussPts
MatrixDouble gaussPts
Matrix of integration points.
Definition: ForcesAndSourcesCore.hpp:109
MoFEM::ForcesAndSourcesCore::getNodesFieldData
MoFEMErrorCode getNodesFieldData(EntitiesFieldData &data, const int bit_number) const
Get data on nodes.
Definition: ForcesAndSourcesCore.cpp:642
MoFEM::ForcesAndSourcesCore::createDataOnElement
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
Definition: ForcesAndSourcesCore.cpp:1305
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
MoFEM::Tools::getTriNormal
static MoFEMErrorCode getTriNormal(const double *coords, double *normal, double *d_normal=nullptr)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:355
QUAD_::npoints
int npoints
Definition: quad.h:29
Generated by Doxygen 1.9.5 and hosted at University of Glasgow