v0.14.0
Loading...
Searching...
No Matches
PostProcOnRefMesh.hpp
Go to the documentation of this file.
1/** \file PostProcOnRefMesh.hpp
2 * \brief Post-process fields on refined mesh
3 *
4 * Create refined mesh, without enforcing continuity between element. Calculate
5 * field values on nodes of that mesh.
6 *
7 * \ingroup mofem_fs_post_proc
8 */
9
10
11
12#ifndef __POSTPROC_ON_REF_MESH_HPP
13#define __POSTPROC_ON_REF_MESH_HPP
14
15/** \brief Set of operators and data structures used for post-processing
16
17This set of functions works that for given problem a new MoAB instance is crated
18only used for post-processing. For each post-processed element in the problem
19an element entity in post-processing mesh is created. Post-processed elements do
20not share nodes and any others entities between them, such that discontinuities
21between element could be shown.
22
23Post-processed entities could be represented by ho-elements, for example 10 node
24tetrahedrons. Moreover each element could be refined such that higher order
25polynomials could be well represented.
26
27* \ingroup mofem_fs_post_proc
28
29*/
31
32 struct CommonData {
33 std::map<std::string, std::vector<VectorDouble>> fieldMap;
34 std::map<std::string, std::vector<MatrixDouble>> gradMap;
35 };
36
39 };
40
41 /**
42 * \brief operator to post-process (save gradients on refined post-processing
43 * mesh) field gradient \ingroup mofem_fs_post_proc
44 *
45 * \todo Implamentation of setting values to fieldMap for Hcurl and Hdiv not
46 * implemented
47 *
48 */
51
52 moab::Interface &postProcMesh;
53 std::vector<EntityHandle> &mapGaussPts;
55 const std::string tagName;
56 Vec V;
57
58 OpGetFieldValues(moab::Interface &post_proc_mesh,
59 std::vector<EntityHandle> &map_gauss_pts,
60 const std::string field_name, const std::string tag_name,
61 CommonData &common_data, Vec v = PETSC_NULL)
64 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
65 commonData(common_data), tagName(tag_name), V(v) {}
66
67 VectorDouble vAlues;
68 VectorDouble *vAluesPtr;
69
70 MoFEMErrorCode doWork(int side, EntityType type,
71 EntitiesFieldData::EntData &data);
72 };
73
74 /**
75 * \brief operator to post-process (save gradients on refined post-processing
76 * mesh) field gradient \ingroup mofem_fs_post_proc
77 *
78 * \todo Implementation for Hdiv and Hcurl to be implemented
79 *
80 */
83
84 moab::Interface &postProcMesh;
85 std::vector<EntityHandle> &mapGaussPts;
87 const std::string tagName;
88 Vec V;
90
91 OpGetFieldGradientValues(moab::Interface &post_proc_mesh,
92 std::vector<EntityHandle> &map_gauss_pts,
93 const std::string field_name,
94 const std::string tag_name,
95 CommonData &common_data, Vec v = PETSC_NULL,
96 int space_dim = 3)
99 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
100 commonData(common_data), tagName(tag_name), V(v),
101 spaceDim(space_dim) {}
102
103 VectorDouble vAlues;
104 VectorDouble *vAluesPtr;
105
106 MoFEMErrorCode doWork(int side, EntityType type,
107 EntitiesFieldData::EntData &data);
108 };
109};
110
111/**
112 * \brief Generic post-processing class
113 *
114 * Generate refined mesh and save data on vertices
115 *
116 * \ingroup mofem_fs_post_proc
117 */
118template <class ELEMENT> struct PostProcTemplateOnRefineMesh : public ELEMENT {
119
120 moab::Core coreMesh;
121
122 moab::Interface &postProcMesh;
123 boost::shared_ptr<WrapMPIComm> wrapRefMeshComm;
124
125 std::vector<EntityHandle> mapGaussPts;
126
128 : ELEMENT(m_field), postProcMesh(coreMesh) {}
129
131 ParallelComm *pcomm_post_proc_mesh =
132 ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
133 if (pcomm_post_proc_mesh != NULL)
134 delete pcomm_post_proc_mesh;
135 }
136
138 THROW_MESSAGE("not implemented");
139 }
140
141 /** \brief Add operator to post-process L2, H1, Hdiv, Hcurl field value
142
143 \param field_name
144 \param v If vector is given, values from vector are used to set tags on mesh
145
146 Note:
147 Name of the tag to store values on post-process mesh is the same as field
148 name
149
150 * \ingroup mofem_fs_post_proc
151
152 */
153 MoFEMErrorCode addFieldValuesPostProc(const std::string field_name,
154 Vec v = PETSC_NULL) {
156 ELEMENT::getOpPtrVector().push_back(
159 getCommonData(), v));
161 }
162
163 /** \brief Add operator to post-process L2 or H1 field value
164
165 \param field_name
166 \param tag_name to store results on post-process mesh
167 \param v If vector is given, values from vector are used to set tags on mesh
168
169 * \ingroup mofem_fs_post_proc
170
171 */
172 MoFEMErrorCode addFieldValuesPostProc(const std::string field_name,
173 const std::string tag_name,
174 Vec v = PETSC_NULL) {
176 ELEMENT::getOpPtrVector().push_back(
178 field_name, tag_name,
179 getCommonData(), v));
181 }
182
183 /** \brief Add operator to post-process L2 or H1 field gradient
184
185 \param field_name
186 \param v If vector is given, values from vector are used to set tags on mesh
187
188 Note:
189 Name of the tag to store values on post-process mesh is the same as field
190 name
191
192 * \ingroup mofem_fs_post_proc
193
194 */
195 MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name,
196 Vec v = PETSC_NULL) {
198 ELEMENT::getOpPtrVector().push_back(
201 getCommonData(), v));
203 }
204
205 /** \brief Add operator to post-process L2 or H1 field gradient
206
207 \param field_name
208 \param tag_name to store results on post-process mesh
209 \param v If vector is given, values from vector are used to set tags on mesh
210
211 * \ingroup mofem_fs_post_proc
212
213 */
214 MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name,
215 const std::string tag_name,
216 Vec v = PETSC_NULL) {
218 ELEMENT::getOpPtrVector().push_back(
221 v));
223 }
224
225 /** \brief Add operator to post-process L2 or H1 field gradient
226
227 \param field_name
228 \param space_dim the dimension of the problem
229 \param v If vector is given, values from vector are used to set tags on mesh
230
231 * \ingroup mofem_fs_post_proc
232
233*/
234 MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name,
235 int space_dim,
236 Vec v = PETSC_NULL) {
238 ELEMENT::getOpPtrVector().push_back(
241 getCommonData(), v, space_dim));
243 }
244
245 /**
246 * \brief wrote results in (MOAB) format, use "file_name.h5m"
247 * @param file_name file name (should always end with .h5m)
248 * @return error code
249
250 * \ingroup mofem_fs_post_proc
251
252 */
253 MoFEMErrorCode writeFile(const std::string file_name,
254 const char *file_type = "MOAB",
255 const char *file_options = "PARALLEL=WRITE_PART") {
257 ParallelComm *pcomm_post_proc_mesh =
258 ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
259 if (pcomm_post_proc_mesh == NULL)
260 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
261 "ParallelComm not allocated");
262 CHKERR postProcMesh.write_file(file_name.c_str(), file_type, file_options);
264 }
265};
266
267template <class VOLUME_ELEMENT>
269 : public PostProcTemplateOnRefineMesh<VOLUME_ELEMENT> {
270
272
275
277 bool ten_nodes_post_proc_tets = true,
278 int nb_ref_levels = -1)
279 : PostProcTemplateOnRefineMesh<VOLUME_ELEMENT>(m_field),
280 tenNodesPostProcTets(ten_nodes_post_proc_tets),
282
283 // Gauss pts set on refined mesh
284 int getRule(int order) { return -1; };
285
288
290 return commonData;
291 }
292
293 /** \brief Generate reference mesh on single element
294
295 Each element is subdivided on smaller elements, i.e. a reference mesh on
296 single element is created. Nodes of such reference mesh are used as
297 integration points at which field values are calculated and to
298 each node a "moab" tag is attached to store those values.
299
300 */
303
304 auto get_nb_of_ref_levels_from_options = [this] {
305 if (nbOfRefLevels == -1) {
306 int max_level = 0;
307 PetscBool flg = PETSC_TRUE;
308 PetscOptionsGetInt(PETSC_NULL, PETSC_NULL,
309 "-my_max_post_proc_ref_level", &max_level, &flg);
310 return max_level;
311 } else {
312 return nbOfRefLevels;
313 }
314 return 0;
315 };
316
317 auto generate_for_hex = [&]() {
319
320 moab::Core core_ref;
321 moab::Interface &moab_ref = core_ref;
322
323 auto create_reference_element = [&moab_ref]() {
325 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0,
326 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1};
327 EntityHandle nodes[8];
328 for (int nn = 0; nn < 8; nn++)
329 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
330 EntityHandle hex;
331 CHKERR moab_ref.create_element(MBHEX, nodes, 8, hex);
333 };
334
335 auto add_ho_nodes = [&]() {
337 Range hexes;
338 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, true);
339 EntityHandle meshset;
340 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
341 CHKERR moab_ref.add_entities(meshset, hexes);
342 CHKERR moab_ref.convert_entities(meshset, true, true, true);
343 CHKERR moab_ref.delete_entities(&meshset, 1);
345 };
346
347 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
349 Range hexes;
350 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, true);
351 Range hexes_nodes;
352 CHKERR moab_ref.get_connectivity(hexes, hexes_nodes, false);
353 auto &gauss_pts = levelGaussPtsOnRefMeshHexes[0];
354 gauss_pts.resize(hexes_nodes.size(), 4, false);
355 size_t gg = 0;
356 for (auto node : hexes_nodes) {
357 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
358 little_map[node] = gg;
359 ++gg;
360 }
361 gauss_pts = trans(gauss_pts);
363 };
364
365 auto set_ref_hexes = [&](std::map<EntityHandle, int> &little_map) {
367 Range hexes;
368 CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, true);
369 size_t hh = 0;
370 auto &ref_hexes = levelRefHexes[0];
371 for (auto hex : hexes) {
372 const EntityHandle *conn;
373 int num_nodes;
374 CHKERR moab_ref.get_connectivity(hex, conn, num_nodes, false);
375 if (ref_hexes.size2() != num_nodes) {
376 ref_hexes.resize(hexes.size(), num_nodes);
377 }
378 for (int nn = 0; nn != num_nodes; ++nn) {
379 ref_hexes(hh, nn) = little_map[conn[nn]];
380 }
381 ++hh;
382 }
384 };
385
386 auto set_shape_functions = [&]() {
388 auto &gauss_pts = levelGaussPtsOnRefMeshHexes[0];
389 auto &shape_functions = levelShapeFunctionsHexes[0];
390 const auto nb_gauss_pts = gauss_pts.size2();
391 shape_functions.resize(nb_gauss_pts, 8);
392 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
393 const double ksi = gauss_pts(0, gg);
394 const double zeta = gauss_pts(1, gg);
395 const double eta = gauss_pts(2, gg);
396 shape_functions(gg, 0) = N_MBHEX0(ksi, zeta, eta);
397 shape_functions(gg, 1) = N_MBHEX1(ksi, zeta, eta);
398 shape_functions(gg, 2) = N_MBHEX2(ksi, zeta, eta);
399 shape_functions(gg, 3) = N_MBHEX3(ksi, zeta, eta);
400 shape_functions(gg, 4) = N_MBHEX4(ksi, zeta, eta);
401 shape_functions(gg, 5) = N_MBHEX5(ksi, zeta, eta);
402 shape_functions(gg, 6) = N_MBHEX6(ksi, zeta, eta);
403 shape_functions(gg, 7) = N_MBHEX7(ksi, zeta, eta);
404 }
406 };
407
408 levelRefHexes.resize(1);
410 levelShapeFunctionsHexes.resize(1);
411
412 CHKERR create_reference_element();
414 CHKERR add_ho_nodes();
415 std::map<EntityHandle, int> little_map;
416 CHKERR set_gauss_pts(little_map);
417 CHKERR set_ref_hexes(little_map);
418 CHKERR set_shape_functions();
419
421 };
422
423 auto generate_for_tet = [&]() {
425
426 const int max_level = get_nb_of_ref_levels_from_options();
427
428 moab::Core core_ref;
429 moab::Interface &moab_ref = core_ref;
430
431 auto create_reference_element = [&moab_ref]() {
433 constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
434 EntityHandle nodes[4];
435 for (int nn = 0; nn < 4; nn++) {
436 CHKERR
437 moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
438 }
439 EntityHandle tet;
440 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
442 };
443
444 MoFEM::CoreTmp<-1> m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
445 MoFEM::Interface &m_field_ref = m_core_ref;
446
447 auto refine_ref_tetrahedron = [this, &m_field_ref, max_level]() {
449 // seed ref mofem database by setting bit ref level to reference
450 // tetrahedron
451 CHKERR
452 m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
453 0, 3, BitRefLevel().set(0));
454 for (int ll = 0; ll != max_level; ++ll) {
455 MOFEM_TAG_AND_LOG_C("WORLD", Sev::verbose, "PostProc",
456 "Refine Level %d", ll);
457 Range edges;
458 CHKERR m_field_ref.getInterface<BitRefManager>()
459 ->getEntitiesByTypeAndRefLevel(
460 BitRefLevel().set(ll), BitRefLevel().set(), MBEDGE, edges);
461 Range tets;
462 CHKERR m_field_ref.getInterface<BitRefManager>()
463 ->getEntitiesByTypeAndRefLevel(
464 BitRefLevel().set(ll), BitRefLevel(ll).set(), MBTET, tets);
465 // refine mesh
466 MeshRefinement *m_ref;
467 CHKERR m_field_ref.getInterface(m_ref);
468 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
469 edges, BitRefLevel().set(ll + 1));
470 CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
471 }
473 };
474
475 auto get_ref_gauss_pts_and_shape_functions = [this, max_level, &moab_ref,
476 &m_field_ref]() {
478 for (int ll = 0; ll != max_level + 1; ++ll) {
479 Range tets;
480 CHKERR
481 m_field_ref.getInterface<BitRefManager>()
482 ->getEntitiesByTypeAndRefLevel(
483 BitRefLevel().set(ll), BitRefLevel().set(ll), MBTET, tets);
485 EntityHandle meshset;
486 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
487 CHKERR moab_ref.add_entities(meshset, tets);
488 CHKERR moab_ref.convert_entities(meshset, true, false, false);
489 CHKERR moab_ref.delete_entities(&meshset, 1);
490 }
491 Range elem_nodes;
492 CHKERR moab_ref.get_connectivity(tets, elem_nodes, false);
493
494 auto &gauss_pts = levelGaussPtsOnRefMeshTets[ll];
495 gauss_pts.resize(elem_nodes.size(), 4, false);
496 std::map<EntityHandle, int> little_map;
497 Range::iterator nit = elem_nodes.begin();
498 for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
499 CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
500 little_map[*nit] = gg;
501 }
502 gauss_pts = trans(gauss_pts);
503
504 auto &ref_tets = levelRefTets[ll];
505 Range::iterator tit = tets.begin();
506 for (int tt = 0; tit != tets.end(); ++tit, ++tt) {
507 const EntityHandle *conn;
508 int num_nodes;
509 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
510 if (tt == 0) {
511 ref_tets.resize(tets.size(), num_nodes);
512 }
513 for (int nn = 0; nn != num_nodes; ++nn) {
514 ref_tets(tt, nn) = little_map[conn[nn]];
515 }
516 }
517
518 auto &shape_functions = levelShapeFunctionsTets[ll];
519 shape_functions.resize(elem_nodes.size(), 4);
520 CHKERR ShapeMBTET(&*shape_functions.data().begin(), &gauss_pts(0, 0),
521 &gauss_pts(1, 0), &gauss_pts(2, 0),
522 elem_nodes.size());
523 }
525 };
526
527 levelRefTets.resize(max_level + 1);
528 levelGaussPtsOnRefMeshTets.resize(max_level + 1);
529 levelShapeFunctionsTets.resize(max_level + 1);
530
531 CHKERR create_reference_element();
532 CHKERR refine_ref_tetrahedron();
533 CHKERR get_ref_gauss_pts_and_shape_functions();
534
536 };
537
538 CHKERR generate_for_hex();
539 CHKERR generate_for_tet();
540
542 }
543
544 size_t getMaxLevel() const {
545 auto get_element_max_dofs_order = [&]() {
546 int max_order = 0;
547 auto dofs_vec = this->getDataVectorDofsPtr();
548 for (auto &dof : *dofs_vec) {
549 const int dof_order = dof->getDofOrder();
550 max_order = (max_order < dof_order) ? dof_order : max_order;
551 };
552 return max_order;
553 };
554 const auto dof_max_order = get_element_max_dofs_order();
555 return (dof_max_order > 0) ? (dof_max_order - 1) / 2 : 0;
556 };
557
558 /** \brief Set integration points
559
560 If reference mesh is generated on single elements. This function maps
561 reference coordinates into physical coordinates and create element
562 on post-processing mesh.
563
564 */
565 MoFEMErrorCode setGaussPts(int order) {
567
568 auto type = type_from_handle(this->getFEEntityHandle());
569
570 auto set_gauss_pts = [&](auto &level_gauss_pts_on_ref_mesh, auto &level_ref,
571 auto &level_shape_functions,
572
573 auto start_vert_handle, auto start_ele_handle,
574 auto &verts_array, auto &conn, auto &ver_count,
575 auto &ele_count
576
577 ) {
579
580 auto level =
581 std::min(getMaxLevel(), level_gauss_pts_on_ref_mesh.size() - 1);
582
583 auto &level_ref_gauss_pts = level_gauss_pts_on_ref_mesh[level];
584 auto &level_ref_ele = level_ref[level];
585 auto &shape_functions = level_shape_functions[level];
586 T::gaussPts.resize(level_ref_gauss_pts.size1(),
587 level_ref_gauss_pts.size2(), false);
588 noalias(T::gaussPts) = level_ref_gauss_pts;
589
590 EntityHandle fe_ent = T::numeredEntFiniteElementPtr->getEnt();
591 {
592 const EntityHandle *conn;
593 int num_nodes;
594 T::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes, true);
595 T::coords.resize(3 * num_nodes, false);
596 CHKERR T::mField.get_moab().get_coords(conn, num_nodes, &T::coords[0]);
597 }
598
599 const int num_nodes = level_ref_gauss_pts.size2();
600 T::mapGaussPts.resize(level_ref_gauss_pts.size2());
601
604 &*shape_functions.data().begin());
606 &verts_array[0][ver_count], &verts_array[1][ver_count],
607 &verts_array[2][ver_count]);
608 for (int gg = 0; gg != num_nodes; ++gg, ++ver_count) {
609
610 T::mapGaussPts[gg] = start_vert_handle + ver_count;
611
612 auto set_float_precision = [](const double x) {
613 if (std::abs(x) < std::numeric_limits<float>::epsilon())
614 return 0.;
615 else
616 return x;
617 };
618
619 t_coords(i) = 0;
620 auto t_ele_coords = getFTensor1FromArray<3, 3>(T::coords);
621 for (int nn = 0; nn != CN::VerticesPerEntity(type); ++nn) {
622 t_coords(i) += t_n * t_ele_coords(i);
623 ++t_ele_coords;
624 ++t_n;
625 }
626
627 for (auto ii : {0, 1, 2})
628 t_coords(ii) = set_float_precision(t_coords(ii));
629
630 ++t_coords;
631 }
632
633 Tag th;
634 int def_in_the_loop = -1;
635 CHKERR T::postProcMesh.tag_get_handle(
636 "NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER, th,
637 MB_TAG_CREAT | MB_TAG_SPARSE, &def_in_the_loop);
638
639 commonData.tEts.clear();
640 const int num_el = level_ref_ele.size1();
641 const int num_nodes_on_ele = level_ref_ele.size2();
642 auto start_e = start_ele_handle + ele_count;
643 commonData.tEts = Range(start_e, start_e + num_el - 1);
644 for (auto tt = 0; tt != level_ref_ele.size1(); ++tt, ++ele_count) {
645 for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
646 conn[num_nodes_on_ele * ele_count + nn] =
647 T::mapGaussPts[level_ref_ele(tt, nn)];
648 }
649 }
650
651 const int n_in_the_loop = T::nInTheLoop;
652 CHKERR T::postProcMesh.tag_clear_data(th, commonData.tEts,
653 &n_in_the_loop);
654
656 };
657
658 switch (type) {
659 case MBTET:
660 return set_gauss_pts(levelGaussPtsOnRefMeshTets, levelRefTets,
662
665
666 );
667 case MBHEX:
668 return set_gauss_pts(levelGaussPtsOnRefMeshHexes, levelRefHexes,
670
673
674 );
675 default:
676 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
677 "Element type not implemented");
678 }
679
681 }
682
683 /** \brief Clear operators list
684
685 Clear operators list, user can use the same mesh instance to post-process
686 different problem or the same problem with different set of post-processed
687 fields.
688
689 */
690 MoFEMErrorCode clearOperators() {
692 T::getOpPtrVector().clear();
694 }
695
696 MoFEMErrorCode preProcess() {
698 moab::Interface &moab = T::coreMesh;
699 ParallelComm *pcomm_post_proc_mesh =
700 ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
701 if (pcomm_post_proc_mesh != NULL)
702 delete pcomm_post_proc_mesh;
703
704 CHKERR T::postProcMesh.delete_mesh();
705
706 auto alloc_vertices_and_elements_on_post_proc_mesh = [&]() {
708
709 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
710
711 auto miit =
712 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
713 boost::make_tuple(this->getFEName(), this->getLoFERank()));
714 auto hi_miit =
715 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
716 boost::make_tuple(this->getFEName(), this->getHiFERank()));
717
718 const int number_of_ents_in_the_loop = this->getLoopSize();
719 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
720 SETERRQ(this->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
721 "Wrong size of indicices. Inconsistent size number of iterated "
722 "elements iterated by problem and from range.");
723 }
724
725 int nb_tet_vertices = 0;
726 int nb_tets = 0;
727 int nb_hex_vertices = 0;
728 int nb_hexes = 0;
729
730 for (; miit != hi_miit; ++miit) {
731 auto type = (*miit)->getEntType();
732
733 // Set pointer to element. So that getDataVectorDofsPtr in getMaxLevel
734 // can work
735 this->numeredEntFiniteElementPtr = *miit;
736
737 bool add = true;
738 if (this->exeTestHook) {
739 add = this->exeTestHook(this);
740 }
741
742 if (add) {
743
744 auto level = getMaxLevel();
745
746 switch (type) {
747 case MBTET:
748 level = std::min(level, levelGaussPtsOnRefMeshTets.size() - 1);
749 nb_tet_vertices += levelGaussPtsOnRefMeshTets[level].size2();
750 nb_tets += levelRefTets[level].size1();
751 break;
752 case MBHEX:
753 level = std::min(level, levelGaussPtsOnRefMeshHexes.size() - 1);
754 nb_hex_vertices += levelGaussPtsOnRefMeshHexes[level].size2();
755 nb_hexes += levelRefHexes[level].size1();
756 break;
757 default:
758 SETERRQ(this->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
759 "Element type not implemented");
760 break;
761 }
762 }
763 }
764
765 ReadUtilIface *iface;
766 CHKERR this->postProcMesh.query_interface(iface);
767
768 if (nb_tets) {
769 CHKERR iface->get_node_coords(
770 3, nb_tet_vertices, 0, startingVertTetHandle, verticesOnTetArrays);
771 CHKERR iface->get_element_connect(nb_tets, levelRefTets[0].size2(),
772 MBTET, 0, startingEleTetHandle,
773 tetConn);
774 }
775
776 if (nb_hexes) {
777 CHKERR iface->get_node_coords(
778 3, nb_hex_vertices, 0, startingVertHexHandle, verticesOnHexArrays);
779 CHKERR iface->get_element_connect(nb_hexes, levelRefHexes[0].size2(),
780 MBHEX, 0, startingEleHexHandle,
781 hexConn);
782 }
783
784 countTet = 0;
785 countVertTet = 0;
786 countHex = 0;
787 countVertHex = 0;
788
790 };
791
792 CHKERR alloc_vertices_and_elements_on_post_proc_mesh();
793
795 }
796
797 MoFEMErrorCode postProcess() {
799
800 auto update_elements = [&]() {
802 ReadUtilIface *iface;
803 CHKERR this->postProcMesh.query_interface(iface);
804
805 if (countTet) {
806 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
807 << "Update tets " << countTet;
808
809 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
810 << "Nb nodes on tets " << levelRefTets[0].size2();
811
812 CHKERR iface->update_adjacencies(startingEleTetHandle, countTet,
813 levelRefTets[0].size2(), tetConn);
814 }
815 if (countHex) {
816 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
817 << "Update hexes " << countHex;
818 CHKERR iface->update_adjacencies(startingEleHexHandle, countHex,
819 levelRefHexes[0].size2(), hexConn);
820 }
822 };
823
824 auto resolve_shared_ents = [&]() {
826 ParallelComm *pcomm_post_proc_mesh =
827 ParallelComm::get_pcomm(&(T::postProcMesh), MYPCOMM_INDEX);
828 if (pcomm_post_proc_mesh == NULL) {
829 // T::wrapRefMeshComm =
830 // boost::make_shared<WrapMPIComm>(T::mField.get_comm(), false);
831 pcomm_post_proc_mesh = new ParallelComm(
833 PETSC_COMM_WORLD /*(T::wrapRefMeshComm)->get_comm()*/);
834 }
835
836 Range edges;
837 CHKERR T::postProcMesh.get_entities_by_type(0, MBEDGE, edges, false);
838 CHKERR T::postProcMesh.delete_entities(edges);
839 Range faces;
840 CHKERR T::postProcMesh.get_entities_by_dimension(0, 2, faces, false);
841 CHKERR T::postProcMesh.delete_entities(faces);
842
843 Range ents;
844 CHKERR T::postProcMesh.get_entities_by_dimension(0, 3, ents, false);
845
846 int rank = T::mField.get_comm_rank();
847 CHKERR T::postProcMesh.tag_clear_data(pcomm_post_proc_mesh->part_tag(),
848 ents, &rank);
849
850 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
851
853 };
854
855 CHKERR resolve_shared_ents();
856 CHKERR update_elements();
857
859 }
860
861 /** \brief Add operator to post-process Hdiv field
862 */
863 MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name) {
865 T::getOpPtrVector().push_back(
868 }
869
870 struct OpHdivFunctions : public VOLUME_ELEMENT::UserDataOperator {
871
872 moab::Interface &postProcMesh;
873 std::vector<EntityHandle> &mapGaussPts;
874
875 OpHdivFunctions(moab::Interface &post_proc_mesh,
876 std::vector<EntityHandle> &map_gauss_pts,
877 const std::string field_name)
878 : VOLUME_ELEMENT::UserDataOperator(field_name,
879 T::UserDataOperator::OPCOL),
880 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
881
882 MoFEMErrorCode doWork(int side, EntityType type,
883 EntitiesFieldData::EntData &data) {
885
886 if (data.getIndices().size() == 0)
888
889 std::vector<Tag> th;
890 th.resize(data.getFieldData().size());
891
892 double def_VAL[9] = {0, 0, 0};
893
894 switch (type) {
895 case MBTRI:
896 for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
897 std::ostringstream ss;
898 ss << "HDIV_FACE_" << side << "_" << dd;
899 CHKERR postProcMesh.tag_get_handle(
900 ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
901 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
902 }
903 break;
904 case MBTET:
905 for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
906 std::ostringstream ss;
907 ss << "HDIV_TET_" << dd;
908 CHKERR postProcMesh.tag_get_handle(
909 ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
910 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
911 }
912 break;
913 default:
914 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
915 "data inconsistency");
916 }
917
918 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
919 for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
920 CHKERR postProcMesh.tag_set_data(th[dd], &mapGaussPts[gg], 1,
921 &data.getVectorN<3>(gg)(dd, 0));
922 }
923 }
924
926 }
927 };
928
929private:
930 std::vector<MatrixDouble> levelShapeFunctionsTets;
931 std::vector<MatrixDouble> levelGaussPtsOnRefMeshTets;
932 std::vector<ublas::matrix<int>> levelRefTets;
933 std::vector<MatrixDouble> levelShapeFunctionsHexes;
934 std::vector<MatrixDouble> levelGaussPtsOnRefMeshHexes;
935 std::vector<ublas::matrix<int>> levelRefHexes;
936
938 std::vector<double *> verticesOnTetArrays;
943
945 std::vector<double *> verticesOnHexArrays;
950};
951
952/** \brief Post processing
953 * \ingroup mofem_fs_post_proc
954 */
957 MoFEM::VolumeElementForcesAndSourcesCore> {
958
961 PostProcTemplateVolumeOnRefinedMesh;
962};
963
964// /** \deprecated Use PostPocOnRefinedMesh instead
965// */
966// DEPRECATED typedef PostProcVolumeOnRefinedMesh PostPocOnRefinedMesh;
967
968/**
969 * \brief Postprocess on prism
970 *
971 * \ingroup mofem_fs_post_proc
972 */
975 MoFEM::FatPrismElementForcesAndSourcesCore> {
976
978
980 bool ten_nodes_post_proc_tets = true)
983 tenNodesPostProcTets(ten_nodes_post_proc_tets) {}
984
985 int getRuleTrianglesOnly(int order) { return -1; };
986 int getRuleThroughThickness(int order) { return -1; };
987
988 struct PointsMap3D {
989 const int kSi;
990 const int eTa;
991 const int zEta;
992 int nN;
993 PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
994 : kSi(ksi), eTa(eta), zEta(zeta), nN(nn) {}
995 };
996
997 typedef multi_index_container<
998 PointsMap3D,
999 indexed_by<ordered_unique<composite_key<
1000 PointsMap3D, member<PointsMap3D, const int, &PointsMap3D::kSi>,
1001 member<PointsMap3D, const int, &PointsMap3D::eTa>,
1002 member<PointsMap3D, const int, &PointsMap3D::zEta>>>>>
1004
1005 // PointsMap3D_multiIndex pointsMap;
1006
1007 MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only);
1008 MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness);
1009 MoFEMErrorCode generateReferenceElementMesh();
1010
1011 // std::map<EntityHandle, EntityHandle> elementsMap;
1012 std::map<EntityHandle, std::vector<EntityHandle>> elementsMap;
1013 std::map<EntityHandle, std::vector<PointsMap3D_multiIndex>>
1015
1016 MoFEMErrorCode preProcess();
1017 MoFEMErrorCode postProcess();
1018
1021
1023 return commonData;
1024 }
1025};
1026
1027/**
1028 * \brief Postprocess on face
1029 *
1030 * \ingroup mofem_fs_post_proc
1031 */
1033 MoFEM::FaceElementForcesAndSourcesCore> {
1034
1036
1038 bool six_node_post_proc_tris = true)
1040 m_field),
1041 sixNodePostProcTris(six_node_post_proc_tris), counterTris(0),
1042 counterQuads(0) {}
1043
1044 // Gauss pts set on refined mesh
1045 int getRule(int order) { return -1; };
1046
1047 MoFEMErrorCode generateReferenceElementMesh();
1048 MoFEMErrorCode setGaussPts(int order);
1049
1050 std::map<EntityHandle, EntityHandle> elementsMap;
1051
1052 MoFEMErrorCode preProcess();
1053 MoFEMErrorCode postProcess();
1054
1057
1059 return commonData;
1060 }
1061
1062 template <int RANK>
1064 : public FaceElementForcesAndSourcesCore::UserDataOperator {
1065
1066 moab::Interface &postProcMesh;
1067 std::vector<EntityHandle> &mapGaussPts;
1068 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideOpFe;
1069 const std::string feVolName;
1070 boost::shared_ptr<MatrixDouble> matPtr;
1071 const std::string tagName;
1072 const std::string fieldName;
1073 const bool saveOnTag;
1074
1076 moab::Interface &post_proc_mesh,
1077 std::vector<EntityHandle> &map_gauss_pts, const std::string field_name,
1078 const std::string tag_name,
1079 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> side_fe,
1080 const std::string vol_fe_name, boost::shared_ptr<MatrixDouble> mat_ptr,
1081 bool save_on_tag)
1084 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
1085 sideOpFe(side_fe), feVolName(vol_fe_name), matPtr(mat_ptr),
1086 tagName(tag_name), fieldName(field_name), saveOnTag(save_on_tag) {}
1087
1088 MoFEMErrorCode doWork(int side, EntityType type,
1089 EntitiesFieldData::EntData &data);
1090 };
1091
1094
1096 const std::string field_name, const std::string vol_fe_name = "dFE",
1097 boost::shared_ptr<MatrixDouble> grad_mat_ptr = nullptr,
1098 bool save_on_tag = true);
1099
1100 MoFEMErrorCode addFieldValuesPostProcOnSkin(
1101 const std::string field_name, const std::string vol_fe_name = "dFE",
1102 boost::shared_ptr<MatrixDouble> mat_ptr = nullptr,
1103 bool save_on_tag = true);
1104
1105private:
1106 MatrixDouble gaussPtsTri; ///< Gauss points coordinates on reference triangle
1107 MatrixDouble gaussPtsQuad; ///< Gauss points coordinates on reference quad
1108
1109 EntityHandle *triConn; ///< Connectivity for created tri elements
1110 EntityHandle *quadConn; ///< Connectivity for created quad elements
1112 startingVertTriHandle; ///< Starting handle for vertices on triangles
1114 startingVertQuadHandle; ///< Starting handle for vertices on quads
1115 std::vector<double *>
1116 verticesOnTriArrays; /// pointers to memory allocated by MoAB for
1117 /// storing X, Y, and Z coordinates
1118 std::vector<double *>
1119 verticesOnQuadArrays; /// pointers to memory allocated by MoAB for
1120 /// storing X, Y, and Z coordinates
1121
1123 startingEleTriHandle; ///< Starting handle for triangles post proc
1124 EntityHandle startingEleQuadHandle; ///< Starting handle for quads post proc
1125
1126 int numberOfTriangles; ///< Number of triangles to create
1127 int numberOfQuads; ///< NUmber of quads to create
1130};
1131
1132/**
1133 * @brief Postprocess 2d face elements
1134 *
1135 */
1137
1139};
1140
1142
1143/**
1144 * \brief Postprocess on edge
1145 *
1146 * \ingroup mofem_fs_post_proc
1147 */
1149 : public PostProcTemplateOnRefineMesh<EdgeEleBasePostProc> {
1150
1152
1154 bool six_node_post_proc_tris = true)
1156 sixNodePostProcTris(six_node_post_proc_tris) {}
1157
1158 // Gauss pts set on refined mesh
1159 int getRule(int order) { return -1; };
1160
1161 MoFEMErrorCode generateReferenceElementMesh();
1162 MoFEMErrorCode setGaussPts(int order);
1163
1164 std::map<EntityHandle, EntityHandle> elementsMap;
1165
1166 MoFEMErrorCode preProcess();
1167 MoFEMErrorCode postProcess();
1168
1171
1173 return commonData;
1174 }
1175};
1176
1177/**
1178 * @brief Post post-proc data at points from hash maps
1179 *
1180 * @tparam DIM1 dimension of vector in data_map_vec and first column of
1181 * data_map_may
1182 * @tparam DIM2 dimension of second column in data_map_mat
1183 */
1184template <int DIM1, int DIM2>
1185struct OpPostProcMap : public ForcesAndSourcesCore::UserDataOperator {
1186
1187 using DataMapVec = std::map<std::string, boost::shared_ptr<VectorDouble>>;
1188 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
1189
1190 /**
1191 * @brief Construct a new OpPostProcMap object
1192 *
1193 * @param post_proc_mesh postprocessing mesh
1194 * @param map_gauss_pts map of gauss points to nodes of postprocessing mesh
1195 * @param data_map_scalar hash map of scalar values (string is name of the
1196 * tag)
1197 * @param data_map_vec hash map of vector values
1198 * @param data_map_mat hash map of second order tensor values
1199 * @param data_symm_map_mat hash map of symmetric second order tensor values
1200 */
1201 OpPostProcMap(moab::Interface &post_proc_mesh,
1202 std::vector<EntityHandle> &map_gauss_pts,
1203 DataMapVec data_map_scalar, DataMapMat data_map_vec,
1204 DataMapMat data_map_mat, DataMapMat data_symm_map_mat)
1205 : ForcesAndSourcesCore::UserDataOperator(
1206 NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE),
1207 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
1208 dataMapScalar(data_map_scalar), dataMapVec(data_map_vec),
1209 dataMapMat(data_map_mat), dataMapSymmMat(data_symm_map_mat) {
1210 // Operator is only executed for vertices
1211 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1212 }
1213 MoFEMErrorCode doWork(int side, EntityType type,
1214 EntitiesFieldData::EntData &data);
1215
1216private:
1217 moab::Interface &postProcMesh;
1218 std::vector<EntityHandle> &mapGaussPts;
1223};
1224
1225template <int DIM1, int DIM2>
1226MoFEMErrorCode
1228 EntitiesFieldData::EntData &data) {
1230
1231 std::array<double, 9> def;
1232 std::fill(def.begin(), def.end(), 0);
1233
1234 auto get_tag = [&](const std::string name, size_t size) {
1235 Tag th;
1236 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
1237 MB_TAG_CREAT | MB_TAG_SPARSE,
1238 def.data());
1239 return th;
1240 };
1241
1242 MatrixDouble3by3 mat(3, 3);
1243
1244 auto set_vector_3d = [&](auto &t) -> MatrixDouble3by3 & {
1245 mat.clear();
1246 for (size_t r = 0; r != DIM1; ++r)
1247 mat(0, r) = t(r);
1248 return mat;
1249 };
1250
1251 auto set_matrix_3d = [&](auto &t) -> MatrixDouble3by3 & {
1252 mat.clear();
1253 for (size_t r = 0; r != DIM1; ++r)
1254 for (size_t c = 0; c != DIM2; ++c)
1255 mat(r, c) = t(r, c);
1256 return mat;
1257 };
1258
1259 auto set_matrix_symm_3d = [&](auto &t) -> MatrixDouble3by3 & {
1260 mat.clear();
1261 for (size_t r = 0; r != DIM1; ++r)
1262 for (size_t c = 0; c != DIM1; ++c)
1263 mat(r, c) = t(r, c);
1264 return mat;
1265 };
1266
1267 auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
1268 mat.clear();
1269 mat(0, 0) = t;
1270 return mat;
1271 };
1272
1273 auto set_float_precision = [](const double x) {
1274 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1275 return 0.;
1276 else
1277 return x;
1278 };
1279
1280 auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
1281 for (auto &v : mat.data())
1282 v = set_float_precision(v);
1283 return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1284 &*mat.data().begin());
1285 };
1286
1287 for (auto &m : dataMapScalar) {
1288 auto th = get_tag(m.first, 1);
1289 auto t_scl = getFTensor0FromVec(*m.second);
1290 auto nb_integration_pts = getGaussPts().size2();
1291 size_t gg = 0;
1292 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1293 CHKERR set_tag(th, gg, set_scalar(t_scl));
1294 ++t_scl;
1295 }
1296 }
1297
1298 for (auto &m : dataMapVec) {
1299 auto th = get_tag(m.first, 3);
1300 auto t_vec = getFTensor1FromMat<DIM1>(*m.second);
1301 auto nb_integration_pts = getGaussPts().size2();
1302 size_t gg = 0;
1303 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1304 CHKERR set_tag(th, gg, set_vector_3d(t_vec));
1305 ++t_vec;
1306 }
1307 }
1308
1309 for (auto &m : dataMapMat) {
1310 auto th = get_tag(m.first, 9);
1311 auto t_mat = getFTensor2FromMat<DIM1, DIM2>(*m.second);
1312 auto nb_integration_pts = getGaussPts().size2();
1313 size_t gg = 0;
1314 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1315 CHKERR set_tag(th, gg, set_matrix_3d(t_mat));
1316 ++t_mat;
1317 }
1318 }
1319
1320 for (auto &m : dataMapSymmMat) {
1321 auto th = get_tag(m.first, 9);
1322 auto t_mat = getFTensor2SymmetricFromMat<DIM1>(*m.second);
1323 auto nb_integration_pts = getGaussPts().size2();
1324 size_t gg = 0;
1325 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1326 CHKERR set_tag(th, gg, set_matrix_symm_3d(t_mat));
1327 ++t_mat;
1328 }
1329 }
1330
1332}
1333
1334#endif //__POSTPROC_ON_REF_MESH_HPP
1335
1336/**
1337 * \defgroup mofem_fs_post_proc Post Process
1338 * \ingroup user_modules
1339 **/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition: LogManager.hpp:370
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ NOSPACE
Definition: definitions.h:83
#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
constexpr int order
#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 N_MBHEX4(x, y, z)
Definition: fem_tools.h:75
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
#define N_MBHEX0(x, y, z)
Definition: fem_tools.h:71
#define N_MBHEX6(x, y, z)
Definition: fem_tools.h:77
#define N_MBHEX2(x, y, z)
Definition: fem_tools.h:73
#define N_MBHEX1(x, y, z)
Definition: fem_tools.h:72
FTensor::Index< 'm', SPACE_DIM > m
double eta
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, const std::string tag_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field value.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, int space_dim, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, const std::string tag_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
constexpr int nb_ref_levels
Three levels of refinement.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
constexpr int DIM2
Definition: level_set.cpp:22
constexpr int DIM1
Definition: level_set.cpp:21
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
double zeta
Viscous hardening.
Definition: plastic.cpp:177
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
Deprecated interface functions.
@ OPCOL
operator doWork function is executed on FE columns
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Post post-proc data at points from hash maps.
DataMapVec dataMapScalar
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
moab::Interface & postProcMesh
OpPostProcMap(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, DataMapVec data_map_scalar, DataMapMat data_map_vec, DataMapMat data_map_mat, DataMapMat data_symm_map_mat)
Construct a new OpPostProcMap object.
std::vector< EntityHandle > & mapGaussPts
DataMapMat dataMapVec
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
DataMapMat dataMapMat
DataMapMat dataMapSymmMat
Range tEts
std::map< std::string, std::vector< VectorDouble > > fieldMap
std::map< std::string, std::vector< MatrixDouble > > gradMap
operator to post-process (save gradients on refined post-processing mesh) field gradient
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::vector< EntityHandle > & mapGaussPts
moab::Interface & postProcMesh
const std::string tagName
VectorDouble vAlues
Vec V
CommonData & commonData
VectorDouble * vAluesPtr
int spaceDim
OpGetFieldGradientValues(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name, const std::string tag_name, CommonData &common_data, Vec v=PETSC_NULL, int space_dim=3)
operator to post-process (save gradients on refined post-processing mesh) field gradient
const std::string tagName
VectorDouble vAlues
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpGetFieldValues(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name, const std::string tag_name, CommonData &common_data, Vec v=PETSC_NULL)
std::vector< EntityHandle > & mapGaussPts
VectorDouble * vAluesPtr
moab::Interface & postProcMesh
Vec V
CommonData & commonData
Set of operators and data structures used for post-processing.
Postprocess on edge.
MoFEMErrorCode postProcess()
function is run at the end of loop
bool sixNodePostProcTris
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode setGaussPts(int order)
int getRule(int order)
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
PostProcEdgeOnRefinedMesh(MoFEM::Interface &m_field, bool six_node_post_proc_tris=true)
CommonData commonData
std::map< EntityHandle, EntityHandle > elementsMap
MoFEMErrorCode generateReferenceElementMesh()
const bool saveOnTag
const std::string feVolName
boost::shared_ptr< MatrixDouble > matPtr
moab::Interface & postProcMesh
OpGetFieldValuesOnSkinImpl(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name, const std::string tag_name, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > side_fe, const std::string vol_fe_name, boost::shared_ptr< MatrixDouble > mat_ptr, bool save_on_tag)
const std::string tagName
const std::string fieldName
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideOpFe
Postprocess 2d face elements.
Postprocess on face.
EntityHandle * quadConn
Connectivity for created quad elements.
int numberOfTriangles
Number of triangles to create.
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()
std::vector< double * > verticesOnTriArrays
MoFEMErrorCode preProcess()
function is run at the beginning of loop
int numberOfQuads
NUmber of quads to create.
std::map< EntityHandle, EntityHandle > elementsMap
std::vector< double * > verticesOnQuadArrays
EntityHandle startingVertQuadHandle
Starting handle for vertices on quads.
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
CommonData commonData
PostProcFaceOnRefinedMesh(MoFEM::Interface &m_field, bool six_node_post_proc_tris=true)
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)
MatrixDouble gaussPtsQuad
Gauss points coordinates on reference quad.
EntityHandle startingVertTriHandle
Starting handle for vertices on triangles.
bool sixNodePostProcTris
EntityHandle startingEleQuadHandle
Starting handle for quads post proc.
int counterTris
MoFEMErrorCode postProcess()
function is run at the end of loop
int getRule(int order)
int counterQuads
MatrixDouble gaussPtsTri
Gauss points coordinates on reference triangle.
EntityHandle startingEleTriHandle
Starting handle for triangles post proc.
EntityHandle * triConn
Connectivity for created tri elements.
PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
const int zEta
const int eTa
const int kSi
int nN
Postprocess on prism.
int getRuleThroughThickness(int order)
CommonData commonData
std::map< EntityHandle, std::vector< EntityHandle > > elementsMap
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
PostProcFatPrismOnRefinedMesh(MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets=true)
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
int getRuleTrianglesOnly(int order)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
std::map< EntityHandle, std::vector< PointsMap3D_multiIndex > > pointsMapVectorMap
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode postProcess()
function is run at the end of loop
bool tenNodesPostProcTets
Generic post-processing class.
std::vector< EntityHandle > mapGaussPts
boost::shared_ptr< WrapMPIComm > wrapRefMeshComm
virtual ~PostProcTemplateOnRefineMesh()
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
moab::Interface & postProcMesh
moab::Core coreMesh
PostProcTemplateOnRefineMesh(MoFEM::Interface &m_field)
moab::Interface & postProcMesh
OpHdivFunctions(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::vector< EntityHandle > & mapGaussPts
int countVertHex
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
PostProcCommonOnRefMesh::CommonDataForVolume CommonData
std::vector< MatrixDouble > levelGaussPtsOnRefMeshTets
std::vector< ublas::matrix< int > > levelRefTets
size_t getMaxLevel() const
EntityHandle startingEleHexHandle
EntityHandle startingVertHexHandle
std::vector< MatrixDouble > levelShapeFunctionsHexes
int countVertTet
PostProcTemplateVolumeOnRefinedMesh(MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets=true, int nb_ref_levels=-1)
std::vector< double * > verticesOnHexArrays
MoFEMErrorCode setGaussPts(int order)
Set integration points.
CommonData commonData
int nbOfRefLevels
std::vector< double * > verticesOnTetArrays
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
bool tenNodesPostProcTets
std::vector< MatrixDouble > levelGaussPtsOnRefMeshHexes
MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name)
Add operator to post-process Hdiv field.
EntityHandle startingVertTetHandle
PostProcTemplateOnRefineMesh< VOLUME_ELEMENT > T
MoFEMErrorCode preProcess()
std::vector< MatrixDouble > levelShapeFunctionsTets
int countHex
EntityHandle * hexConn
std::vector< ublas::matrix< int > > levelRefHexes
MoFEMErrorCode postProcess()
EntityHandle * tetConn
int countTet
MoFEMErrorCode clearOperators()
Clear operators list.
EntityHandle startingEleTetHandle
int getRule(int order)
Post processing.