v0.15.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 <>
69
70template <>
75
76template <>
81
82template <>
87
88template <>
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 std::vector<double> tag_data_vec;
247 auto data_ptr = [&](auto tag) {
248 const void *tag_data;
249 auto core_ent = this->getFEEntityHandle();
250 CHK_MOAB_THROW(calc_mesh.tag_get_by_ptr(tag, &core_ent, 1, &tag_data),
251 "get tag data");
252 if (data_type(tag) == MB_TYPE_DOUBLE) {
253 tag_data_vec.resize(length(tag));
254 std::copy(static_cast<const double *>(tag_data),
255 static_cast<const double *>(tag_data) + length(tag),
256 tag_data_vec.begin());
257 for (auto &v : tag_data_vec)
258 v = std::abs(v) < std::numeric_limits<float>::epsilon() ? 0. : v;
259 return static_cast<const void *>(tag_data_vec.data());
260 }
261 return tag_data;
262 };
263
264 for (auto tag : tagsToTransfer) {
265 Tag tag_postproc;
266 CHKERR getPostProcMesh().tag_get_handle(
267 name(tag).c_str(), length(tag), data_type(tag), tag_postproc,
268 type(tag) | MB_TAG_CREAT, default_value(tag));
269 CHKERR getPostProcMesh().tag_clear_data(tag_postproc, postProcElements,
270 data_ptr(tag));
271 }
272
274};
275
276template <typename E>
278 MoFEM::Interface &m_field, std::string opts_prefix)
279 : E(m_field), optionsPrefix(opts_prefix) {}
280
281template <typename E>
283 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
284 std::string opts_prefix)
285 : PostProcBrokenMeshInMoabBase(m_field, opts_prefix) {
286 coreMeshPtr = core_mesh_ptr;
287}
288
289template <typename E>
291 ParallelComm *pcomm_post_proc_mesh =
292 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
293 if (pcomm_post_proc_mesh != NULL)
294 delete pcomm_post_proc_mesh;
295}
296
297template <typename E> int PostProcBrokenMeshInMoabBase<E>::getRule(int order) {
298 return -1;
299};
300
301template <typename E>
304
305 auto type = type_from_handle(this->getFEEntityHandle());
306
308
309 try {
310 ref_ele = refElementsMap.at(type);
311 } catch (const out_of_range &e) {
312 SETERRQ(
313 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
314 "Generation of reference elements for type <%s> is not implemented",
315 moab::CN::EntityTypeName(type));
316 }
317
318 auto set_gauss_pts = [&](auto &level_gauss_pts_on_ref_mesh, auto &level_ref,
319 auto &level_shape_functions,
320
321 auto start_vert_handle, auto start_ele_handle,
322 auto &verts_array, auto &conn, auto &ver_count,
323 auto &ele_count
324
325 ) {
327
328 size_t level = getMaxLevel();
329 level = std::min(level, level_gauss_pts_on_ref_mesh.size() - 1);
330
331 auto &level_ref_gauss_pts = level_gauss_pts_on_ref_mesh[level];
332 auto &level_ref_ele = level_ref[level];
333 auto &shape_functions = level_shape_functions[level];
334 E::gaussPts.resize(level_ref_gauss_pts.size1(), level_ref_gauss_pts.size2(),
335 false);
336 noalias(E::gaussPts) = level_ref_gauss_pts;
337
338 const auto fe_ent = E::numeredEntFiniteElementPtr->getEnt();
339 auto get_fe_coords = [&]() {
340 const EntityHandle *conn;
341 int num_nodes;
343 E::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes, true),
344 "error get connectivity");
345 VectorDouble coords(num_nodes * 3);
347 E::mField.get_moab().get_coords(conn, num_nodes, &*coords.begin()),
348 "error get coordinates");
349 return coords;
350 };
351
352 auto coords = get_fe_coords();
353
354 const int num_nodes = level_ref_gauss_pts.size2();
355 mapGaussPts.resize(level_ref_gauss_pts.size2());
356
357 FTensor::Index<'i', 3> i;
359 &*shape_functions.data().begin());
361 &verts_array[0][ver_count], &verts_array[1][ver_count],
362 &verts_array[2][ver_count]);
363 for (int gg = 0; gg != num_nodes; ++gg, ++ver_count) {
364
365 mapGaussPts[gg] = start_vert_handle + ver_count;
366
367 auto set_float_precision = [](const double x) {
368 if (std::abs(x) < std::numeric_limits<float>::epsilon())
369 return 0.;
370 else
371 return x;
372 };
373
374 t_coords(i) = 0;
375 auto t_ele_coords = getFTensor1FromArray<3, 3>(coords);
376 for (int nn = 0; nn != CN::VerticesPerEntity(type); ++nn) {
377 t_coords(i) += t_n * t_ele_coords(i);
378 ++t_ele_coords;
379 ++t_n;
380 }
381
382 for (auto ii : {0, 1, 2})
383 t_coords(ii) = set_float_precision(t_coords(ii));
384
385 ++t_coords;
386 }
387
388 Tag th;
389 int def_in_the_loop = -1;
390 CHKERR getPostProcMesh().tag_get_handle(
391 "NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER, th, MB_TAG_CREAT | MB_TAG_SPARSE,
392 &def_in_the_loop);
393
394 postProcElements.clear();
395 const int num_el = level_ref_ele.size1();
396 const int num_nodes_on_ele = level_ref_ele.size2();
397 auto start_e = start_ele_handle + ele_count;
398 postProcElements = Range(start_e, start_e + num_el - 1);
399 for (auto tt = 0; tt != level_ref_ele.size1(); ++tt, ++ele_count) {
400 for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
401 conn[num_nodes_on_ele * ele_count + nn] =
402 mapGaussPts[level_ref_ele(tt, nn)];
403 }
404 }
405
406 const int n_in_the_loop = E::nInTheLoop;
407 CHKERR getPostProcMesh().tag_clear_data(th, postProcElements,
408 &n_in_the_loop);
409 CHKERR transferTags();
410
412 };
413
414 CHKERR set_gauss_pts(
415
416 ref_ele->levelGaussPtsOnRefMesh, ref_ele->levelRef,
417 ref_ele->levelShapeFunctions,
418
419 ref_ele->startingVertEleHandle, ref_ele->startingEleHandle,
420 ref_ele->verticesOnEleArrays, ref_ele->eleConn, ref_ele->countVertEle,
421 ref_ele->countEle
422
423 );
424
426};
427
428template <typename E>
431
432auto get_ref_ele = [&](const EntityType type) {
433 PostProcGenerateRefMeshPtr ref_ele_ptr;
434
435 auto it = refElementsMap.find(type);
436 if (it != refElementsMap.end()) {
437 ref_ele_ptr = it->second;
438 } else {
439 switch (type) {
440 case MBTET:
441 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTET>>();
442 break;
443 case MBHEX:
444 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBHEX>>();
445 break;
446 case MBTRI:
447 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
448 break;
449 case MBQUAD:
450 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBQUAD>>();
451 break;
452 case MBEDGE:
453 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBEDGE>>();
454 break;
455 default:
456 MOFEM_LOG("SELF", Sev::error)
457 << "Generation of reference elements for type < "
458 << moab::CN::EntityTypeName(type) << " > is not implemented";
459 CHK_THROW_MESSAGE(MOFEM_NOT_IMPLEMENTED, "Element not implemented");
460 }
461
462 CHK_THROW_MESSAGE(ref_ele_ptr->getOptions(optionsPrefix), "getOptions");
463 CHK_THROW_MESSAGE(ref_ele_ptr->generateReferenceElementMesh(),
464 "Error when generating reference element");
465
466 refElementsMap[type] = ref_ele_ptr;
467 }
468
469 return ref_ele_ptr;
470};
471
472 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
473
474 auto miit =
475 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
476 boost::make_tuple(this->getFEName(), this->getLoFERank()));
477 auto hi_miit =
478 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
479 boost::make_tuple(this->getFEName(), this->getHiFERank()));
480
481 const int number_of_ents_in_the_loop = this->getLoopSize();
482 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
483 SETERRQ(E::mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
484 "Wrong size of indicies. Inconsistent size number of iterated "
485 "elements iterated by problem and from range.");
486 }
487
488 for (auto &m : refElementsMap) {
489 m.second->nbVertices = 0;
490 m.second->nbEles = 0;
491 m.second->countEle = 0;
492 m.second->countVertEle = 0;
493 }
494
495 for (; miit != hi_miit; ++miit) {
496 auto type = (*miit)->getEntType();
497 auto ref_ele = get_ref_ele(type);
498
499 // Set pointer to element. So that getDataVectorDofsPtr in getMaxLevel
500 // can work
501 E::numeredEntFiniteElementPtr = *miit;
502 bool add = true;
503 if (E::exeTestHook) {
504 add = E::exeTestHook(this);
505 }
506
507 if (add) {
508 size_t level = getMaxLevel();
509 level = std::min(level, ref_ele->levelGaussPtsOnRefMesh.size() - 1);
510 ref_ele->nbVertices += ref_ele->levelGaussPtsOnRefMesh[level].size2();
511 ref_ele->nbEles += ref_ele->levelRef[level].size1();
512 }
513 }
514
515 auto alloc_vertices_and_elements_on_post_proc_mesh = [&]() {
517
518 ReadUtilIface *iface;
519 CHKERR getPostProcMesh().query_interface(iface);
520
521 for (auto &m : refElementsMap) {
522 if (m.second->nbEles) {
523 CHKERR iface->get_node_coords(3, m.second->nbVertices, 0,
524 m.second->startingVertEleHandle,
525 m.second->verticesOnEleArrays);
526 CHKERR iface->get_element_connect(
527 m.second->nbEles, m.second->levelRef[0].size2(), m.first, 0,
528 m.second->startingEleHandle, m.second->eleConn);
529
530 m.second->countEle = 0;
531 m.second->countVertEle = 0;
532 }
533 }
534
536 };
537
538 CHKERR alloc_vertices_and_elements_on_post_proc_mesh();
539
541}
542
543template <typename E>
546
547 auto update_elements = [&]() {
548 ReadUtilIface *iface;
549 CHKERR getPostProcMesh().query_interface(iface);
551
552 Range ents;
553 for (auto &m : refElementsMap) {
554 if (m.second->nbEles) {
555 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
556 << "Update < " << moab::CN::EntityTypeName(m.first)
557 << " number of processed " << m.second->countEle;
558 CHKERR iface->update_adjacencies(
559 m.second->startingEleHandle, m.second->countEle,
560 m.second->levelRef[0].size2(), m.second->eleConn);
561 ents.merge(Range(m.second->startingEleHandle,
562 m.second->startingEleHandle + m.second->countEle - 1));
563 }
564 }
565
567 };
568
569 auto remove_obsolete_entities = [&]() {
571 Range ents, adj;
572 for (auto &m : refElementsMap) {
573 if (m.second->nbEles) {
574 ents.merge(Range(m.second->startingEleHandle,
575 m.second->startingEleHandle + m.second->countEle - 1));
576 const int dim = moab::CN::Dimension(m.first);
577 for (auto d = 1; d != dim; ++d) {
578 CHKERR getPostProcMesh().get_adjacencies(ents, d, false, adj,
579 moab::Interface::UNION);
580 }
581 }
582 }
583 CHKERR getPostProcMesh().delete_entities(adj);
585 };
586
587 auto set_proc_tags = [&]() {
589 ParallelComm *pcomm_post_proc_mesh =
590 ParallelComm::get_pcomm(&(getPostProcMesh()), MYPCOMM_INDEX);
591 if (pcomm_post_proc_mesh) {
592 Range ents;
593 for (auto &m : refElementsMap) {
594 if (m.second->nbEles) {
595 ents.merge(
596 Range(m.second->startingEleHandle,
597 m.second->startingEleHandle + m.second->countEle - 1));
598 }
599 }
600
601 int rank = E::mField.get_comm_rank();
602 CHKERR getPostProcMesh().tag_clear_data(pcomm_post_proc_mesh->part_tag(),
603 ents, &rank);
604 }
606 };
607
608 CHKERR update_elements();
609 CHKERR remove_obsolete_entities();
610 CHKERR set_proc_tags();
611
613}
614
615template <typename E>
618
619 ParallelComm *pcomm_post_proc_mesh =
620 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
621 if (pcomm_post_proc_mesh != NULL)
622 delete pcomm_post_proc_mesh;
623 CHKERR getPostProcMesh().delete_mesh();
624 pcomm_post_proc_mesh =
625 new ParallelComm(&(getPostProcMesh()), PETSC_COMM_WORLD);
626
627 CHKERR preProcPostProc();
628
630}
631
632template <typename E>
635
636 ParallelComm *pcomm_post_proc_mesh =
637 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
638 if(pcomm_post_proc_mesh == nullptr)
639 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "PComm not allocated");
640
641 CHKERR postProcPostProc();
642 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
643
645}
646
648 return mapGaussPts;
649}
650
652 if (!coreMeshPtr)
653 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Core mesh not set");
654 moab::Interface &post_proc_mesh_interface = *coreMeshPtr;
655 return post_proc_mesh_interface;
656}
657
658template <typename E>
660 return postProcElements;
661}
662
663template <typename E>
665PostProcBrokenMeshInMoabBase<E>::writeFile(const std::string file_name) {
667 ParallelComm *pcomm_post_proc_mesh =
668 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
669 if (pcomm_post_proc_mesh == NULL)
670 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
671 "ParallelComm not allocated");
672 CHKERR getPostProcMesh().write_file(file_name.c_str(), "MOAB",
673 "PARALLEL=WRITE_PART");
675};
676
677template <typename E> struct PostProcBrokenMeshInMoab;
678
679template <>
685
686template <>
692
693template <>
699
700/**
701 * @brief Post post-proc data at points from hash maps
702 *
703 * @tparam DIM1 dimension of vector in data_map_vec and first column of
704 * data_map_may
705 * @tparam DIM2 dimension of second column in data_map_mat
706 */
707template <int DIM1, int DIM2,
709struct OpPostProcMapInMoab : public O {
710
711 using DataMapVec = std::map<std::string, boost::shared_ptr<VectorDouble>>;
712 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
713
714 /**
715 * @brief Construct a new OpPostProcMapInMoab object
716 *
717 * @param post_proc_mesh postprocessing mesh
718 * @param map_gauss_pts map of gauss points to nodes of postprocessing mesh
719 * @param data_map_scalar hash map of scalar values (string is name of the
720 * tag)
721 * @param data_map_vec hash map of vector values
722 * @param data_map_mat hash map of second order tensor values
723 * @param data_symm_map_mat hash map of symmetric second order tensor values
724 */
725 OpPostProcMapInMoab(moab::Interface &post_proc_mesh,
726 std::vector<EntityHandle> &map_gauss_pts,
727 DataMapVec data_map_scalar, DataMapMat data_map_vec,
728 DataMapMat data_map_mat, DataMapMat data_symm_map_mat)
729 : O(NOSPACE, O::OPSPACE), postProcMesh(post_proc_mesh),
730 mapGaussPts(map_gauss_pts), dataMapScalar(data_map_scalar),
731 dataMapVec(data_map_vec), dataMapMat(data_map_mat),
732 dataMapSymmMat(data_symm_map_mat) {
733 // Operator is only executed for vertices
734 std::fill(&O::doEntities[MBEDGE], &O::doEntities[MBMAXTYPE], false);
735 }
736 MoFEMErrorCode doWork(int side, EntityType type,
738
739private:
740 moab::Interface &postProcMesh;
741 std::vector<EntityHandle> &mapGaussPts;
746};
747
748template <int DIM1, int DIM2, typename O>
753
754 std::array<double, 9> def;
755 std::fill(def.begin(), def.end(), 0);
756
757 auto get_tag = [&](const std::string name, size_t size) {
758 Tag th;
759 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
760 MB_TAG_CREAT | MB_TAG_SPARSE,
761 def.data());
762 return th;
763 };
764
765 MatrixDouble3by3 mat(3, 3);
766
767 auto set_vector_3d = [&](auto &t) -> MatrixDouble3by3 & {
768 mat.clear();
769 for (size_t r = 0; r != DIM1; ++r)
770 mat(0, r) = t(r);
771 return mat;
772 };
773
774 auto set_matrix_3d = [&](auto &t) -> MatrixDouble3by3 & {
775 mat.clear();
776 for (size_t r = 0; r != DIM1; ++r)
777 for (size_t c = 0; c != DIM2; ++c)
778 mat(r, c) = t(r, c);
779 return mat;
780 };
781
782 auto set_matrix_symm_3d = [&](auto &t) -> MatrixDouble3by3 & {
783 mat.clear();
784 for (size_t r = 0; r != DIM1; ++r)
785 for (size_t c = 0; c != DIM1; ++c)
786 mat(r, c) = t(r, c);
787 return mat;
788 };
789
790 auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
791 mat.clear();
792 mat(0, 0) = t;
793 return mat;
794 };
795
796 auto set_float_precision = [](const double x) {
797 if (std::abs(x) < std::numeric_limits<float>::epsilon())
798 return 0.;
799 else
800 return x;
801 };
802
803 auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
804 for (auto &v : mat.data())
805 v = set_float_precision(v);
806 return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
807 &*mat.data().begin());
808 };
809
810 for (auto &m : dataMapScalar) {
811 if (m.second) {
812 auto th = get_tag(m.first, 1);
813 auto t_scl = getFTensor0FromVec(*m.second);
814 auto nb_integration_pts = O::getGaussPts().size2();
815 for (int gg = 0; gg != nb_integration_pts; ++gg) {
816 CHKERR set_tag(th, gg, set_scalar(t_scl));
817 ++t_scl;
818 }
819 }
820 }
821
822 for (auto &m : dataMapVec) {
823 if (m.second) {
824 auto th = get_tag(m.first, 3);
825 auto t_vec = getFTensor1FromMat<DIM1>(*m.second);
826 auto nb_integration_pts = O::getGaussPts().size2();
827 for (int gg = 0; gg != nb_integration_pts; ++gg) {
828 CHKERR set_tag(th, gg, set_vector_3d(t_vec));
829 ++t_vec;
830 }
831 }
832 }
833
834 for (auto &m : dataMapMat) {
835 if (m.second) {
836 auto th = get_tag(m.first, 9);
837 auto t_mat = getFTensor2FromMat<DIM1, DIM2>(*m.second);
838 auto nb_integration_pts = O::getGaussPts().size2();
839 for (int gg = 0; gg != nb_integration_pts; ++gg) {
840 CHKERR set_tag(th, gg, set_matrix_3d(t_mat));
841 ++t_mat;
842 }
843 }
844 }
845
846 for (auto &m : dataMapSymmMat) {
847 if (m.second) {
848 auto th = get_tag(m.first, 9);
849 auto t_mat = getFTensor2SymmetricFromMat<DIM1>(*m.second);
850 auto nb_integration_pts = O::getGaussPts().size2();
851 for (int gg = 0; gg != nb_integration_pts; ++gg) {
852 CHKERR set_tag(th, gg, set_matrix_symm_3d(t_mat));
853 ++t_mat;
854 }
855 }
856 }
857
859}
860
861template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseBeginImpl;
862
863template <typename E>
865 : protected PostProcBrokenMeshInMoabBase<E> {
866
868 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
869 std::string opts_prefix = "")
870 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr, opts_prefix) {}
871
874 ParallelComm *pcomm_post_proc_mesh =
875 ParallelComm::get_pcomm(&this->getPostProcMesh(), MYPCOMM_INDEX);
876 if (pcomm_post_proc_mesh != NULL)
877 delete pcomm_post_proc_mesh;
878 CHKERR this->getPostProcMesh().delete_mesh();
879 pcomm_post_proc_mesh =
880 new ParallelComm(&(this->getPostProcMesh()), PETSC_COMM_WORLD);
882 }
883
884 inline FEMethod *getFEMethod() { return this; }
885
888};
889
890template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseContImpl;
891
892template <typename E>
894 : public PostProcBrokenMeshInMoabBase<E> {
895
897 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
898 std::string opts_prefix = "")
899 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr, opts_prefix) {}
900
903
904protected:
906};
907
908template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseEndImpl;
909
910template <typename E>
912 : protected PostProcBrokenMeshInMoabBase<E> {
913
915 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr,
916 std::string opts_prefix = "")
917 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr, opts_prefix) {}
918
923 ParallelComm *pcomm_post_proc_mesh =
924 ParallelComm::get_pcomm(&this->getPostProcMesh(), MYPCOMM_INDEX);
925 if (pcomm_post_proc_mesh == nullptr)
926 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
927 "PComm not allocated");
928 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
930 }
931
932 inline FEMethod *getFEMethod() { return this; }
933
935};
936
937/**
938 * @brief Enable to run stack of post-processing elements. Use this to begin stack.
939 *
940 * See @ref scalar_check_approximation.cpp
941 *
942 */
945
946/**
947 * @brief Enable to run stack of post-processing elements.
948 *
949 * See @ref scalar_check_approximation.cpp
950 *
951 * @tparam E
952 */
953template <typename E>
956
957/**
958 * @brief Enable to run stack of post-processing elements. Use this to end stack.
959 *
960 */
963
964
965} // namespace MoFEM
966
967#endif //__POSTPROCBROKENMESHINMOABBASE_HPP
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOSPACE
Definition definitions.h:83
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
#define MOFEM_LOG(channel, severity)
Log.
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.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
auto type_from_handle(const EntityHandle h)
get type from entity handle
boost::shared_ptr< PostProcGenerateRefMeshBase > PostProcGenerateRefMeshPtr
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
constexpr double t
plate stiffness
Definition plate.cpp:58
FTensor::Index< 'm', 3 > m
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
Post post-proc data at points from hash maps.
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
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.
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
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="")
PostProcBrokenMeshInMoabBaseEndImpl(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="")
auto & getMapGaussPts()
Get vector of vectors associated to integration points.
MoFEMErrorCode setTagsToTransfer(std::vector< Tag > tags_to_transfer)
Set tags to be transferred to post-processing mesh.
PostProcBrokenMeshInMoabBase(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr, std::string opts_prefix="")
auto & getPostProcMesh()
Get postprocessing mesh.
std::map< EntityType, PostProcGenerateRefMeshPtr > refElementsMap
auto & getPostProcElements()
Get postprocessing elements.
virtual int getMaxLevel() const
Determine refinement level based on fields approx ordre.
MoFEMErrorCode preProcess()
Generate vertices and elements.
virtual MoFEMErrorCode generateReferenceElementMesh()=0
virtual ~PostProcGenerateRefMeshBase()=default
std::vector< ublas::matrix< int > > levelRef
virtual MoFEMErrorCode getOptions(std::string prefix)
Element for postprocessing. Uses MoAB to generate post-processing mesh.