v0.14.0
Loading...
Searching...
No Matches
PostProcBrokenMeshInMoabBase.hpp
Go to the documentation of this file.
1/**
2 * @file PostProcBrokenMeshInMoabBase.hpp
3 * @author your name (you@domain.com)
4 * @brief
5 * @version 0.1
6 * @date 2022-09-14
7 *
8 * @copyright Copyright (c) 2022
9 *
10 */
11
12#ifndef __POSTPROCBROKENMESHINMOABBASE_HPP
13#define __POSTPROCBROKENMESHINMOABBASE_HPP
14
15namespace MoFEM {
16
17/**
18 * Each element is subdivided on smaller elements, i.e. a reference mesh on
19 * single element is created. Nodes of such reference mesh are used as
20 * integration points at which field values are calculated and to
21 * each node a "moab" tag is attached to store those values.
22 */
24
25 std::vector<MatrixDouble>
26 levelShapeFunctions; //< values of shape functions on refences element
27 // refinement levels
28 std::vector<MatrixDouble>
29 levelGaussPtsOnRefMesh; //< gauss points at refinement levels
30 std::vector<ublas::matrix<int>> levelRef; //< connectivity at refinement level
31
32 EntityHandle startingVertEleHandle; //< starting handle of first vertex
33 std::vector<double *> verticesOnEleArrays; //< array of vertices coordinate
34 EntityHandle startingEleHandle; //< starting handle of first element
35 EntityHandle *eleConn; //< array of elements connectivity
36
37 int countEle; //< count elements
38 int countVertEle; //< count vertices on the mesh
39
40 int nbVertices; //< numer of vertices on the mesh
41 int nbEles; //< number of elements on the mesh
42
44 virtual ~PostProcGenerateRefMeshBase() = default;
45
46 virtual MoFEMErrorCode getOptions(std::string prefix); //< get options from command
48
49 PetscBool hoNodes; //< if true mid nodes are added
50 int defMaxLevel; //< default max number of refinement levels
51};
52
54 boost::shared_ptr<PostProcGenerateRefMeshBase>;
55
56/**
57 * @brief Element for postprocessing. Uses MoAB to generate post-processing
58 * mesh.
59 *
60 * @tparam T Finite Element Implementation
61 */
62template <EntityType T> struct PostProcGenerateRefMesh;
63
64template <>
67 MoFEMErrorCode generateReferenceElementMesh();
68};
69
70template <>
73 MoFEMErrorCode generateReferenceElementMesh();
74};
75
76template <>
79 MoFEMErrorCode generateReferenceElementMesh();
80};
81
82template <>
85 MoFEMErrorCode generateReferenceElementMesh();
86};
87
88template <>
91 MoFEMErrorCode generateReferenceElementMesh();
92};
93
94template <typename E> struct PostProcBrokenMeshInMoabBase : public E {
95
97 std::string opts_prefix = "");
99 boost::shared_ptr<moab::Core> core_mesh_ptr,
100 std::string opts_prefix = "");
101
103
104 /**
105 * @brief Get vector of vectors associated to integration points
106 *
107 * @return std::vector<EntityHandle>&
108 */
109 inline auto &getMapGaussPts();
110
111 /**
112 * @brief Get postprocessing mesh
113 *
114 * @return moab::Interface&
115 */
116 inline auto &getPostProcMesh();
117
118 /**
119 * @brief Get postprocessing elements
120 *
121 * @return auto&
122 */
123 inline auto &getPostProcElements();
124
125 /**
126 * \brief wrote results in (MOAB) format, use "file_name.h5m"
127 * @param file_name file name (should always end with .h5m)
128 * @return error code
129 * \ingroup mofem_fs_post_proc
130 */
131 MoFEMErrorCode writeFile(const std::string file_name);
132
133 /**
134 * @brief Set tags to be transferred to post-processing mesh
135 *
136 * @param tags_to_transfer
137 * @return MoFEMErrorCode
138 */
139 MoFEMErrorCode setTagsToTransfer(std::vector<Tag> tags_to_transfer);
140
141protected:
142
144
145 /**
146 * @brief Generate vertices and elements
147 *
148 * @return MoFEMErrorCode
149 */
151
153
154 int getRule(int order);
155
156 /**
157 * @brief Determine refinement level based on fields approx ordre.
158 *
159 * level = (order - 1) / 2
160 *
161 * @return int
162 */
163 virtual int getMaxLevel() const;
164
165 boost::shared_ptr<moab::Core> coreMeshPtr = boost::make_shared<moab::Core>();
166
167 std::vector<EntityHandle> mapGaussPts;
169
170 std::map<EntityType, PostProcGenerateRefMeshPtr>
171 refElementsMap; ///< Storing data about element types, and data for
172 ///< ReadUtilIface to create element entities on
173 ///< post-process mesh.
174 std::vector<Tag> tagsToTransfer; ///< Set of tags on mesh to transfer to
175 ///< postprocessing mesh
176 std::string optionsPrefix = ""; ///< Prefix for options
177
180
181 virtual MoFEMErrorCode transferTags(); ///< transfer tags from mesh to
182 ///< post-process mesh
183};
184
185template <typename E>
187 std::vector<Tag> tags_to_transfer) {
189 tagsToTransfer.swap(tags_to_transfer);
191}
192
193template <typename E> int PostProcBrokenMeshInMoabBase<E>::getMaxLevel() const {
194 auto get_element_max_dofs_order = [&]() {
195 int max_order = 0;
196 auto dofs_vec = E::getDataVectorDofsPtr();
197 for (auto &dof : *dofs_vec) {
198 const int dof_order = dof->getDofOrder();
199 max_order = (max_order < dof_order) ? dof_order : max_order;
200 };
201 return max_order;
202 };
203 const auto dof_max_order = get_element_max_dofs_order();
204 return (dof_max_order > 0) ? (dof_max_order - 1) / 2 : 0;
205};
206
207template <typename E>
210
211 auto &calc_mesh = this->mField.get_moab();
212
213 auto name = [&](auto tag) {
214 std::string name;
215 CHK_MOAB_THROW(calc_mesh.tag_get_name(tag, name), "get name");
216 return name;
217 };
218
219 auto data_type = [&](auto tag) {
220 moab::DataType data_type;
221 CHK_MOAB_THROW(calc_mesh.tag_get_data_type(tag, data_type),
222 "get data type");
223 return data_type;
224 };
225
226 auto type = [&](auto tag) {
227 moab::TagType type;
228 CHK_MOAB_THROW(calc_mesh.tag_get_type(tag, type), "get tag type");
229 return type;
230 };
231
232 auto length = [&](auto tag) {
233 int length;
234 CHK_MOAB_THROW(calc_mesh.tag_get_length(tag, length), "get length ");
235 return length;
236 };
237
238 auto default_value = [&](auto tag) {
239 const void *def_val;
240 int size;
241 CHK_MOAB_THROW(calc_mesh.tag_get_default_value(tag, def_val, size),
242 "get default tag value");
243 return def_val;
244 };
245
246 auto data_ptr = [&](auto tag) {
247 const void *tag_data;
248 auto core_ent = this->getFEEntityHandle();
249 CHK_MOAB_THROW(calc_mesh.tag_get_by_ptr(tag, &core_ent, 1, &tag_data),
250 "get tag data");
251 return tag_data;
252 };
253
254 for (auto tag : tagsToTransfer) {
255 Tag tag_postproc;
256 CHKERR getPostProcMesh().tag_get_handle(
257 name(tag).c_str(), length(tag), data_type(tag), tag_postproc,
258 type(tag) | MB_TAG_CREAT, default_value(tag));
259 CHKERR getPostProcMesh().tag_clear_data(tag_postproc, postProcElements,
260 data_ptr(tag));
261 }
262
264};
265
266template <typename E>
268 MoFEM::Interface &m_field, std::string opts_prefix)
269 : E(m_field), optionsPrefix(opts_prefix) {}
270
271template <typename E>
273 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
274 std::string opts_prefix)
275 : PostProcBrokenMeshInMoabBase(m_field, opts_prefix) {
276 coreMeshPtr = core_mesh_ptr;
277}
278
279template <typename E>
281 ParallelComm *pcomm_post_proc_mesh =
282 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
283 if (pcomm_post_proc_mesh != NULL)
284 delete pcomm_post_proc_mesh;
285}
286
287template <typename E> int PostProcBrokenMeshInMoabBase<E>::getRule(int order) {
288 return -1;
289};
290
291template <typename E>
294
295 auto type = type_from_handle(this->getFEEntityHandle());
296
298
299 try {
300 ref_ele = refElementsMap.at(type);
301 } catch (const out_of_range &e) {
302 SETERRQ1(
303 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
304 "Generation of reference elements for type <%s> is not implemented",
305 moab::CN::EntityTypeName(type));
306 }
307
308 auto set_gauss_pts = [&](auto &level_gauss_pts_on_ref_mesh, auto &level_ref,
309 auto &level_shape_functions,
310
311 auto start_vert_handle, auto start_ele_handle,
312 auto &verts_array, auto &conn, auto &ver_count,
313 auto &ele_count
314
315 ) {
317
318 size_t level = getMaxLevel();
319 level = std::min(level, level_gauss_pts_on_ref_mesh.size() - 1);
320
321 auto &level_ref_gauss_pts = level_gauss_pts_on_ref_mesh[level];
322 auto &level_ref_ele = level_ref[level];
323 auto &shape_functions = level_shape_functions[level];
324 E::gaussPts.resize(level_ref_gauss_pts.size1(), level_ref_gauss_pts.size2(),
325 false);
326 noalias(E::gaussPts) = level_ref_gauss_pts;
327
328 const auto fe_ent = E::numeredEntFiniteElementPtr->getEnt();
329 auto get_fe_coords = [&]() {
330 const EntityHandle *conn;
331 int num_nodes;
333 E::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes, true),
334 "error get connectivity");
335 VectorDouble coords(num_nodes * 3);
337 E::mField.get_moab().get_coords(conn, num_nodes, &*coords.begin()),
338 "error get coordinates");
339 return coords;
340 };
341
342 auto coords = get_fe_coords();
343
344 const int num_nodes = level_ref_gauss_pts.size2();
345 mapGaussPts.resize(level_ref_gauss_pts.size2());
346
349 &*shape_functions.data().begin());
351 &verts_array[0][ver_count], &verts_array[1][ver_count],
352 &verts_array[2][ver_count]);
353 for (int gg = 0; gg != num_nodes; ++gg, ++ver_count) {
354
355 mapGaussPts[gg] = start_vert_handle + ver_count;
356
357 auto set_float_precision = [](const double x) {
358 if (std::abs(x) < std::numeric_limits<float>::epsilon())
359 return 0.;
360 else
361 return x;
362 };
363
364 t_coords(i) = 0;
365 auto t_ele_coords = getFTensor1FromArray<3, 3>(coords);
366 for (int nn = 0; nn != CN::VerticesPerEntity(type); ++nn) {
367 t_coords(i) += t_n * t_ele_coords(i);
368 ++t_ele_coords;
369 ++t_n;
370 }
371
372 for (auto ii : {0, 1, 2})
373 t_coords(ii) = set_float_precision(t_coords(ii));
374
375 ++t_coords;
376 }
377
378 Tag th;
379 int def_in_the_loop = -1;
380 CHKERR getPostProcMesh().tag_get_handle(
381 "NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER, th, MB_TAG_CREAT | MB_TAG_SPARSE,
382 &def_in_the_loop);
383
384 postProcElements.clear();
385 const int num_el = level_ref_ele.size1();
386 const int num_nodes_on_ele = level_ref_ele.size2();
387 auto start_e = start_ele_handle + ele_count;
388 postProcElements = Range(start_e, start_e + num_el - 1);
389 for (auto tt = 0; tt != level_ref_ele.size1(); ++tt, ++ele_count) {
390 for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
391 conn[num_nodes_on_ele * ele_count + nn] =
392 mapGaussPts[level_ref_ele(tt, nn)];
393 }
394 }
395
396 const int n_in_the_loop = E::nInTheLoop;
397 CHKERR getPostProcMesh().tag_clear_data(th, postProcElements,
398 &n_in_the_loop);
399 CHKERR transferTags();
400
402 };
403
404 CHKERR set_gauss_pts(
405
406 ref_ele->levelGaussPtsOnRefMesh, ref_ele->levelRef,
407 ref_ele->levelShapeFunctions,
408
409 ref_ele->startingVertEleHandle, ref_ele->startingEleHandle,
410 ref_ele->verticesOnEleArrays, ref_ele->eleConn, ref_ele->countVertEle,
411 ref_ele->countEle
412
413 );
414
416};
417
418template <typename E>
421
422auto get_ref_ele = [&](const EntityType type) {
423 PostProcGenerateRefMeshPtr ref_ele_ptr;
424
425 auto it = refElementsMap.find(type);
426 if (it != refElementsMap.end()) {
427 ref_ele_ptr = it->second;
428 } else {
429 switch (type) {
430 case MBTET:
431 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTET>>();
432 break;
433 case MBHEX:
434 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBHEX>>();
435 break;
436 case MBTRI:
437 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
438 break;
439 case MBQUAD:
440 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBQUAD>>();
441 break;
442 case MBEDGE:
443 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBEDGE>>();
444 break;
445 default:
446 MOFEM_LOG("SELF", Sev::error)
447 << "Generation of reference elements for type < "
448 << moab::CN::EntityTypeName(type) << " > is not implemented";
449 CHK_THROW_MESSAGE(MOFEM_NOT_IMPLEMENTED, "Element not implemented");
450 }
451
452 CHK_THROW_MESSAGE(ref_ele_ptr->getOptions(optionsPrefix), "getOptions");
453 CHK_THROW_MESSAGE(ref_ele_ptr->generateReferenceElementMesh(),
454 "Error when generating reference element");
455
456 refElementsMap[type] = ref_ele_ptr;
457 }
458
459 return ref_ele_ptr;
460};
461
462 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
463
464 auto miit =
465 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
466 boost::make_tuple(this->getFEName(), this->getLoFERank()));
467 auto hi_miit =
468 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
469 boost::make_tuple(this->getFEName(), this->getHiFERank()));
470
471 const int number_of_ents_in_the_loop = this->getLoopSize();
472 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
473 SETERRQ(E::mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
474 "Wrong size of indicies. Inconsistent size number of iterated "
475 "elements iterated by problem and from range.");
476 }
477
478 for (auto &m : refElementsMap) {
479 m.second->nbVertices = 0;
480 m.second->nbEles = 0;
481 m.second->countEle = 0;
482 m.second->countVertEle = 0;
483 }
484
485 for (; miit != hi_miit; ++miit) {
486 auto type = (*miit)->getEntType();
487 auto ref_ele = get_ref_ele(type);
488
489 // Set pointer to element. So that getDataVectorDofsPtr in getMaxLevel
490 // can work
491 E::numeredEntFiniteElementPtr = *miit;
492 bool add = true;
493 if (E::exeTestHook) {
494 add = E::exeTestHook(this);
495 }
496
497 if (add) {
498 size_t level = getMaxLevel();
499 level = std::min(level, ref_ele->levelGaussPtsOnRefMesh.size() - 1);
500 ref_ele->nbVertices += ref_ele->levelGaussPtsOnRefMesh[level].size2();
501 ref_ele->nbEles += ref_ele->levelRef[level].size1();
502 }
503 }
504
505 auto alloc_vertices_and_elements_on_post_proc_mesh = [&]() {
507
508 ReadUtilIface *iface;
509 CHKERR getPostProcMesh().query_interface(iface);
510
511 for (auto &m : refElementsMap) {
512 if (m.second->nbEles) {
513 CHKERR iface->get_node_coords(3, m.second->nbVertices, 0,
514 m.second->startingVertEleHandle,
515 m.second->verticesOnEleArrays);
516 CHKERR iface->get_element_connect(
517 m.second->nbEles, m.second->levelRef[0].size2(), m.first, 0,
518 m.second->startingEleHandle, m.second->eleConn);
519
520 m.second->countEle = 0;
521 m.second->countVertEle = 0;
522 }
523 }
524
526 };
527
528 CHKERR alloc_vertices_and_elements_on_post_proc_mesh();
529
531}
532
533template <typename E>
536
537 auto update_elements = [&]() {
538 ReadUtilIface *iface;
539 CHKERR getPostProcMesh().query_interface(iface);
541
542 Range ents;
543 for (auto &m : refElementsMap) {
544 if (m.second->nbEles) {
545 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
546 << "Update < " << moab::CN::EntityTypeName(m.first)
547 << " number of processed " << m.second->countEle;
548 CHKERR iface->update_adjacencies(
549 m.second->startingEleHandle, m.second->countEle,
550 m.second->levelRef[0].size2(), m.second->eleConn);
551 ents.merge(Range(m.second->startingEleHandle,
552 m.second->startingEleHandle + m.second->countEle - 1));
553 }
554 }
555
557 };
558
559 auto remove_obsolete_entities = [&]() {
561 Range ents, adj;
562 for (auto &m : refElementsMap) {
563 if (m.second->nbEles) {
564 ents.merge(Range(m.second->startingEleHandle,
565 m.second->startingEleHandle + m.second->countEle - 1));
566 const int dim = moab::CN::Dimension(m.first);
567 for (auto d = 1; d != dim; ++d) {
568 CHKERR getPostProcMesh().get_adjacencies(ents, d, false, adj,
569 moab::Interface::UNION);
570 }
571 }
572 }
573 CHKERR getPostProcMesh().delete_entities(adj);
575 };
576
577 auto set_proc_tags = [&]() {
579 ParallelComm *pcomm_post_proc_mesh =
580 ParallelComm::get_pcomm(&(getPostProcMesh()), MYPCOMM_INDEX);
581 if (pcomm_post_proc_mesh) {
582 Range ents;
583 for (auto &m : refElementsMap) {
584 if (m.second->nbEles) {
585 ents.merge(
586 Range(m.second->startingEleHandle,
587 m.second->startingEleHandle + m.second->countEle - 1));
588 }
589 }
590
591 int rank = E::mField.get_comm_rank();
592 CHKERR getPostProcMesh().tag_clear_data(pcomm_post_proc_mesh->part_tag(),
593 ents, &rank);
594 }
596 };
597
598 CHKERR update_elements();
599 CHKERR remove_obsolete_entities();
600 CHKERR set_proc_tags();
601
603}
604
605template <typename E>
608
609 ParallelComm *pcomm_post_proc_mesh =
610 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
611 if (pcomm_post_proc_mesh != NULL)
612 delete pcomm_post_proc_mesh;
613 CHKERR getPostProcMesh().delete_mesh();
614 pcomm_post_proc_mesh =
615 new ParallelComm(&(getPostProcMesh()), PETSC_COMM_WORLD);
616
617 CHKERR preProcPostProc();
618
620}
621
622template <typename E>
625
626 ParallelComm *pcomm_post_proc_mesh =
627 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
628 if(pcomm_post_proc_mesh == nullptr)
629 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "PComm not allocated");
630
631 CHKERR postProcPostProc();
632 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
633
635}
636
638 return mapGaussPts;
639}
640
642 if (!coreMeshPtr)
643 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Core mesh not set");
644 moab::Interface &post_proc_mesh_interface = *coreMeshPtr;
645 return post_proc_mesh_interface;
646}
647
648template <typename E>
650 return postProcElements;
651}
652
653template <typename E>
655PostProcBrokenMeshInMoabBase<E>::writeFile(const std::string file_name) {
657 ParallelComm *pcomm_post_proc_mesh =
658 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
659 if (pcomm_post_proc_mesh == NULL)
660 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
661 "ParallelComm not allocated");
662 CHKERR getPostProcMesh().write_file(file_name.c_str(), "MOAB",
663 "PARALLEL=WRITE_PART");
665};
666
667template <typename E> struct PostProcBrokenMeshInMoab;
668
669template <>
671 : public PostProcBrokenMeshInMoabBase<VolumeElementForcesAndSourcesCore> {
673 VolumeElementForcesAndSourcesCore>::PostProcBrokenMeshInMoabBase;
674};
675
676template <>
678 : public PostProcBrokenMeshInMoabBase<FaceElementForcesAndSourcesCore> {
680 FaceElementForcesAndSourcesCore>::PostProcBrokenMeshInMoabBase;
681};
682
683template <>
685 : public PostProcBrokenMeshInMoabBase<EdgeElementForcesAndSourcesCore> {
687 EdgeElementForcesAndSourcesCore>::PostProcBrokenMeshInMoabBase;
688};
689
690/**
691 * @brief Post post-proc data at points from hash maps
692 *
693 * @tparam DIM1 dimension of vector in data_map_vec and first column of
694 * data_map_may
695 * @tparam DIM2 dimension of second column in data_map_mat
696 */
697template <int DIM1, int DIM2>
699
700 using DataMapVec = std::map<std::string, boost::shared_ptr<VectorDouble>>;
701 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
702
703 /**
704 * @brief Construct a new OpPostProcMapInMoab object
705 *
706 * @param post_proc_mesh postprocessing mesh
707 * @param map_gauss_pts map of gauss points to nodes of postprocessing mesh
708 * @param data_map_scalar hash map of scalar values (string is name of the
709 * tag)
710 * @param data_map_vec hash map of vector values
711 * @param data_map_mat hash map of second order tensor values
712 * @param data_symm_map_mat hash map of symmetric second order tensor values
713 */
714 OpPostProcMapInMoab(moab::Interface &post_proc_mesh,
715 std::vector<EntityHandle> &map_gauss_pts,
716 DataMapVec data_map_scalar, DataMapMat data_map_vec,
717 DataMapMat data_map_mat, DataMapMat data_symm_map_mat)
720 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
721 dataMapScalar(data_map_scalar), dataMapVec(data_map_vec),
722 dataMapMat(data_map_mat), dataMapSymmMat(data_symm_map_mat) {
723 // Operator is only executed for vertices
724 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
725 }
726 MoFEMErrorCode doWork(int side, EntityType type,
728
729private:
730 moab::Interface &postProcMesh;
731 std::vector<EntityHandle> &mapGaussPts;
736};
737
738template <int DIM1, int DIM2>
743
744 std::array<double, 9> def;
745 std::fill(def.begin(), def.end(), 0);
746
747 auto get_tag = [&](const std::string name, size_t size) {
748 Tag th;
749 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
750 MB_TAG_CREAT | MB_TAG_SPARSE,
751 def.data());
752 return th;
753 };
754
755 MatrixDouble3by3 mat(3, 3);
756
757 auto set_vector_3d = [&](auto &t) -> MatrixDouble3by3 & {
758 mat.clear();
759 for (size_t r = 0; r != DIM1; ++r)
760 mat(0, r) = t(r);
761 return mat;
762 };
763
764 auto set_matrix_3d = [&](auto &t) -> MatrixDouble3by3 & {
765 mat.clear();
766 for (size_t r = 0; r != DIM1; ++r)
767 for (size_t c = 0; c != DIM2; ++c)
768 mat(r, c) = t(r, c);
769 return mat;
770 };
771
772 auto set_matrix_symm_3d = [&](auto &t) -> MatrixDouble3by3 & {
773 mat.clear();
774 for (size_t r = 0; r != DIM1; ++r)
775 for (size_t c = 0; c != DIM1; ++c)
776 mat(r, c) = t(r, c);
777 return mat;
778 };
779
780 auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
781 mat.clear();
782 mat(0, 0) = t;
783 return mat;
784 };
785
786 auto set_float_precision = [](const double x) {
787 if (std::abs(x) < std::numeric_limits<float>::epsilon())
788 return 0.;
789 else
790 return x;
791 };
792
793 auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
794 for (auto &v : mat.data())
795 v = set_float_precision(v);
796 return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
797 &*mat.data().begin());
798 };
799
800 for (auto &m : dataMapScalar) {
801 if (m.second) {
802 auto th = get_tag(m.first, 1);
803 auto t_scl = getFTensor0FromVec(*m.second);
804 auto nb_integration_pts = getGaussPts().size2();
805 for (int gg = 0; gg != nb_integration_pts; ++gg) {
806 CHKERR set_tag(th, gg, set_scalar(t_scl));
807 ++t_scl;
808 }
809 }
810 }
811
812 for (auto &m : dataMapVec) {
813 if (m.second) {
814 auto th = get_tag(m.first, 3);
815 auto t_vec = getFTensor1FromMat<DIM1>(*m.second);
816 auto nb_integration_pts = getGaussPts().size2();
817 for (int gg = 0; gg != nb_integration_pts; ++gg) {
818 CHKERR set_tag(th, gg, set_vector_3d(t_vec));
819 ++t_vec;
820 }
821 }
822 }
823
824 for (auto &m : dataMapMat) {
825 if (m.second) {
826 auto th = get_tag(m.first, 9);
827 auto t_mat = getFTensor2FromMat<DIM1, DIM2>(*m.second);
828 auto nb_integration_pts = getGaussPts().size2();
829 for (int gg = 0; gg != nb_integration_pts; ++gg) {
830 CHKERR set_tag(th, gg, set_matrix_3d(t_mat));
831 ++t_mat;
832 }
833 }
834 }
835
836 for (auto &m : dataMapSymmMat) {
837 if (m.second) {
838 auto th = get_tag(m.first, 9);
839 auto t_mat = getFTensor2SymmetricFromMat<DIM1>(*m.second);
840 auto nb_integration_pts = getGaussPts().size2();
841 for (int gg = 0; gg != nb_integration_pts; ++gg) {
842 CHKERR set_tag(th, gg, set_matrix_symm_3d(t_mat));
843 ++t_mat;
844 }
845 }
846 }
847
849}
850
851template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseBeginImpl;
852
853template <typename E>
855 : protected PostProcBrokenMeshInMoabBase<E> {
856
858 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
859 std::string opts_prefix = "")
860 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr, opts_prefix) {}
861
864 ParallelComm *pcomm_post_proc_mesh =
865 ParallelComm::get_pcomm(&this->getPostProcMesh(), MYPCOMM_INDEX);
866 if (pcomm_post_proc_mesh != NULL)
867 delete pcomm_post_proc_mesh;
868 CHKERR this->getPostProcMesh().delete_mesh();
869 pcomm_post_proc_mesh =
870 new ParallelComm(&(this->getPostProcMesh()), PETSC_COMM_WORLD);
872 }
873
874 inline FEMethod *getFEMethod() { return this; }
875
878};
879
880template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseContImpl;
881
882template <typename E>
884 : public PostProcBrokenMeshInMoabBase<E> {
885
887 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
888 std::string opts_prefix = "")
889 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr, opts_prefix) {}
890
893
894protected:
896};
897
898template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseEndImpl;
899
900template <typename E>
902 : protected PostProcBrokenMeshInMoabBase<E> {
903
905 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
906 std::string opts_prefix = "")
907 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr, opts_prefix) {}
908
913 ParallelComm *pcomm_post_proc_mesh =
914 ParallelComm::get_pcomm(&this->getPostProcMesh(), MYPCOMM_INDEX);
915 if (pcomm_post_proc_mesh == nullptr)
916 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
917 "PComm not allocated");
918 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
920 }
921
922 inline FEMethod *getFEMethod() { return this; }
923
925};
926
927/**
928 * @brief Enable to run stack of post-processing elements. Use this to begin stack.
929 *
930 * See @ref scalar_check_approximation.cpp
931 *
932 */
935
936/**
937 * @brief Enable to run stack of post-processing elements.
938 *
939 * See @ref scalar_check_approximation.cpp
940 *
941 * @tparam E
942 */
943template <typename E>
946
947/**
948 * @brief Enable to run stack of post-processing elements. Use this to end stack.
949 *
950 */
953
954
955} // namespace MoFEM
956
957#endif //__POSTPROCBROKENMESHINMOABBASE_HPP
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#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
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ 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
constexpr int order
const int dim
FTensor::Index< 'm', SPACE_DIM > m
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1869
boost::shared_ptr< PostProcGenerateRefMeshBase > PostProcGenerateRefMeshPtr
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
constexpr double t
plate stiffness
Definition: plate.cpp:59
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
@ OPSPACE
operator do Work is execute on space data
structure to get information form mofem into EntitiesFieldData
Post post-proc data at points from hash maps.
DataMapMat dataMapMat
DataMapMat dataMapVec
DataMapMat dataMapSymmMat
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
moab::Interface & postProcMesh
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
std::vector< EntityHandle > & mapGaussPts
DataMapVec dataMapScalar
OpPostProcMapInMoab(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 OpPostProcMapInMoab object.
MoFEMErrorCode operator()()
MoFEMErrorCode preProcess()
MoFEMErrorCode postProcess()
FEMethod * getFEMethod()
PostProcBrokenMeshInMoabBaseBeginImpl(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr, std::string opts_prefix="")
PostProcBrokenMeshInMoabBaseContImpl(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr, std::string opts_prefix="")
MoFEMErrorCode postProcess()
MoFEMErrorCode preProcess()
MoFEMErrorCode operator()()
FEMethod * getFEMethod()
MoFEMErrorCode preProcess()
MoFEMErrorCode postProcess()
PostProcBrokenMeshInMoabBaseEndImpl(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr, std::string opts_prefix="")
auto & getPostProcElements()
Get postprocessing elements.
MoFEMErrorCode preProcess()
Generate vertices and elements.
virtual MoFEMErrorCode postProcPostProc()
MoFEMErrorCode setGaussPts(int order)
std::vector< Tag > tagsToTransfer
std::string optionsPrefix
Prefix for options.
PostProcBrokenMeshInMoabBase(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr, std::string opts_prefix="")
PostProcBrokenMeshInMoabBase(MoFEM::Interface &m_field, std::string opts_prefix="")
virtual MoFEMErrorCode transferTags()
auto & getMapGaussPts()
Get vector of vectors associated to integration points.
auto & getPostProcMesh()
Get postprocessing mesh.
MoFEMErrorCode postProcess()
boost::shared_ptr< moab::Core > coreMeshPtr
std::vector< EntityHandle > mapGaussPts
virtual int getMaxLevel() const
Determine refinement level based on fields approx ordre.
Range postProcElements
virtual ~PostProcBrokenMeshInMoabBase()
MoFEMErrorCode setTagsToTransfer(std::vector< Tag > tags_to_transfer)
Set tags to be transferred to post-processing mesh.
std::map< EntityType, PostProcGenerateRefMeshPtr > refElementsMap
int getRule(int order)
virtual MoFEMErrorCode preProcPostProc()
EntityHandle startingEleHandle
int countVertEle
virtual MoFEMErrorCode generateReferenceElementMesh()=0
virtual ~PostProcGenerateRefMeshBase()=default
EntityHandle * eleConn
std::vector< ublas::matrix< int > > levelRef
std::vector< double * > verticesOnEleArrays
PetscBool hoNodes
int defMaxLevel
virtual MoFEMErrorCode getOptions(std::string prefix)
std::vector< MatrixDouble > levelShapeFunctions
PostProcGenerateRefMeshBase()
int countEle
int nbVertices
std::vector< MatrixDouble > levelGaussPtsOnRefMesh
int nbEles
EntityHandle startingVertEleHandle
Element for postprocessing. Uses MoAB to generate post-processing mesh.