v0.13.2
Loading...
Searching...
No Matches
PostProcOnRefMesh.cpp
Go to the documentation of this file.
1/** \file PostProcOnRefMesh.cpp
2 * \brief Postprocess fields on refined mesh made for 10 Node tets
3 *
4 * Create refined mesh, without enforcing continuity between element. Calculate
5 * field values on nodes of that mesh.
6 * \ingroup mofem_fs_post_proc
7 */
8
9
10
11#include <MoFEM.hpp>
12using namespace MoFEM;
13using namespace boost::numeric;
14#include <PostProcOnRefMesh.hpp>
15
17 int side, EntityType type, EntitiesFieldData::EntData &data) {
19
20 if (data.getFieldData().size() == 0)
22
23 if (V) {
24 vAlues.resize(data.getFieldData().size());
25 double *a;
26 CHKERR VecGetArray(V, &a);
27 VectorDofs::iterator it, hi_it;
28 it = data.getFieldDofs().begin();
29 hi_it = data.getFieldDofs().end();
30 auto row_dofs = getFEMethod()->getRowDofsPtr();
31 for (int ii = 0; it != hi_it; it++, ii++) {
32 int local_idx = row_dofs->find((*it)->getLocalUniqueId())
33 ->get()
34 ->getPetscLocalDofIdx();
35 vAlues[ii] = a[local_idx];
36 }
37 CHKERR VecRestoreArray(V, &a);
39 } else {
40 vAluesPtr = &data.getFieldData();
41 }
42
43 auto dof_ptr = data.getFieldDofs()[0];
44 int rank = dof_ptr->getNbOfCoeffs();
45
46 int tag_length = rank;
47 FieldSpace space = dof_ptr->getSpace();
48 switch (space) {
49 case L2:
50 case H1:
51 break;
52 case HCURL:
53 case HDIV:
54 tag_length *= 3;
55 break;
56 default:
57 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
58 "field with that space is not implemented");
59 }
60
61 if (tag_length > 1 && tag_length < 3)
62 tag_length = 3;
63 else if (tag_length > 3 && tag_length < 9)
64 tag_length = 9;
65
66 double def_VAL[tag_length];
67 bzero(def_VAL, tag_length * sizeof(double));
68 Tag th;
69 CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
70 MB_TYPE_DOUBLE, th,
71 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
72
73 // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
74 // no need for L2
75 const void *tags_ptr[mapGaussPts.size()];
76 int nb_gauss_pts = data.getN().size1();
77 if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
78 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
79 "data inconsistency %d!=%d", mapGaussPts.size(), nb_gauss_pts);
80 }
81
82 auto set_float_precision = [](const double x) {
83 if (std::abs(x) < std::numeric_limits<float>::epsilon())
84 return 0.;
85 else
86 return x;
87 };
88
89 switch (space) {
90 case H1:
91 commonData.fieldMap[rowFieldName].resize(nb_gauss_pts);
92 if (type == MBVERTEX) {
93 for (int gg = 0; gg < nb_gauss_pts; gg++) {
94 rval = postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
95 (commonData.fieldMap[rowFieldName])[gg].resize(rank);
96 (commonData.fieldMap[rowFieldName])[gg].clear();
97 }
98 }
99 CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
100 tags_ptr);
101 for (int gg = 0; gg < nb_gauss_pts; gg++) {
102 for (int rr = 0; rr < rank; rr++) {
103 const double val =
104 cblas_ddot((vAluesPtr->size() / rank), &(data.getN(gg)[0]), 1,
105 &((*vAluesPtr)[rr]), rank);
106 (commonData.fieldMap[rowFieldName])[gg][rr] =
107 ((double *)tags_ptr[gg])[rr] += set_float_precision(val);
108 }
109 }
110 break;
111 case L2:
112 commonData.fieldMap[rowFieldName].resize(nb_gauss_pts);
113 for (int gg = 0; gg < nb_gauss_pts; gg++) {
114 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
115 (commonData.fieldMap[rowFieldName])[gg].resize(rank);
116 (commonData.fieldMap[rowFieldName])[gg].clear();
117 }
118 CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
119 tags_ptr);
120 for (int gg = 0; gg < nb_gauss_pts; gg++) {
121 bzero((double *)tags_ptr[gg], sizeof(double) * tag_length);
122 for (int rr = 0; rr < rank; rr++) {
123 const double val =
124 cblas_ddot((vAluesPtr->size() / rank), &(data.getN(gg)[0]), 1,
125 &((*vAluesPtr)[rr]), rank);
126 (commonData.fieldMap[rowFieldName])[gg][rr] =
127 ((double *)tags_ptr[gg])[rr] = set_float_precision(val);
128 }
129 }
130 break;
131 case HCURL:
132 // FIXME: fieldMap not set
133 if (type == MBEDGE && side == 0) {
134 for (int gg = 0; gg < nb_gauss_pts; gg++) {
135 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
136 }
137 }
138 CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
139 tags_ptr);
140 {
142 auto t_n_hcurl = data.getFTensor1N<3>();
143 for (int gg = 0; gg != nb_gauss_pts; gg++) {
144 double *ptr = &((double *)tags_ptr[gg])[0];
145 int ll = 0;
146 for (; ll != static_cast<int>(vAluesPtr->size() / rank); ll++) {
147 FTensor::Tensor1<double *, 3> t_tag_val(ptr, &ptr[1], &ptr[2], 3);
148 for (int rr = 0; rr != rank; rr++) {
149 const double dof_val = (*vAluesPtr)[ll * rank + rr];
150 t_tag_val(i) += set_float_precision(dof_val) * t_n_hcurl(i);
151 ++t_tag_val;
152 }
153 ++t_n_hcurl;
154 }
155 for (; ll != static_cast<int>(data.getN().size2() / 3); ll++) {
156 ++t_n_hcurl;
157 }
158 }
159 }
160 break;
161 case HDIV:
162 // FIXME: fieldMap not set
163 if (moab::CN::Dimension(type) == 2 && side == 0) {
164 for (int gg = 0; gg < nb_gauss_pts; gg++) {
165 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
166 }
167 }
168 CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
169 tags_ptr);
170 {
172 auto t_n_hdiv = data.getFTensor1N<3>();
173 for (int gg = 0; gg != nb_gauss_pts; gg++) {
174 double *ptr = &((double *)tags_ptr[gg])[0];
175 int ll = 0;
176 for (; ll != static_cast<int>(vAluesPtr->size() / rank); ll++) {
177 FTensor::Tensor1<double *, 3> t_tag_val(ptr, &ptr[1], &ptr[2], 3);
178 for (int rr = 0; rr != rank; rr++) {
179 const double dof_val = (*vAluesPtr)[ll * rank + rr];
180 t_tag_val(i) += set_float_precision(dof_val) * t_n_hdiv(i);
181 ++t_tag_val;
182 }
183 ++t_n_hdiv;
184 }
185 for (; ll != static_cast<int>(data.getN().size2() / 3); ll++) {
186 ++t_n_hdiv;
187 }
188 }
189 }
190 break;
191 default:
192 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
193 "field with that space is not implemented");
194 }
195
197}
198
200 int side, EntityType type, EntitiesFieldData::EntData &data) {
202
203 if (data.getFieldData().size() == 0)
205 if (V) {
206 vAlues.resize(data.getFieldData().size());
207 double *a;
208 CHKERR VecGetArray(V, &a);
209 VectorDofs::iterator it, hi_it;
210 it = data.getFieldDofs().begin();
211 hi_it = data.getFieldDofs().end();
212 auto row_dofs = getFEMethod()->getRowDofsPtr();
213 for (int ii = 0; it != hi_it; it++, ii++) {
214 int local_idx = row_dofs->find((*it)->getLocalUniqueId())
215 ->get()
216 ->getPetscLocalDofIdx();
217 vAlues[ii] = a[local_idx];
218 }
219 CHKERR VecRestoreArray(V, &a);
220 vAluesPtr = &vAlues;
221 } else {
222 vAluesPtr = &data.getFieldData();
223 }
224
225 auto dof_ptr = data.getFieldDofs()[0];
226 int rank = dof_ptr->getNbOfCoeffs();
227
228 FieldSpace space = dof_ptr->getSpace();
229 int space_dim = spaceDim;
230 if (space == HCURL || space == HDIV)
231 space_dim = 3;
232
233 int tag_length = rank * space_dim;
234 // FieldSpace space = dof_ptr->getSpace();
235 switch (space) {
236 case L2:
237 case H1:
238 break;
239 case HCURL:
240 case HDIV:
241 tag_length *= space_dim;
242 break;
243 default:
244 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
245 "field with that space is not implemented");
246 }
247
248 double def_VAL[tag_length];
249 bzero(def_VAL, tag_length * sizeof(double));
250 Tag th;
251 CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
252 MB_TYPE_DOUBLE, th,
253 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
254
255 // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
256 // no need for L2
257 const void *tags_ptr[mapGaussPts.size()];
258 int nb_gauss_pts = data.getN().size1();
259 if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
260 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
261 }
262
263 auto set_float_precision = [](const double x) {
264 if (std::abs(x) < std::numeric_limits<float>::epsilon())
265 return 0.;
266 else
267 return x;
268 };
269
270 auto clear_vals = [&]() {
272 for (int gg = 0; gg < nb_gauss_pts; gg++) {
273 (commonData.gradMap[rowFieldName])[gg].resize(rank, space_dim);
274 (commonData.gradMap[rowFieldName])[gg].clear();
275 CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
276 }
278 };
279
280 auto set_vals = [&]() {
282 CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
283 tags_ptr);
284 for (int gg = 0; gg < nb_gauss_pts; gg++) {
285 for (int rr = 0; rr < rank; rr++) {
286 for (int dd = 0; dd < space_dim; dd++) {
287 for (unsigned int dof = 0; dof < (vAluesPtr->size() / rank); dof++) {
288 const double val =
289 data.getDiffN(gg)(dof, dd) * (*vAluesPtr)[rank * dof + rr];
290 (commonData.gradMap[rowFieldName])[gg](rr, dd) += val;
291 }
292 ((double *)tags_ptr[gg])[rank * rr + dd] =
293 (commonData.gradMap[rowFieldName])[gg](rr, dd);
294 }
295 }
296 }
298 };
299
300 commonData.gradMap[rowFieldName].resize(nb_gauss_pts);
301
302 switch (space) {
303 case H1:
304 if (type == MBVERTEX)
305 clear_vals();
306 CHKERR set_vals();
307 break;
308 case L2:
309 clear_vals();
310 CHKERR set_vals();
311 break;
312 default:
313 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
314 "field with that space is not implemented");
315 }
316
318}
319
323}
324
326 int order_triangles_only) {
328 // if(gaussPtsTrianglesOnly.size1()==0 || gaussPtsTrianglesOnly.size2()==0) {
329 // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"post-process mesh not
330 // generated");
331 // }
332
333 int max_ref_level_thick = 0;
334 PetscBool flg = PETSC_TRUE;
335 PetscOptionsGetInt(PETSC_NULL, PETSC_NULL,
336 "-my_max_post_proc_ref_level_prism_thick",
337 &max_ref_level_thick, &flg);
338 int const level_thick = max_ref_level_thick + 1;
339 const EntityHandle *conn;
340 int num_nodes;
341 EntityHandle prism;
342 const int nb_on_triangle = 6;
343 const int nb_through_thickness = 3;
344
345 std::vector<PointsMap3D_multiIndex> pointsMapVector;
346 PointsMap3D_multiIndex pointsMap;
347
348 if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
349 elementsMap.end()) {
350 pointsMapVector = pointsMapVectorMap[numeredEntFiniteElementPtr->getEnt()];
351 } else {
352 pointsMapVector.resize(0);
353 gaussPtsTrianglesOnly.resize(3, nb_on_triangle, false);
354 gaussPtsTrianglesOnly.clear();
355 gaussPtsThroughThickness.resize(2, nb_through_thickness * level_thick,
356 false);
357 gaussPtsThroughThickness.clear();
358
359 MatrixDouble gauss_pts_triangles_only, gauss_pts_through_thickness;
360
361 double inc_thick = 1.0 / (double)level_thick;
362 for (int ll = 0; ll < level_thick; ll++) {
363
364 gauss_pts_triangles_only.resize(3, 3, false);
365 gauss_pts_triangles_only.clear();
366 gauss_pts_triangles_only(0, 0) = 0;
367 gauss_pts_triangles_only(1, 0) = 0;
368 gauss_pts_triangles_only(0, 1) = 1;
369 gauss_pts_triangles_only(1, 1) = 0;
370 gauss_pts_triangles_only(0, 2) = 0;
371 gauss_pts_triangles_only(1, 2) = 1;
372 gauss_pts_through_thickness.resize(2, 2, false);
373 gauss_pts_through_thickness(0, 0) = ll * inc_thick;
374 gauss_pts_through_thickness(0, 1) = (ll + 1) * inc_thick;
375
376 ublas::vector<EntityHandle> prism_conn(6);
377 VectorDouble coords(3);
378 for (int ggf = 0; ggf != 3; ggf++) {
379 double ksi = gauss_pts_triangles_only(0, ggf);
380 double eta = gauss_pts_triangles_only(1, ggf);
381 coords[0] = ksi;
382 coords[1] = eta;
383 for (int ggt = 0; ggt != 2; ggt++) {
384 double zeta = gauss_pts_through_thickness(0, ggt);
385 coords[2] = zeta;
386 int side = ggt * 3 + ggf;
387 CHKERR postProcMesh.create_vertex(&coords[0], prism_conn[side]);
388 }
389 }
390 CHKERR postProcMesh.create_element(MBPRISM, &prism_conn[0], 6, prism);
391
392 elementsMap[numeredEntFiniteElementPtr->getEnt()].push_back(prism);
393 // Range faces;
394 // CHKERR postProcMesh.get_adjacencies(&prism,1,2,true,faces);
395 Range edges;
396 CHKERR postProcMesh.get_adjacencies(&prism, 1, 1, true, edges);
397 EntityHandle meshset;
398 CHKERR postProcMesh.create_meshset(MESHSET_SET, meshset);
399 CHKERR postProcMesh.add_entities(meshset, &prism, 1);
400 // CHKERR postProcMesh.add_entities(meshset,faces);
401 CHKERR postProcMesh.add_entities(meshset, edges);
402 if (tenNodesPostProcTets) {
403 CHKERR postProcMesh.convert_entities(meshset, true, false, false);
404 }
405 CHKERR postProcMesh.delete_entities(&meshset, 1);
406 CHKERR postProcMesh.delete_entities(edges);
407 // CHKERR postProcMesh.delete_entities(faces);
408
409 CHKERR mField.get_moab().get_connectivity(
410 numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
411 MatrixDouble coords_prism_global;
412 coords_prism_global.resize(num_nodes, 3, false);
413 CHKERR mField.get_moab().get_coords(conn, num_nodes,
414 &coords_prism_global(0, 0));
415
416 CHKERR postProcMesh.get_connectivity(prism, conn, num_nodes, false);
417 MatrixDouble coords_prism_local;
418 coords_prism_local.resize(num_nodes, 3, false);
419 CHKERR postProcMesh.get_coords(conn, num_nodes,
420 &coords_prism_local(0, 0));
421
422 const double eps = 1e-6;
423 int ggf = 0, ggt = 0;
424
425 pointsMap.clear();
426
427 for (int nn = 0; nn != num_nodes; nn++) {
428 double ksi = coords_prism_local(nn, 0);
429 double eta = coords_prism_local(nn, 1);
430 double zeta = coords_prism_local(nn, 2);
431 pointsMap.insert(PointsMap3D(ksi * 100, eta * 100, zeta * 100, nn));
432 if (fabs(zeta) < eps) {
433 gaussPtsTrianglesOnly(0, ggf) = ksi;
434 gaussPtsTrianglesOnly(1, ggf) = eta;
435 ggf++;
436 }
437 if (fabs(ksi) < eps && fabs(eta) < eps) {
438 gaussPtsThroughThickness(0, ggt + nb_through_thickness * ll) = zeta;
439 ggt++;
440 }
441 }
442
443 pointsMapVector.push_back(pointsMap);
444
445 for (int nn = 0; nn != num_nodes; nn++) {
446 const double ksi = coords_prism_local(nn, 0);
447 const double eta = coords_prism_local(nn, 1);
448 const double zeta = coords_prism_local(nn, 2);
449 const double n0 = N_MBTRI0(ksi, eta);
450 const double n1 = N_MBTRI1(ksi, eta);
451 const double n2 = N_MBTRI2(ksi, eta);
452 const double e0 = N_MBEDGE0(zeta);
453 const double e1 = N_MBEDGE1(zeta);
454 double coords_global[3];
455 for (int dd = 0; dd != 3; dd++) {
456 coords_global[dd] = n0 * e0 * coords_prism_global(0, dd) +
457 n1 * e0 * coords_prism_global(1, dd) +
458 n2 * e0 * coords_prism_global(2, dd) +
459 n0 * e1 * coords_prism_global(3, dd) +
460 n1 * e1 * coords_prism_global(4, dd) +
461 n2 * e1 * coords_prism_global(5, dd);
462 }
463 CHKERR postProcMesh.set_coords(&conn[nn], 1, coords_global);
464 }
465 }
466
467 pointsMapVectorMap[numeredEntFiniteElementPtr->getEnt()] = pointsMapVector;
468 }
469
470 mapGaussPts.clear();
471 mapGaussPts.resize(nb_through_thickness * nb_on_triangle * level_thick);
472 fill(mapGaussPts.begin(), mapGaussPts.end(), 0);
473 int gg = 0;
474 int ll;
475
476 for (unsigned int ggf = 0; ggf != gaussPtsTrianglesOnly.size2(); ggf++) {
477 const double ksi = gaussPtsTrianglesOnly(0, ggf);
478 const double eta = gaussPtsTrianglesOnly(1, ggf);
479 for (unsigned int ggt = 0; ggt != gaussPtsThroughThickness.size2();
480 ggt++, gg++) {
481 ll = ggt / 3;
482 prism = elementsMap[numeredEntFiniteElementPtr->getEnt()][ll];
483 pointsMap = pointsMapVector[ll];
484 CHKERR postProcMesh.get_connectivity(prism, conn, num_nodes, false);
485
486 const double zeta = gaussPtsThroughThickness(0, ggt);
487 PointsMap3D_multiIndex::iterator it;
488 it = pointsMap.find(boost::make_tuple(ksi * 100, eta * 100, zeta * 100));
489 if (it != pointsMap.end()) {
490 mapGaussPts[gg] = conn[it->nN];
491 }
492 }
493 }
494
495 int g = 0;
496
498}
499
501 int order_thickness) {
503 if (gaussPtsThroughThickness.size1() == 0 ||
504 gaussPtsThroughThickness.size2() == 0) {
505 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
506 "post-process mesh not generated");
507 }
509}
510
514}
515
518 ParallelComm *pcomm_post_proc_mesh =
519 ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
520 if (pcomm_post_proc_mesh == NULL) {
521 // wrapRefMeshComm = boost::make_shared<WrapMPIComm>(mField.get_comm(), false);
522 pcomm_post_proc_mesh =
523 new ParallelComm(&postProcMesh, PETSC_COMM_WORLD/*wrapRefMeshComm->get_comm()*/);
524 }
525
526 Range prims;
527 CHKERR postProcMesh.get_entities_by_type(0, MBPRISM, prims, false);
528 // std::cerr << "total prims size " << prims.size() << std::endl;
529 int rank = mField.get_comm_rank();
530 Range::iterator pit = prims.begin();
531 for (; pit != prims.end(); pit++) {
532 CHKERR postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(), &*pit, 1,
533 &rank);
534 }
535 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
537}
538
541
542 auto generate_tri = [&](auto &gauss_pts) {
544 gauss_pts.resize(3, 3, false);
545 gauss_pts.clear();
546 gauss_pts(0, 0) = 0;
547 gauss_pts(1, 0) = 0;
548 gauss_pts(0, 1) = 1;
549 gauss_pts(1, 1) = 0;
550 gauss_pts(0, 2) = 0;
551 gauss_pts(1, 2) = 1;
552
553 moab::Core core_ref;
554 moab::Interface &moab_ref = core_ref;
555 const EntityHandle *conn;
556 int num_nodes;
557 EntityHandle tri_conn[3];
558 MatrixDouble coords(6, 3);
559 for (int gg = 0; gg != 3; gg++) {
560 coords(gg, 0) = gauss_pts(0, gg);
561 coords(gg, 1) = gauss_pts(1, gg);
562 coords(gg, 2) = 0;
563 CHKERR moab_ref.create_vertex(&coords(gg, 0), tri_conn[gg]);
564 }
565
566 EntityHandle tri;
567 CHKERR moab_ref.create_element(MBTRI, tri_conn, 3, tri);
568
569 if (sixNodePostProcTris) {
570 Range edges;
571 CHKERR moab_ref.get_adjacencies(&tri, 1, 1, true, edges);
572 EntityHandle meshset;
573 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
574 CHKERR moab_ref.add_entities(meshset, &tri, 1);
575 CHKERR moab_ref.add_entities(meshset, edges);
576 CHKERR moab_ref.convert_entities(meshset, true, false, false);
577 CHKERR moab_ref.get_connectivity(tri, conn, num_nodes, false);
578 CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
579 gauss_pts.resize(3, num_nodes, false);
580 gauss_pts.clear();
581 for (int nn = 0; nn < num_nodes; nn++) {
582 gauss_pts(0, nn) = coords(nn, 0);
583 gauss_pts(1, nn) = coords(nn, 1);
584 }
585 } else {
586 CHKERR moab_ref.get_connectivity(tri, conn, num_nodes, false);
587 CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
588 gauss_pts.resize(3, num_nodes, false);
589 gauss_pts.clear();
590 for (int nn = 0; nn < 3; nn++) {
591 gauss_pts(0, nn) = coords(nn, 0);
592 gauss_pts(1, nn) = coords(nn, 1);
593 }
594 }
596 };
597
598 auto generate_quad = [&](auto &gauss_pts) {
600
601 gauss_pts.resize(3, 4, false);
602 gauss_pts.clear();
603
604 gauss_pts(0, 0) = 0;
605 gauss_pts(1, 0) = 0;
606 gauss_pts(0, 1) = 1;
607 gauss_pts(1, 1) = 0;
608 gauss_pts(0, 2) = 1;
609 gauss_pts(1, 2) = 1;
610 gauss_pts(0, 3) = 0;
611 gauss_pts(1, 3) = 1;
612
613 moab::Core core_ref;
614 moab::Interface &moab_ref = core_ref;
615
616 const EntityHandle *conn;
617 int num_nodes;
618 std::array<EntityHandle, 4> quad_conn;
619 MatrixDouble coords(8, 3);
620 for (int gg = 0; gg != 4; gg++) {
621 coords(gg, 0) = gauss_pts(0, gg);
622 coords(gg, 1) = gauss_pts(1, gg);
623 coords(gg, 2) = 0;
624 CHKERR moab_ref.create_vertex(&coords(gg, 0), quad_conn[gg]);
625 }
626
627 EntityHandle quad;
628 CHKERR moab_ref.create_element(MBQUAD, quad_conn.data(), 4, quad);
629
630 if (sixNodePostProcTris) {
631 Range edges;
632 CHKERR moab_ref.get_adjacencies(&quad, 1, 1, true, edges);
633 EntityHandle meshset;
634 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
635 CHKERR moab_ref.add_entities(meshset, &quad, 1);
636 CHKERR moab_ref.add_entities(meshset, edges);
637 CHKERR moab_ref.convert_entities(meshset, true, false, false);
638 CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
639 CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
640 gauss_pts.resize(3, num_nodes, false);
641 gauss_pts.clear();
642 for (int nn = 0; nn != num_nodes; nn++) {
643 gauss_pts(0, nn) = coords(nn, 0);
644 gauss_pts(1, nn) = coords(nn, 1);
645 }
646 } else {
647 CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
648 CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
649 gauss_pts.resize(3, num_nodes, false);
650 gauss_pts.clear();
651 for (int nn = 0; nn != 4; nn++) {
652 gauss_pts(0, nn) = coords(nn, 0);
653 gauss_pts(1, nn) = coords(nn, 1);
654 }
655 }
656
658 };
659
660 CHKERR generate_tri(gaussPtsTri);
661 CHKERR generate_quad(gaussPtsQuad);
662
664}
665
668 if (gaussPtsTri.size1() == 0 && gaussPtsQuad.size1() == 0)
669 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
670 "post-process mesh not generated");
671
672 auto create_tri = [&](auto &gauss_pts) {
673 std::array<EntityHandle, 3> tri_conn;
674 MatrixDouble3by3 coords_tri(3, 3);
675 CHKERR mField.get_moab().get_connectivity(
676 numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
677 CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_tri(0, 0));
678 const int num_nodes_on_ele = gauss_pts.size2();
679
680 for (int gg = 0; gg != num_nodes_on_ele; gg++) {
681
682 const double ksi = gauss_pts(0, gg);
683 const double eta = gauss_pts(1, gg);
684 const double n0 = N_MBTRI0(ksi, eta);
685 const double n1 = N_MBTRI1(ksi, eta);
686 const double n2 = N_MBTRI2(ksi, eta);
687 const double x =
688 n0 * coords_tri(0, 0) + n1 * coords_tri(1, 0) + n2 * coords_tri(2, 0);
689 const double y =
690 n0 * coords_tri(0, 1) + n1 * coords_tri(1, 1) + n2 * coords_tri(2, 1);
691 const double z =
692 n0 * coords_tri(0, 2) + n1 * coords_tri(1, 2) + n2 * coords_tri(2, 2);
693
694 verticesOnTriArrays[0][counterTris * num_nodes_on_ele + gg] = x;
695 verticesOnTriArrays[1][counterTris * num_nodes_on_ele + gg] = y;
696 verticesOnTriArrays[2][counterTris * num_nodes_on_ele + gg] = z;
697 }
698
699 mapGaussPts.resize(num_nodes_on_ele);
700 for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
701 triConn[num_nodes_on_ele * counterTris + nn] =
702 num_nodes_on_ele * counterTris + nn + startingVertTriHandle;
703 mapGaussPts[nn] = triConn[num_nodes_on_ele * counterTris + nn];
704 }
705
706 const auto tri = startingEleTriHandle + counterTris;
707
708 return tri;
709 };
710
711 auto create_quad = [&](auto &gauss_pts) {
712 std::array<EntityHandle, 4> quad_conn;
713 MatrixDouble coords_quad(4, 3);
714 CHKERR mField.get_moab().get_connectivity(
715 numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
716 CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_quad(0, 0));
717 const int num_nodes_on_ele = gauss_pts.size2();
718 for (int gg = 0; gg != num_nodes_on_ele; ++gg) {
719 double ksi = gauss_pts(0, gg);
720 double eta = gauss_pts(1, gg);
721 double n0 = N_MBQUAD0(ksi, eta);
722 double n1 = N_MBQUAD1(ksi, eta);
723 double n2 = N_MBQUAD2(ksi, eta);
724 double n3 = N_MBQUAD3(ksi, eta);
725 double x = n0 * coords_quad(0, 0) + n1 * coords_quad(1, 0) +
726 n2 * coords_quad(2, 0) + n3 * coords_quad(3, 0);
727 double y = n0 * coords_quad(0, 1) + n1 * coords_quad(1, 1) +
728 n2 * coords_quad(2, 1) + n3 * coords_quad(3, 1);
729 double z = n0 * coords_quad(0, 2) + n1 * coords_quad(1, 2) +
730 n2 * coords_quad(2, 2) + n3 * coords_quad(3, 2);
731 verticesOnQuadArrays[0][counterQuads * num_nodes_on_ele + gg] = x;
732 verticesOnQuadArrays[1][counterQuads * num_nodes_on_ele + gg] = y;
733 verticesOnQuadArrays[2][counterQuads * num_nodes_on_ele + gg] = z;
734 }
735
736 mapGaussPts.resize(num_nodes_on_ele);
737
738 for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
739 quadConn[num_nodes_on_ele * counterQuads + nn] =
740 num_nodes_on_ele * counterQuads + nn + startingVertQuadHandle;
741 mapGaussPts[nn] = quadConn[num_nodes_on_ele * counterQuads + nn];
742 }
743
744 const auto quad = startingEleQuadHandle + counterQuads;
745 return quad;
746 };
747
748 EntityHandle tri;
749
750 if (elementsMap.size() == getLoopSize()) {
751 // Note "at" that will trigger error if element is not there.
752 tri = elementsMap.at(numeredEntFiniteElementPtr->getEnt());
753 switch (numeredEntFiniteElementPtr->getEntType()) {
754 case MBTRI:
755 gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(), false);
756 noalias(gaussPts) = gaussPtsTri;
757 break;
758 case MBQUAD:
759 gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(), false);
760 noalias(gaussPts) = gaussPtsQuad;
761 break;
762 default:
763 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
764 "Not implemented element type");
765 }
766
767 // Set values which map nodes with integration points on the prism
768 const EntityHandle *conn;
769 int num_nodes;
770 CHKERR postProcMesh.get_connectivity(tri, conn, num_nodes, false);
771 mapGaussPts.resize(num_nodes);
772 for (int nn = 0; nn != num_nodes; nn++)
773 mapGaussPts[nn] = conn[nn];
774
775 } else {
776 switch (numeredEntFiniteElementPtr->getEntType()) {
777 case MBTRI:
778 gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(), false);
779 noalias(gaussPts) = gaussPtsTri;
780 tri = create_tri(gaussPtsTri);
781 ++counterTris;
782 break;
783 case MBQUAD:
784 gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(), false);
785 noalias(gaussPts) = gaussPtsQuad;
786 tri = create_quad(gaussPtsQuad);
787 ++counterQuads;
788 break;
789 default:
790 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
791 "Not implemented element type");
792 }
793
794 elementsMap[numeredEntFiniteElementPtr->getEnt()] = tri;
795 }
796
798}
799
802
803 ReadUtilIface *iface;
804 CHKERR postProcMesh.query_interface(iface);
805
806 const int number_of_ents_in_the_loop = getLoopSize();
807 if (elementsMap.size() != number_of_ents_in_the_loop) {
808
809 elementsMap.clear();
810 postProcMesh.delete_mesh();
811
812 auto get_number_of_computational_elements = [&]() {
813 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
814
815 auto miit =
816 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
817 boost::make_tuple(this->getFEName(), this->getLoFERank()));
818 auto hi_miit =
819 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
820 boost::make_tuple(this->getFEName(), this->getHiFERank()));
821
822 const int number_of_ents_in_the_loop = this->getLoopSize();
823 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
825 "Wrong size of indicices. Inconsistent size number of iterated "
826 "elements iterated by problem and from range.");
827 }
828
829 std::array<int, MBMAXTYPE> nb_elemms_by_type;
830 std::fill(nb_elemms_by_type.begin(), nb_elemms_by_type.end(), 0);
831
832 for (; miit != hi_miit; ++miit) {
833
834 bool add = true;
835 if (exeTestHook) {
836 numeredEntFiniteElementPtr = *miit;
837 add = exeTestHook(this);
838 }
839
840 if (add) {
841 auto type = (*miit)->getEntType();
842 ++nb_elemms_by_type[type];
843 }
844 }
845
846 return nb_elemms_by_type;
847 };
848
849 auto nb_computational_elements_by_type =
850 get_number_of_computational_elements();
851
852 const int numberOfTriangles = nb_computational_elements_by_type[MBTRI];
853 const int numberOfQuads = nb_computational_elements_by_type[MBQUAD];
854
855 // Here we create vertices using ReadUtilface
856 const int total_number_of_nodes_on_tri =
857 numberOfTriangles * gaussPtsTri.size2();
858 const int total_number_of_nodes_on_quad =
859 numberOfQuads * gaussPtsQuad.size2();
860
861 if (total_number_of_nodes_on_tri) {
862 CHKERR iface->get_node_coords(3, total_number_of_nodes_on_tri, 0,
863 startingVertTriHandle, verticesOnTriArrays);
864 CHKERR iface->get_element_connect(numberOfTriangles, gaussPtsTri.size2(),
865 MBTRI, 0, startingEleTriHandle,
866 triConn);
867 }
868
869 if (total_number_of_nodes_on_quad) {
870 CHKERR iface->get_node_coords(3, total_number_of_nodes_on_quad, 0,
871 startingVertQuadHandle,
872 verticesOnQuadArrays);
873 CHKERR iface->get_element_connect(numberOfQuads, gaussPtsQuad.size2(),
874 MBQUAD, 0, startingEleQuadHandle,
875 quadConn);
876 }
877 }
878
879 counterTris = 0;
880 counterQuads = 0;
881
883}
884
887
888 auto update_elements = [&]() {
890 ReadUtilIface *iface;
891 CHKERR postProcMesh.query_interface(iface);
892
893 if (counterTris)
894 CHKERR iface->update_adjacencies(startingEleTriHandle, counterTris,
895 gaussPtsTri.size2(), triConn);
896 if (counterQuads)
897 CHKERR iface->update_adjacencies(startingEleQuadHandle, counterQuads,
898 gaussPtsQuad.size2(), quadConn);
900 };
901
902 auto resolve_shared = [&]() {
904 ParallelComm *pcomm_post_proc_mesh =
905 ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
906 if (pcomm_post_proc_mesh == NULL) {
907 // wrapRefMeshComm =
908 // boost::make_shared<WrapMPIComm>(mField.get_comm(), false);
909 pcomm_post_proc_mesh = new ParallelComm(
910 &postProcMesh, PETSC_COMM_WORLD /* wrapRefMeshComm->get_comm()*/);
911 }
912
913 Range faces;
914 CHKERR postProcMesh.get_entities_by_dimension(0, 2, faces, false);
915 int rank = mField.get_comm_rank();
916 CHKERR postProcMesh.tag_clear_data(pcomm_post_proc_mesh->part_tag(), faces,
917 &rank);
918 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
920 };
921
922 CHKERR update_elements();
923 CHKERR resolve_shared();
924
926}
927
928template <int RANK>
931 int side, EntityType type, EntitiesFieldData::EntData &data) {
933
934 if (type != MBVERTEX)
936
937 CHKERR loopSideVolumes(feVolName, *sideOpFe);
938
939 // quit if tag is not needed
940 if (!saveOnTag)
942
943 auto &m_field = getPtrFE()->mField;
944 auto field_ptr = m_field.get_field_structure(fieldName);
945 const int rank = field_ptr->getNbOfCoeffs();
946 FieldSpace space = field_ptr->getSpace();
947
948 // auto dof_ptr = data.getFieldDofs()[0];
949 // int rank = dof_ptr->getNbOfCoeffs();
950 // FieldSpace space = dof_ptr->getSpace();
951
952 int full_size = rank * RANK;
953 if (space == HDIV)
954 full_size *= 3;
955 // for paraview
956 int tag_length = full_size > 3 && full_size < 9 ? 9 : full_size;
957 double def_VAL[tag_length];
958 bzero(def_VAL, tag_length * sizeof(double));
959 Tag th;
960 CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
961 MB_TYPE_DOUBLE, th,
962 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
963
964 // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
965 // no need for L2
966 const void *tags_ptr[mapGaussPts.size()];
967 int nb_gauss_pts = data.getN().size1();
968 if (mapGaussPts.size() != (unsigned int)nb_gauss_pts ||
969 nb_gauss_pts != matPtr->size2()) {
970 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
971 }
972 switch (space) {
973 case H1:
974 case L2:
975 case HDIV:
976
977 CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
978 tags_ptr);
979
980 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
981 const double *my_ptr2 = static_cast<const double *>(tags_ptr[gg]);
982 double *my_ptr = const_cast<double *>(my_ptr2);
983 for (int rr = 0; rr != RANK; ++rr) {
984 for (int dd = 0; dd != rank; ++dd) {
985 my_ptr[rank * rr + dd] = (*matPtr)(rank * rr + dd, gg);
986 }
987 }
988 }
989
990 break;
991 default:
992 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
993 "field with that space is not implemented");
994 }
995
997}
998
1000 const std::string field_name, const std::string vol_fe_name,
1001 boost::shared_ptr<MatrixDouble> grad_mat_ptr, bool save_on_tag) {
1003
1004 if (!grad_mat_ptr)
1005 grad_mat_ptr = boost::make_shared<MatrixDouble>();
1006
1007 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> my_side_fe =
1008 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1009
1010 if (mField.check_field("MESH_NODE_POSITIONS"))
1011 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *my_side_fe, true, false, false,
1012 false);
1013
1014 // check number of coefficients
1015 auto field_ptr = mField.get_field_structure(field_name);
1016 const int nb_coefficients = field_ptr->getNbOfCoeffs();
1017
1018 switch (nb_coefficients) {
1019 case 1:
1020 my_side_fe->getOpPtrVector().push_back(
1022 break;
1023 case 3:
1024 my_side_fe->getOpPtrVector().push_back(
1026 break;
1027 case 6:
1028 my_side_fe->getOpPtrVector().push_back(
1030 grad_mat_ptr));
1031 break;
1032 default:
1033 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1034 "field with that number of coefficients is not implemented");
1035 }
1036
1037 FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
1039 postProcMesh, mapGaussPts, field_name, field_name + "_GRAD",
1040 my_side_fe, vol_fe_name, grad_mat_ptr, save_on_tag));
1041
1043}
1044
1046 const std::string field_name, const std::string vol_fe_name,
1047 boost::shared_ptr<MatrixDouble> mat_ptr, bool save_on_tag) {
1049
1050 if (!mat_ptr)
1051 mat_ptr = boost::make_shared<MatrixDouble>();
1052
1053 auto my_side_fe =
1054 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1055
1056 // check number of coefficients
1057 auto field_ptr = mField.get_field_structure(field_name);
1058 const int nb_coefficients = field_ptr->getNbOfCoeffs();
1059 FieldSpace space = field_ptr->getSpace();
1060
1061 switch (space) {
1062 case L2:
1063 case H1:
1064 switch (nb_coefficients) {
1065 case 1:
1066 my_side_fe->getOpPtrVector().push_back(
1068 break;
1069 case 3:
1070 my_side_fe->getOpPtrVector().push_back(
1072 break;
1073 case 6:
1074 my_side_fe->getOpPtrVector().push_back(
1076 break;
1077 case 9:
1078 my_side_fe->getOpPtrVector().push_back(
1080 break;
1081 default:
1082 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1083 "field with that number of coefficients is not implemented");
1084 }
1085 break;
1086 case HDIV:
1087 switch (nb_coefficients) {
1088 case 1:
1089 my_side_fe->getOpPtrVector().push_back(
1091 break;
1092 case 3:
1093 my_side_fe->getOpPtrVector().push_back(
1095 break;
1096 default:
1097 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1098 "field with that number of coefficients is not implemented");
1099 }
1100 break;
1101 default:
1102 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1103 "field with that space is not implemented.");
1104 }
1105
1106 FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
1107 new OpGetFieldValuesOnSkinImpl<1>(postProcMesh, mapGaussPts, field_name,
1108 field_name, my_side_fe, vol_fe_name,
1109 mat_ptr, save_on_tag));
1110
1112}
1113
1116
1117 gaussPts.resize(3, 2, false);
1118 gaussPts.clear();
1119 gaussPts(0, 0) = 0;
1120 gaussPts(1, 0) = 0;
1121 gaussPts(0, 1) = 1;
1122 gaussPts(1, 1) = 0;
1123 mapGaussPts.resize(gaussPts.size2());
1124
1125 moab::Core core_ref;
1126 moab::Interface &moab_ref = core_ref;
1127 const EntityHandle *conn;
1128 int num_nodes;
1129 EntityHandle edge_conn[2];
1130 MatrixDouble coords(2, 3);
1131 for (int gg = 0; gg != 2; gg++) {
1132 coords(gg, 0) = gaussPts(0, gg);
1133 coords(gg, 1) = gaussPts(1, gg);
1134 coords(gg, 2) = 0;
1135 CHKERR moab_ref.create_vertex(&coords(gg, 0), edge_conn[gg]);
1136 }
1137
1138 EntityHandle edge;
1139 CHKERR moab_ref.create_element(MBEDGE, edge_conn, 2, edge);
1141}
1142
1145 ParallelComm *pcomm_post_proc_mesh =
1146 ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
1147 if (pcomm_post_proc_mesh != NULL) {
1148 delete pcomm_post_proc_mesh;
1149 }
1151}
1152
1155 ParallelComm *pcomm_post_proc_mesh =
1156 ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
1157 if (pcomm_post_proc_mesh == NULL) {
1158 // wrapRefMeshComm = boost::make_shared<WrapMPIComm>(mField.get_comm(), false);
1159 pcomm_post_proc_mesh = new ParallelComm(
1160 &postProcMesh, PETSC_COMM_WORLD /*wrapRefMeshComm->get_comm()*/);
1161 }
1162
1163 Range edges;
1164 CHKERR postProcMesh.get_entities_by_type(0, MBEDGE, edges, false);
1165 int rank = mField.get_comm_rank();
1166 auto set_edges_rank = [&](const auto rank, const auto &edges) {
1167 std::vector<EntityHandle> ranks(edges.size(), rank);
1168 CHKERR postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(), edges,
1169 &*ranks.begin());
1170 };
1171 set_edges_rank(rank, edges);
1172
1173 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
1175}
1176
1179 if (gaussPts.size1() == 0) {
1180 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1181 "post-process mesh not generated");
1182 }
1183
1184 const EntityHandle *conn;
1185 int num_nodes;
1186 EntityHandle edge;
1187
1188 if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
1189 elementsMap.end()) {
1190 edge = elementsMap[numeredEntFiniteElementPtr->getEnt()];
1191 } else {
1192 ublas::vector<EntityHandle> edge_conn(2);
1193 MatrixDouble coords_edge(2, 3);
1194 VectorDouble coords(3);
1195 CHKERR mField.get_moab().get_connectivity(
1196 numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
1197 CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_edge(0, 0));
1198 for (int gg = 0; gg != 2; gg++) {
1199 double ksi = gaussPts(0, gg);
1200 // double eta = gaussPts(1, gg);
1201 double n0 = N_MBEDGE0(ksi);
1202 double n1 = N_MBEDGE1(ksi);
1203
1204 double x = n0 * coords_edge(0, 0) + n1 * coords_edge(1, 0);
1205 double y = n0 * coords_edge(0, 1) + n1 * coords_edge(1, 1);
1206 double z = n0 * coords_edge(0, 2) + n1 * coords_edge(1, 2);
1207
1208 coords[0] = x;
1209 coords[1] = y;
1210 coords[2] = z;
1211
1212 CHKERR postProcMesh.create_vertex(&coords[0], edge_conn[gg]);
1213 }
1214 CHKERR postProcMesh.create_element(MBEDGE, &edge_conn[0], 2, edge);
1215 elementsMap[numeredEntFiniteElementPtr->getEnt()] = edge;
1216 }
1217
1218 {
1219 CHKERR postProcMesh.get_connectivity(edge, conn, num_nodes, false);
1220 mapGaussPts.resize(num_nodes);
1221 for (int nn = 0; nn < num_nodes; nn++) {
1222 mapGaussPts[nn] = conn[nn];
1223 }
1224 }
1225
1227}
Post-process fields on refined mesh.
constexpr double a
static const double eps
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ 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 MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#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
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
#define N_MBTRI1(x, y)
triangle shape function
Definition: fem_tools.h:47
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:60
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:105
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:57
#define N_MBTRI0(x, y)
triangle shape function
Definition: fem_tools.h:46
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:59
#define N_MBTRI2(x, y)
triangle shape function
Definition: fem_tools.h:48
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:58
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:106
constexpr double eta
FTensor::Index< 'i', SPACE_DIM > i
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
double zeta
Definition: plastic.cpp:107
constexpr auto field_name
constexpr double g
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
auto getRowDofsPtr() const
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Get vector field for H-div approximation.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Get field gradients at integration pts for symmetric tensorial field rank 2.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
std::map< std::string, std::vector< VectorDouble > > fieldMap
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
const std::string tagName
VectorDouble vAlues
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::vector< EntityHandle > & mapGaussPts
VectorDouble * vAluesPtr
moab::Interface & postProcMesh
Vec V
CommonData & commonData
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode setGaussPts(int order)
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode addFieldValuesPostProcOnSkin(const std::string field_name, const std::string vol_fe_name="dFE", boost::shared_ptr< MatrixDouble > mat_ptr=nullptr, bool save_on_tag=true)
MoFEMErrorCode setGaussPts(int order)
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode addFieldValuesGradientPostProcOnSkin(const std::string field_name, const std::string vol_fe_name="dFE", boost::shared_ptr< MatrixDouble > grad_mat_ptr=nullptr, bool save_on_tag=true)
MoFEMErrorCode postProcess()
function is run at the end of loop
multi_index_container< PointsMap3D, indexed_by< ordered_unique< composite_key< PointsMap3D, member< PointsMap3D, const int, &PointsMap3D::kSi >, member< PointsMap3D, const int, &PointsMap3D::eTa >, member< PointsMap3D, const int, &PointsMap3D::zEta > > > > > PointsMap3D_multiIndex
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode postProcess()
function is run at the end of loop