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