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 */
959 MoFEM::VolumeElementForcesAndSourcesCore> {
960
963 PostProcTemplateVolumeOnRefinedMesh;
964};
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
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
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()
Post-processing function executed at loop completion.
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
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()
Pre-processing function executed at loop initialization.
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()
Post-processing function executed at loop completion.
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)
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)
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
std::map< EntityHandle, std::vector< PointsMap3D_multiIndex > > pointsMapVectorMap
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
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.
PostProcCommonOnRefMesh::CommonDataForVolume CommonData
std::vector< MatrixDouble > levelGaussPtsOnRefMeshTets
std::vector< ublas::matrix< int > > levelRefTets
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.
PostProcTemplateOnRefineMesh< VOLUME_ELEMENT > T
std::vector< MatrixDouble > levelShapeFunctionsTets
std::vector< ublas::matrix< int > > levelRefHexes
MoFEMErrorCode clearOperators()
Clear operators list.