v0.15.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
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_NULLPTR)
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_NULLPTR,
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
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_NULLPTR) {
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_NULLPTR) {
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_NULLPTR) {
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_NULLPTR) {
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_NULLPTR) {
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
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_NULLPTR, PETSC_NULLPTR,
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
602 FTensor::Index<'i', 3> i;
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:
662
665
666 );
667 break;
668 case MBHEX:
671
674
675 );
676 break;
677 default:
678 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
679 "Element type not implemented");
680 }
681
683 }
684
685 /** \brief Clear operators list
686
687 Clear operators list, user can use the same mesh instance to post-process
688 different problem or the same problem with different set of post-processed
689 fields.
690
691 */
692 MoFEMErrorCode clearOperators() {
694 T::getOpPtrVector().clear();
696 }
697
698 MoFEMErrorCode preProcess() {
700 moab::Interface &moab = T::coreMesh;
701 ParallelComm *pcomm_post_proc_mesh =
702 ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
703 if (pcomm_post_proc_mesh != NULL)
704 delete pcomm_post_proc_mesh;
705
706 CHKERR T::postProcMesh.delete_mesh();
707
708 auto alloc_vertices_and_elements_on_post_proc_mesh = [&]() {
710
711 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
712
713 auto miit =
714 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
715 boost::make_tuple(this->getFEName(), this->getLoFERank()));
716 auto hi_miit =
717 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
718 boost::make_tuple(this->getFEName(), this->getHiFERank()));
719
720 const int number_of_ents_in_the_loop = this->getLoopSize();
721 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
722 SETERRQ(this->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
723 "Wrong size of indicices. Inconsistent size number of iterated "
724 "elements iterated by problem and from range.");
725 }
726
727 int nb_tet_vertices = 0;
728 int nb_tets = 0;
729 int nb_hex_vertices = 0;
730 int nb_hexes = 0;
731
732 for (; miit != hi_miit; ++miit) {
733 auto type = (*miit)->getEntType();
734
735 // Set pointer to element. So that getDataVectorDofsPtr in getMaxLevel
736 // can work
737 this->numeredEntFiniteElementPtr = *miit;
738
739 bool add = true;
740 if (this->exeTestHook) {
741 add = this->exeTestHook(this);
742 }
743
744 if (add) {
745
746 auto level = getMaxLevel();
747
748 switch (type) {
749 case MBTET:
750 level = std::min(level, levelGaussPtsOnRefMeshTets.size() - 1);
751 nb_tet_vertices += levelGaussPtsOnRefMeshTets[level].size2();
752 nb_tets += levelRefTets[level].size1();
753 break;
754 case MBHEX:
755 level = std::min(level, levelGaussPtsOnRefMeshHexes.size() - 1);
756 nb_hex_vertices += levelGaussPtsOnRefMeshHexes[level].size2();
757 nb_hexes += levelRefHexes[level].size1();
758 break;
759 default:
760 SETERRQ(this->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
761 "Element type not implemented");
762 break;
763 }
764 }
765 }
766
767 ReadUtilIface *iface;
768 CHKERR this->postProcMesh.query_interface(iface);
769
770 if (nb_tets) {
771 CHKERR iface->get_node_coords(
772 3, nb_tet_vertices, 0, startingVertTetHandle, verticesOnTetArrays);
773 CHKERR iface->get_element_connect(nb_tets, levelRefTets[0].size2(),
774 MBTET, 0, startingEleTetHandle,
775 tetConn);
776 }
777
778 if (nb_hexes) {
779 CHKERR iface->get_node_coords(
780 3, nb_hex_vertices, 0, startingVertHexHandle, verticesOnHexArrays);
781 CHKERR iface->get_element_connect(nb_hexes, levelRefHexes[0].size2(),
782 MBHEX, 0, startingEleHexHandle,
783 hexConn);
784 }
785
786 countTet = 0;
787 countVertTet = 0;
788 countHex = 0;
789 countVertHex = 0;
790
792 };
793
794 CHKERR alloc_vertices_and_elements_on_post_proc_mesh();
795
797 }
798
799 MoFEMErrorCode postProcess() {
801
802 auto update_elements = [&]() {
804 ReadUtilIface *iface;
805 CHKERR this->postProcMesh.query_interface(iface);
806
807 if (countTet) {
808 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
809 << "Update tets " << countTet;
810
811 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
812 << "Nb nodes on tets " << levelRefTets[0].size2();
813
814 CHKERR iface->update_adjacencies(startingEleTetHandle, countTet,
815 levelRefTets[0].size2(), tetConn);
816 }
817 if (countHex) {
818 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
819 << "Update hexes " << countHex;
820 CHKERR iface->update_adjacencies(startingEleHexHandle, countHex,
821 levelRefHexes[0].size2(), hexConn);
822 }
824 };
825
826 auto resolve_shared_ents = [&]() {
828 ParallelComm *pcomm_post_proc_mesh =
829 ParallelComm::get_pcomm(&(T::postProcMesh), MYPCOMM_INDEX);
830 if (pcomm_post_proc_mesh == NULL) {
831 // T::wrapRefMeshComm =
832 // boost::make_shared<WrapMPIComm>(T::mField.get_comm(), false);
833 pcomm_post_proc_mesh = new ParallelComm(
835 PETSC_COMM_WORLD /*(T::wrapRefMeshComm)->get_comm()*/);
836 }
837
838 Range edges;
839 CHKERR T::postProcMesh.get_entities_by_type(0, MBEDGE, edges, false);
840 CHKERR T::postProcMesh.delete_entities(edges);
841 Range faces;
842 CHKERR T::postProcMesh.get_entities_by_dimension(0, 2, faces, false);
843 CHKERR T::postProcMesh.delete_entities(faces);
844
845 Range ents;
846 CHKERR T::postProcMesh.get_entities_by_dimension(0, 3, ents, false);
847
848 int rank = T::mField.get_comm_rank();
849 CHKERR T::postProcMesh.tag_clear_data(pcomm_post_proc_mesh->part_tag(),
850 ents, &rank);
851
852 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
853
855 };
856
857 CHKERR resolve_shared_ents();
858 CHKERR update_elements();
859
861 }
862
863 /** \brief Add operator to post-process Hdiv field
864 */
865 MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name) {
867 T::getOpPtrVector().push_back(
870 }
871
872 struct OpHdivFunctions : public VOLUME_ELEMENT::UserDataOperator {
873
874 moab::Interface &postProcMesh;
875 std::vector<EntityHandle> &mapGaussPts;
876
877 OpHdivFunctions(moab::Interface &post_proc_mesh,
878 std::vector<EntityHandle> &map_gauss_pts,
879 const std::string field_name)
880 : VOLUME_ELEMENT::UserDataOperator(field_name,
881 T::UserDataOperator::OPCOL),
882 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
883
884 MoFEMErrorCode doWork(int side, EntityType type,
885 EntitiesFieldData::EntData &data) {
887
888 if (data.getIndices().size() == 0)
890
891 std::vector<Tag> th;
892 th.resize(data.getFieldData().size());
893
894 double def_VAL[9] = {0, 0, 0};
895
896 switch (type) {
897 case MBTRI:
898 for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
899 std::ostringstream ss;
900 ss << "HDIV_FACE_" << side << "_" << dd;
901 CHKERR postProcMesh.tag_get_handle(
902 ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
903 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
904 }
905 break;
906 case MBTET:
907 for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
908 std::ostringstream ss;
909 ss << "HDIV_TET_" << dd;
910 CHKERR postProcMesh.tag_get_handle(
911 ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
912 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
913 }
914 break;
915 default:
916 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
917 "data inconsistency");
918 }
919
920 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
921 for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
922 CHKERR postProcMesh.tag_set_data(th[dd], &mapGaussPts[gg], 1,
923 &data.getVectorN<3>(gg)(dd, 0));
924 }
925 }
926
928 }
929 };
930
931private:
932 std::vector<MatrixDouble> levelShapeFunctionsTets;
933 std::vector<MatrixDouble> levelGaussPtsOnRefMeshTets;
934 std::vector<ublas::matrix<int>> levelRefTets;
935 std::vector<MatrixDouble> levelShapeFunctionsHexes;
936 std::vector<MatrixDouble> levelGaussPtsOnRefMeshHexes;
937 std::vector<ublas::matrix<int>> levelRefHexes;
938
940 std::vector<double *> verticesOnTetArrays;
945
947 std::vector<double *> verticesOnHexArrays;
952};
953
954/** \brief Post processing
955 * \ingroup mofem_fs_post_proc
956 */
965
966// /** \deprecated Use PostPocOnRefinedMesh instead
967// */
968// DEPRECATED typedef PostProcVolumeOnRefinedMesh PostPocOnRefinedMesh;
969
970/**
971 * \brief Postprocess on prism
972 *
973 * \ingroup mofem_fs_post_proc
974 */
977 MoFEM::FatPrismElementForcesAndSourcesCore> {
978
980
982 bool ten_nodes_post_proc_tets = true)
985 tenNodesPostProcTets(ten_nodes_post_proc_tets) {}
986
987 int getRuleTrianglesOnly(int order) { return -1; };
988 int getRuleThroughThickness(int order) { return -1; };
989
990 struct PointsMap3D {
991 const int kSi;
992 const int eTa;
993 const int zEta;
994 int nN;
995 PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
996 : kSi(ksi), eTa(eta), zEta(zeta), nN(nn) {}
997 };
998
999 typedef multi_index_container<
1000 PointsMap3D,
1001 indexed_by<ordered_unique<composite_key<
1002 PointsMap3D, member<PointsMap3D, const int, &PointsMap3D::kSi>,
1003 member<PointsMap3D, const int, &PointsMap3D::eTa>,
1004 member<PointsMap3D, const int, &PointsMap3D::zEta>>>>>
1006
1007 // PointsMap3D_multiIndex pointsMap;
1008
1009 MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only);
1010 MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness);
1011 MoFEMErrorCode generateReferenceElementMesh();
1012
1013 // std::map<EntityHandle, EntityHandle> elementsMap;
1014 std::map<EntityHandle, std::vector<EntityHandle>> elementsMap;
1015 std::map<EntityHandle, std::vector<PointsMap3D_multiIndex>>
1017
1018 MoFEMErrorCode preProcess();
1019 MoFEMErrorCode postProcess();
1020
1023
1027};
1028
1029/**
1030 * \brief Postprocess on face
1031 *
1032 * \ingroup mofem_fs_post_proc
1033 */
1035 MoFEM::FaceElementForcesAndSourcesCore> {
1036
1038
1040 bool six_node_post_proc_tris = true)
1042 m_field),
1043 sixNodePostProcTris(six_node_post_proc_tris), counterTris(0),
1044 counterQuads(0) {}
1045
1046 // Gauss pts set on refined mesh
1047 int getRule(int order) { return -1; };
1048
1049 MoFEMErrorCode generateReferenceElementMesh();
1050 MoFEMErrorCode setGaussPts(int order);
1051
1052 std::map<EntityHandle, EntityHandle> elementsMap;
1053
1054 MoFEMErrorCode preProcess();
1055 MoFEMErrorCode postProcess();
1056
1059
1063
1064 template <int RANK>
1066 : public FaceElementForcesAndSourcesCore::UserDataOperator {
1067
1068 moab::Interface &postProcMesh;
1069 std::vector<EntityHandle> &mapGaussPts;
1070 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideOpFe;
1071 const std::string feVolName;
1072 boost::shared_ptr<MatrixDouble> matPtr;
1073 const std::string tagName;
1074 const std::string fieldName;
1075 const bool saveOnTag;
1076
1078 moab::Interface &post_proc_mesh,
1079 std::vector<EntityHandle> &map_gauss_pts, const std::string field_name,
1080 const std::string tag_name,
1081 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> side_fe,
1082 const std::string vol_fe_name, boost::shared_ptr<MatrixDouble> mat_ptr,
1083 bool save_on_tag)
1086 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
1087 sideOpFe(side_fe), feVolName(vol_fe_name), matPtr(mat_ptr),
1088 tagName(tag_name), fieldName(field_name), saveOnTag(save_on_tag) {}
1089
1090 MoFEMErrorCode doWork(int side, EntityType type,
1091 EntitiesFieldData::EntData &data);
1092 };
1093
1096
1098 const std::string field_name, const std::string vol_fe_name = "dFE",
1099 boost::shared_ptr<MatrixDouble> grad_mat_ptr = nullptr,
1100 bool save_on_tag = true);
1101
1102 MoFEMErrorCode addFieldValuesPostProcOnSkin(
1103 const std::string field_name, const std::string vol_fe_name = "dFE",
1104 boost::shared_ptr<MatrixDouble> mat_ptr = nullptr,
1105 bool save_on_tag = true);
1106
1107private:
1108 MatrixDouble gaussPtsTri; ///< Gauss points coordinates on reference triangle
1109 MatrixDouble gaussPtsQuad; ///< Gauss points coordinates on reference quad
1110
1111 EntityHandle *triConn; ///< Connectivity for created tri elements
1112 EntityHandle *quadConn; ///< Connectivity for created quad elements
1114 startingVertTriHandle; ///< Starting handle for vertices on triangles
1116 startingVertQuadHandle; ///< Starting handle for vertices on quads
1117 std::vector<double *>
1118 verticesOnTriArrays; /// pointers to memory allocated by MoAB for
1119 /// storing X, Y, and Z coordinates
1120 std::vector<double *>
1121 verticesOnQuadArrays; /// pointers to memory allocated by MoAB for
1122 /// storing X, Y, and Z coordinates
1123
1125 startingEleTriHandle; ///< Starting handle for triangles post proc
1126 EntityHandle startingEleQuadHandle; ///< Starting handle for quads post proc
1127
1128 int numberOfTriangles; ///< Number of triangles to create
1129 int numberOfQuads; ///< NUmber of quads to create
1132};
1133
1134/**
1135 * @brief Postprocess 2d face elements
1136 *
1137 */
1142
1144
1145/**
1146 * \brief Postprocess on edge
1147 *
1148 * \ingroup mofem_fs_post_proc
1149 */
1151 : public PostProcTemplateOnRefineMesh<EdgeEleBasePostProc> {
1152
1154
1156 bool six_node_post_proc_tris = true)
1158 sixNodePostProcTris(six_node_post_proc_tris) {}
1159
1160 // Gauss pts set on refined mesh
1161 int getRule(int order) { return -1; };
1162
1163 MoFEMErrorCode generateReferenceElementMesh();
1164 MoFEMErrorCode setGaussPts(int order);
1165
1166 std::map<EntityHandle, EntityHandle> elementsMap;
1167
1168 MoFEMErrorCode preProcess();
1169 MoFEMErrorCode postProcess();
1170
1173
1177};
1178
1179/**
1180 * @brief Post post-proc data at points from hash maps
1181 *
1182 * @tparam DIM1 dimension of vector in data_map_vec and first column of
1183 * data_map_may
1184 * @tparam DIM2 dimension of second column in data_map_mat
1185 */
1186template <int DIM1, int DIM2>
1187struct OpPostProcMap : public ForcesAndSourcesCore::UserDataOperator {
1188
1189 using DataMapVec = std::map<std::string, boost::shared_ptr<VectorDouble>>;
1190 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
1191
1192 /**
1193 * @brief Construct a new OpPostProcMap object
1194 *
1195 * @param post_proc_mesh postprocessing mesh
1196 * @param map_gauss_pts map of gauss points to nodes of postprocessing mesh
1197 * @param data_map_scalar hash map of scalar values (string is name of the
1198 * tag)
1199 * @param data_map_vec hash map of vector values
1200 * @param data_map_mat hash map of second order tensor values
1201 * @param data_symm_map_mat hash map of symmetric second order tensor values
1202 */
1203 OpPostProcMap(moab::Interface &post_proc_mesh,
1204 std::vector<EntityHandle> &map_gauss_pts,
1205 DataMapVec data_map_scalar, DataMapMat data_map_vec,
1206 DataMapMat data_map_mat, DataMapMat data_symm_map_mat)
1207 : ForcesAndSourcesCore::UserDataOperator(
1208 NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE),
1209 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
1210 dataMapScalar(data_map_scalar), dataMapVec(data_map_vec),
1211 dataMapMat(data_map_mat), dataMapSymmMat(data_symm_map_mat) {
1212 // Operator is only executed for vertices
1213 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1214 }
1215 MoFEMErrorCode doWork(int side, EntityType type,
1216 EntitiesFieldData::EntData &data);
1217
1218private:
1219 moab::Interface &postProcMesh;
1220 std::vector<EntityHandle> &mapGaussPts;
1225};
1226
1227template <int DIM1, int DIM2>
1228MoFEMErrorCode
1229OpPostProcMap<DIM1, DIM2>::doWork(int side, EntityType type,
1230 EntitiesFieldData::EntData &data) {
1232
1233 std::array<double, 9> def;
1234 std::fill(def.begin(), def.end(), 0);
1235
1236 auto get_tag = [&](const std::string name, size_t size) {
1237 Tag th;
1238 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
1239 MB_TAG_CREAT | MB_TAG_SPARSE,
1240 def.data());
1241 return th;
1242 };
1243
1244 MatrixDouble3by3 mat(3, 3);
1245
1246 auto set_vector_3d = [&](auto &t) -> MatrixDouble3by3 & {
1247 mat.clear();
1248 for (size_t r = 0; r != DIM1; ++r)
1249 mat(0, r) = t(r);
1250 return mat;
1251 };
1252
1253 auto set_matrix_3d = [&](auto &t) -> MatrixDouble3by3 & {
1254 mat.clear();
1255 for (size_t r = 0; r != DIM1; ++r)
1256 for (size_t c = 0; c != DIM2; ++c)
1257 mat(r, c) = t(r, c);
1258 return mat;
1259 };
1260
1261 auto set_matrix_symm_3d = [&](auto &t) -> MatrixDouble3by3 & {
1262 mat.clear();
1263 for (size_t r = 0; r != DIM1; ++r)
1264 for (size_t c = 0; c != DIM1; ++c)
1265 mat(r, c) = t(r, c);
1266 return mat;
1267 };
1268
1269 auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
1270 mat.clear();
1271 mat(0, 0) = t;
1272 return mat;
1273 };
1274
1275 auto set_float_precision = [](const double x) {
1276 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1277 return 0.;
1278 else
1279 return x;
1280 };
1281
1282 auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
1283 for (auto &v : mat.data())
1284 v = set_float_precision(v);
1285 return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1286 &*mat.data().begin());
1287 };
1288
1289 for (auto &m : dataMapScalar) {
1290 auto th = get_tag(m.first, 1);
1291 auto t_scl = getFTensor0FromVec(*m.second);
1292 auto nb_integration_pts = getGaussPts().size2();
1293 size_t gg = 0;
1294 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1295 CHKERR set_tag(th, gg, set_scalar(t_scl));
1296 ++t_scl;
1297 }
1298 }
1299
1300 for (auto &m : dataMapVec) {
1301 auto th = get_tag(m.first, 3);
1302 auto t_vec = getFTensor1FromMat<DIM1>(*m.second);
1303 auto nb_integration_pts = getGaussPts().size2();
1304 size_t gg = 0;
1305 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1306 CHKERR set_tag(th, gg, set_vector_3d(t_vec));
1307 ++t_vec;
1308 }
1309 }
1310
1311 for (auto &m : dataMapMat) {
1312 auto th = get_tag(m.first, 9);
1313 auto t_mat = getFTensor2FromMat<DIM1, DIM2>(*m.second);
1314 auto nb_integration_pts = getGaussPts().size2();
1315 size_t gg = 0;
1316 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1317 CHKERR set_tag(th, gg, set_matrix_3d(t_mat));
1318 ++t_mat;
1319 }
1320 }
1321
1322 for (auto &m : dataMapSymmMat) {
1323 auto th = get_tag(m.first, 9);
1324 auto t_mat = getFTensor2SymmetricFromMat<DIM1>(*m.second);
1325 auto nb_integration_pts = getGaussPts().size2();
1326 size_t gg = 0;
1327 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1328 CHKERR set_tag(th, gg, set_matrix_symm_3d(t_mat));
1329 ++t_mat;
1330 }
1331 }
1332
1334}
1335
1336#endif //__POSTPROC_ON_REF_MESH_HPP
1337
1338/**
1339 * \defgroup mofem_fs_post_proc Post Process
1340 * \ingroup user_modules
1341 **/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOSPACE
Definition definitions.h:83
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ 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()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
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
double eta
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, const std::string tag_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2 or H1 field value.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, const std::string tag_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2 or H1 field gradient.
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_NULLPTR)
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:130
constexpr double t
plate stiffness
Definition plate.cpp:58
constexpr auto field_name
FTensor::Index< 'm', 3 > m
Deprecated interface functions.
@ OPCOL
operator doWork function is executed on FE columns
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Post post-proc data at points from hash maps.
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
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
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)
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_NULLPTR, int space_dim=3)
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
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_NULLPTR)
Set of operators and data structures used for post-processing.
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode setGaussPts(int order)
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
PostProcEdgeOnRefinedMesh(MoFEM::Interface &m_field, bool six_node_post_proc_tris=true)
std::map< EntityHandle, EntityHandle > elementsMap
MoFEMErrorCode generateReferenceElementMesh()
Postprocess 2d face elements.
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)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideOpFe
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()
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.
EntityHandle startingEleQuadHandle
Starting handle for quads post proc.
MoFEMErrorCode postProcess()
function is run at the end of loop
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)
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
std::map< EntityHandle, std::vector< EntityHandle > > elementsMap
PostProcFatPrismOnRefinedMesh(MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets=true)
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
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
Generic post-processing class.
std::vector< EntityHandle > mapGaussPts
boost::shared_ptr< WrapMPIComm > wrapRefMeshComm
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
PostProcTemplateOnRefineMesh(MoFEM::Interface &m_field)
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)
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
std::vector< MatrixDouble > levelGaussPtsOnRefMeshTets
std::vector< ublas::matrix< int > > levelRefTets
PostProcTemplateOnRefineMesh< VOLUME_ELEMENT > T
std::vector< MatrixDouble > levelShapeFunctionsHexes
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.
std::vector< double * > verticesOnTetArrays
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
std::vector< MatrixDouble > levelGaussPtsOnRefMeshHexes
MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name)
Add operator to post-process Hdiv field.
std::vector< MatrixDouble > levelShapeFunctionsTets
std::vector< ublas::matrix< int > > levelRefHexes
PostProcCommonOnRefMesh::CommonDataForVolume CommonData
MoFEMErrorCode clearOperators()
Clear operators list.