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
422 auto get_ref_ele = [&](const EntityType type) {
423 PostProcGenerateRefMeshPtr ref_ele_ptr;
424
425 try {
426
427 ref_ele_ptr = refElementsMap.at(type);
428
429 } catch (const out_of_range &e) {
430
431 switch (type) {
432 case MBTET:
433 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTET>>();
434 break;
435 case MBHEX:
436 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBHEX>>();
437 break;
438 case MBTRI:
439 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
440 break;
441 case MBQUAD:
442 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBQUAD>>();
443 break;
444 case MBEDGE:
445 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBEDGE>>();
446 break;
447 default:
448 MOFEM_LOG("SELF", Sev::error)
449 << "Generation of reference elements for type < "
450 << moab::CN::EntityTypeName(type) << " > is not implemented";
451 CHK_THROW_MESSAGE(MOFEM_NOT_IMPLEMENTED, "Element not implemented");
452 }
453
454 CHK_THROW_MESSAGE(ref_ele_ptr->getOptions(optionsPrefix), "getOptions");
455 CHK_THROW_MESSAGE(ref_ele_ptr->generateReferenceElementMesh(),
456 "Error when generating reference element");
457 }
458
459 refElementsMap[type] = ref_ele_ptr;
460
461 return ref_ele_ptr;
462 };
463
464 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
465
466 auto miit =
467 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
468 boost::make_tuple(this->getFEName(), this->getLoFERank()));
469 auto hi_miit =
470 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
471 boost::make_tuple(this->getFEName(), this->getHiFERank()));
472
473 const int number_of_ents_in_the_loop = this->getLoopSize();
474 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
475 SETERRQ(E::mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
476 "Wrong size of indicies. Inconsistent size number of iterated "
477 "elements iterated by problem and from range.");
478 }
479
480 for (auto &m : refElementsMap) {
481 m.second->nbVertices = 0;
482 m.second->nbEles = 0;
483 m.second->countEle = 0;
484 m.second->countVertEle = 0;
485 }
486
487 for (; miit != hi_miit; ++miit) {
488 auto type = (*miit)->getEntType();
489 auto ref_ele = get_ref_ele(type);
490
491 // Set pointer to element. So that getDataVectorDofsPtr in getMaxLevel
492 // can work
493 E::numeredEntFiniteElementPtr = *miit;
494 bool add = true;
495 if (E::exeTestHook) {
496 add = E::exeTestHook(this);
497 }
498
499 if (add) {
500 size_t level = getMaxLevel();
501 level = std::min(level, ref_ele->levelGaussPtsOnRefMesh.size() - 1);
502 ref_ele->nbVertices += ref_ele->levelGaussPtsOnRefMesh[level].size2();
503 ref_ele->nbEles += ref_ele->levelRef[level].size1();
504 }
505 }
506
507 auto alloc_vertices_and_elements_on_post_proc_mesh = [&]() {
509
510 ReadUtilIface *iface;
511 CHKERR getPostProcMesh().query_interface(iface);
512
513 for (auto &m : refElementsMap) {
514 if (m.second->nbEles) {
515 CHKERR iface->get_node_coords(3, m.second->nbVertices, 0,
516 m.second->startingVertEleHandle,
517 m.second->verticesOnEleArrays);
518 CHKERR iface->get_element_connect(
519 m.second->nbEles, m.second->levelRef[0].size2(), m.first, 0,
520 m.second->startingEleHandle, m.second->eleConn);
521
522 m.second->countEle = 0;
523 m.second->countVertEle = 0;
524 }
525 }
526
528 };
529
530 CHKERR alloc_vertices_and_elements_on_post_proc_mesh();
531
533}
534
535template <typename E>
538
539 auto update_elements = [&]() {
540 ReadUtilIface *iface;
541 CHKERR getPostProcMesh().query_interface(iface);
543
544 Range ents;
545 for (auto &m : refElementsMap) {
546 if (m.second->nbEles) {
547 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
548 << "Update < " << moab::CN::EntityTypeName(m.first)
549 << " number of processed " << m.second->countEle;
550 CHKERR iface->update_adjacencies(
551 m.second->startingEleHandle, m.second->countEle,
552 m.second->levelRef[0].size2(), m.second->eleConn);
553 ents.merge(Range(m.second->startingEleHandle,
554 m.second->startingEleHandle + m.second->countEle - 1));
555 }
556 }
557
559 };
560
561 auto remove_obsolete_entities = [&]() {
563 Range ents, adj;
564 for (auto &m : refElementsMap) {
565 if (m.second->nbEles) {
566 ents.merge(Range(m.second->startingEleHandle,
567 m.second->startingEleHandle + m.second->countEle - 1));
568 const int dim = moab::CN::Dimension(m.first);
569 for (auto d = 1; d != dim; ++d) {
570 CHKERR getPostProcMesh().get_adjacencies(ents, d, false, adj,
571 moab::Interface::UNION);
572 }
573 }
574 }
575 CHKERR getPostProcMesh().delete_entities(adj);
577 };
578
579 auto set_proc_tags = [&]() {
581 ParallelComm *pcomm_post_proc_mesh =
582 ParallelComm::get_pcomm(&(getPostProcMesh()), MYPCOMM_INDEX);
583 if (pcomm_post_proc_mesh) {
584 Range ents;
585 for (auto &m : refElementsMap) {
586 if (m.second->nbEles) {
587 ents.merge(
588 Range(m.second->startingEleHandle,
589 m.second->startingEleHandle + m.second->countEle - 1));
590 }
591 }
592
593 int rank = E::mField.get_comm_rank();
594 CHKERR getPostProcMesh().tag_clear_data(pcomm_post_proc_mesh->part_tag(),
595 ents, &rank);
596 }
598 };
599
600 CHKERR update_elements();
601 CHKERR remove_obsolete_entities();
602 CHKERR set_proc_tags();
603
605}
606
607template <typename E>
610
611 ParallelComm *pcomm_post_proc_mesh =
612 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
613 if (pcomm_post_proc_mesh != NULL)
614 delete pcomm_post_proc_mesh;
615 CHKERR getPostProcMesh().delete_mesh();
616 pcomm_post_proc_mesh =
617 new ParallelComm(&(getPostProcMesh()), PETSC_COMM_WORLD);
618
619 CHKERR preProcPostProc();
620
622}
623
624template <typename E>
627
628 ParallelComm *pcomm_post_proc_mesh =
629 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
630 if(pcomm_post_proc_mesh == nullptr)
631 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "PComm not allocated");
632
633 CHKERR postProcPostProc();
634 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
635
637}
638
640 return mapGaussPts;
641}
642
644 if (!coreMeshPtr)
645 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Core mesh not set");
646 moab::Interface &post_proc_mesh_interface = *coreMeshPtr;
647 return post_proc_mesh_interface;
648}
649
650template <typename E>
652 return postProcElements;
653}
654
655template <typename E>
657PostProcBrokenMeshInMoabBase<E>::writeFile(const std::string file_name) {
659 ParallelComm *pcomm_post_proc_mesh =
660 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
661 if (pcomm_post_proc_mesh == NULL)
662 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
663 "ParallelComm not allocated");
664 CHKERR getPostProcMesh().write_file(file_name.c_str(), "MOAB",
665 "PARALLEL=WRITE_PART");
667};
668
669template <typename E> struct PostProcBrokenMeshInMoab;
670
671template <>
673 : public PostProcBrokenMeshInMoabBase<VolumeElementForcesAndSourcesCore> {
675 VolumeElementForcesAndSourcesCore>::PostProcBrokenMeshInMoabBase;
676};
677
678template <>
680 : public PostProcBrokenMeshInMoabBase<FaceElementForcesAndSourcesCore> {
682 FaceElementForcesAndSourcesCore>::PostProcBrokenMeshInMoabBase;
683};
684
685template <>
687 : public PostProcBrokenMeshInMoabBase<EdgeElementForcesAndSourcesCore> {
689 EdgeElementForcesAndSourcesCore>::PostProcBrokenMeshInMoabBase;
690};
691
692/**
693 * @brief Post post-proc data at points from hash maps
694 *
695 * @tparam DIM1 dimension of vector in data_map_vec and first column of
696 * data_map_may
697 * @tparam DIM2 dimension of second column in data_map_mat
698 */
699template <int DIM1, int DIM2>
701
702 using DataMapVec = std::map<std::string, boost::shared_ptr<VectorDouble>>;
703 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
704
705 /**
706 * @brief Construct a new OpPostProcMapInMoab object
707 *
708 * @param post_proc_mesh postprocessing mesh
709 * @param map_gauss_pts map of gauss points to nodes of postprocessing mesh
710 * @param data_map_scalar hash map of scalar values (string is name of the
711 * tag)
712 * @param data_map_vec hash map of vector values
713 * @param data_map_mat hash map of second order tensor values
714 * @param data_symm_map_mat hash map of symmetric second order tensor values
715 */
716 OpPostProcMapInMoab(moab::Interface &post_proc_mesh,
717 std::vector<EntityHandle> &map_gauss_pts,
718 DataMapVec data_map_scalar, DataMapMat data_map_vec,
719 DataMapMat data_map_mat, DataMapMat data_symm_map_mat)
722 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
723 dataMapScalar(data_map_scalar), dataMapVec(data_map_vec),
724 dataMapMat(data_map_mat), dataMapSymmMat(data_symm_map_mat) {
725 // Operator is only executed for vertices
726 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
727 }
728 MoFEMErrorCode doWork(int side, EntityType type,
730
731private:
732 moab::Interface &postProcMesh;
733 std::vector<EntityHandle> &mapGaussPts;
738};
739
740template <int DIM1, int DIM2>
745
746 std::array<double, 9> def;
747 std::fill(def.begin(), def.end(), 0);
748
749 auto get_tag = [&](const std::string name, size_t size) {
750 Tag th;
751 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
752 MB_TAG_CREAT | MB_TAG_SPARSE,
753 def.data());
754 return th;
755 };
756
757 MatrixDouble3by3 mat(3, 3);
758
759 auto set_vector_3d = [&](auto &t) -> MatrixDouble3by3 & {
760 mat.clear();
761 for (size_t r = 0; r != DIM1; ++r)
762 mat(0, r) = t(r);
763 return mat;
764 };
765
766 auto set_matrix_3d = [&](auto &t) -> MatrixDouble3by3 & {
767 mat.clear();
768 for (size_t r = 0; r != DIM1; ++r)
769 for (size_t c = 0; c != DIM2; ++c)
770 mat(r, c) = t(r, c);
771 return mat;
772 };
773
774 auto set_matrix_symm_3d = [&](auto &t) -> MatrixDouble3by3 & {
775 mat.clear();
776 for (size_t r = 0; r != DIM1; ++r)
777 for (size_t c = 0; c != DIM1; ++c)
778 mat(r, c) = t(r, c);
779 return mat;
780 };
781
782 auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
783 mat.clear();
784 mat(0, 0) = t;
785 return mat;
786 };
787
788 auto set_float_precision = [](const double x) {
789 if (std::abs(x) < std::numeric_limits<float>::epsilon())
790 return 0.;
791 else
792 return x;
793 };
794
795 auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
796 for (auto &v : mat.data())
797 v = set_float_precision(v);
798 return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
799 &*mat.data().begin());
800 };
801
802 for (auto &m : dataMapScalar) {
803 if (m.second) {
804 auto th = get_tag(m.first, 1);
805 auto t_scl = getFTensor0FromVec(*m.second);
806 auto nb_integration_pts = getGaussPts().size2();
807 for (int gg = 0; gg != nb_integration_pts; ++gg) {
808 CHKERR set_tag(th, gg, set_scalar(t_scl));
809 ++t_scl;
810 }
811 }
812 }
813
814 for (auto &m : dataMapVec) {
815 if (m.second) {
816 auto th = get_tag(m.first, 3);
817 auto t_vec = getFTensor1FromMat<DIM1>(*m.second);
818 auto nb_integration_pts = getGaussPts().size2();
819 for (int gg = 0; gg != nb_integration_pts; ++gg) {
820 CHKERR set_tag(th, gg, set_vector_3d(t_vec));
821 ++t_vec;
822 }
823 }
824 }
825
826 for (auto &m : dataMapMat) {
827 if (m.second) {
828 auto th = get_tag(m.first, 9);
829 auto t_mat = getFTensor2FromMat<DIM1, DIM2>(*m.second);
830 auto nb_integration_pts = getGaussPts().size2();
831 for (int gg = 0; gg != nb_integration_pts; ++gg) {
832 CHKERR set_tag(th, gg, set_matrix_3d(t_mat));
833 ++t_mat;
834 }
835 }
836 }
837
838 for (auto &m : dataMapSymmMat) {
839 if (m.second) {
840 auto th = get_tag(m.first, 9);
841 auto t_mat = getFTensor2SymmetricFromMat<DIM1>(*m.second);
842 auto nb_integration_pts = getGaussPts().size2();
843 for (int gg = 0; gg != nb_integration_pts; ++gg) {
844 CHKERR set_tag(th, gg, set_matrix_symm_3d(t_mat));
845 ++t_mat;
846 }
847 }
848 }
849
851}
852
853template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseBeginImpl;
854
855template <typename E>
857 : protected PostProcBrokenMeshInMoabBase<E> {
858
860 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
861 std::string opts_prefix = "")
862 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr, opts_prefix) {}
863
866 ParallelComm *pcomm_post_proc_mesh =
867 ParallelComm::get_pcomm(&this->getPostProcMesh(), MYPCOMM_INDEX);
868 if (pcomm_post_proc_mesh != NULL)
869 delete pcomm_post_proc_mesh;
870 CHKERR this->getPostProcMesh().delete_mesh();
871 pcomm_post_proc_mesh =
872 new ParallelComm(&(this->getPostProcMesh()), PETSC_COMM_WORLD);
874 }
875
876 inline FEMethod *getFEMethod() { return this; }
877
880};
881
882template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseContImpl;
883
884template <typename E>
886 : public PostProcBrokenMeshInMoabBase<E> {
887
889 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
890 std::string opts_prefix = "")
891 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr, opts_prefix) {}
892
895
896protected:
898};
899
900template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseEndImpl;
901
902template <typename E>
904 : protected PostProcBrokenMeshInMoabBase<E> {
905
907 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
908 std::string opts_prefix = "")
909 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr, opts_prefix) {}
910
915 ParallelComm *pcomm_post_proc_mesh =
916 ParallelComm::get_pcomm(&this->getPostProcMesh(), MYPCOMM_INDEX);
917 if (pcomm_post_proc_mesh == nullptr)
918 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
919 "PComm not allocated");
920 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
922 }
923
924 inline FEMethod *getFEMethod() { return this; }
925
927};
928
929/**
930 * @brief Enable to run stack of post-processing elements. Use this to begin stack.
931 *
932 * See @ref scalar_check_approximation.cpp
933 *
934 */
937
938/**
939 * @brief Enable to run stack of post-processing elements.
940 *
941 * See @ref scalar_check_approximation.cpp
942 *
943 * @tparam E
944 */
945template <typename E>
948
949/**
950 * @brief Enable to run stack of post-processing elements. Use this to end stack.
951 *
952 */
955
956
957} // namespace MoFEM
958
959#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)
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:1779
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.