v0.13.2
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 MoFEMErrorCode getOptions(); //< get options from command
46
48
49 PetscBool hoNodes; //< if true mid nodes are added
50 int defMaxLevel; //< default max number of refinement levels
51 std::string optPrefix; //< prefix for options
52};
53
55 boost::shared_ptr<PostProcGenerateRefMeshBase>;
56
57/**
58 * @brief Element for postprocessing. Uses MoAB to generate post-processing
59 * mesh.
60 *
61 * @tparam T Finite Element Implementation
62 */
63template <EntityType T> struct PostProcGenerateRefMesh;
64
65template <>
68 MoFEMErrorCode generateReferenceElementMesh();
69};
70
71template <>
74 MoFEMErrorCode generateReferenceElementMesh();
75};
76
77template <>
80 MoFEMErrorCode generateReferenceElementMesh();
81};
82
83template <>
86 MoFEMErrorCode generateReferenceElementMesh();
87};
88
89template <>
92 MoFEMErrorCode generateReferenceElementMesh();
93};
94
95template <typename E> struct PostProcBrokenMeshInMoabBase : public E {
96
98
100
101 /**
102 * @brief Get vector of vectors associated to integration points
103 *
104 * @return std::vector<EntityHandle>&
105 */
106 inline auto &getMapGaussPts();
107
108 /**
109 * @brief Get postprocessing mesh
110 *
111 * @return moab::Interface&
112 */
113 inline auto &getPostProcMesh();
114
115 /**
116 * @brief Get postprocessing elements
117 *
118 * @return auto&
119 */
120 inline auto &getPostProcElements();
121
122 /**
123 * \brief wrote results in (MOAB) format, use "file_name.h5m"
124 * @param file_name file name (should always end with .h5m)
125 * @return error code
126 * \ingroup mofem_fs_post_proc
127 */
128 MoFEMErrorCode writeFile(const std::string file_name);
129
130 /**
131 * @brief Set tags to be transferred to post-processing mesh
132 *
133 * @param tags_to_transfer
134 * @return MoFEMErrorCode
135 */
136 MoFEMErrorCode setTagsToTransfer(std::vector<Tag> tags_to_transfer);
137
138protected:
139
141
142 /**
143 * @brief Generate vertices and elements
144 *
145 * @return MoFEMErrorCode
146 */
148
150
151 int getRule(int order);
152
153 /**
154 * @brief Determine refinement level based on fields approx ordre.
155 *
156 * level = (order - 1) / 2
157 *
158 * @return int
159 */
160 virtual int getMaxLevel() const;
161
162 boost::shared_ptr<moab::Core> coreMeshPtr = boost::make_shared<moab::Core>();
164 boost::shared_ptr<moab::Core> core_mesh_ptr);
165
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
179
180 virtual MoFEMErrorCode transferTags(); ///< transfer tags from mesh to
181 ///< post-process mesh
182};
183
184template <typename E>
186 std::vector<Tag> tags_to_transfer) {
188 tagsToTransfer.swap(tags_to_transfer);
190}
191
192template <typename E> int PostProcBrokenMeshInMoabBase<E>::getMaxLevel() const {
193 auto get_element_max_dofs_order = [&]() {
194 int max_order = 0;
195 auto dofs_vec = E::getDataVectorDofsPtr();
196 for (auto &dof : *dofs_vec) {
197 const int dof_order = dof->getDofOrder();
198 max_order = (max_order < dof_order) ? dof_order : max_order;
199 };
200 return max_order;
201 };
202 const auto dof_max_order = get_element_max_dofs_order();
203 return (dof_max_order > 0) ? (dof_max_order - 1) / 2 : 0;
204};
205
206template <typename E>
209
210 auto &calc_mesh = this->mField.get_moab();
211
212 auto name = [&](auto tag) {
213 std::string name;
214 CHK_MOAB_THROW(calc_mesh.tag_get_name(tag, name), "get name");
215 return name;
216 };
217
218 auto data_type = [&](auto tag) {
219 moab::DataType data_type;
220 CHK_MOAB_THROW(calc_mesh.tag_get_data_type(tag, data_type),
221 "get data type");
222 return data_type;
223 };
224
225 auto type = [&](auto tag) {
226 moab::TagType type;
227 CHK_MOAB_THROW(calc_mesh.tag_get_type(tag, type), "get tag type");
228 return type;
229 };
230
231 auto length = [&](auto tag) {
232 int length;
233 CHK_MOAB_THROW(calc_mesh.tag_get_length(tag, length), "get length ");
234 return length;
235 };
236
237 auto default_value = [&](auto tag) {
238 const void *def_val;
239 int size;
240 CHK_MOAB_THROW(calc_mesh.tag_get_default_value(tag, def_val, size),
241 "get default tag value");
242 return def_val;
243 };
244
245 auto data_ptr = [&](auto tag) {
246 const void *tag_data;
247 auto core_ent = this->getFEEntityHandle();
248 CHK_MOAB_THROW(calc_mesh.tag_get_by_ptr(tag, &core_ent, 1, &tag_data),
249 "get tag data");
250 return tag_data;
251 };
252
253 for (auto tag : tagsToTransfer) {
254 Tag tag_postproc;
255 CHKERR getPostProcMesh().tag_get_handle(
256 name(tag).c_str(), length(tag), data_type(tag), tag_postproc,
257 type(tag) | MB_TAG_CREAT, default_value(tag));
258 CHKERR getPostProcMesh().tag_clear_data(tag_postproc, postProcElements,
259 data_ptr(tag));
260 }
261
263};
264
265template <typename E>
267 MoFEM::Interface &m_field)
268 : E(m_field) {}
269
270template <typename E>
272 MoFEM::Interface &m_field, boost::shared_ptr<moab::Core> core_mesh_ptr)
274 coreMeshPtr = core_mesh_ptr;
275}
276
277template <typename E>
279 ParallelComm *pcomm_post_proc_mesh =
280 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
281 if (pcomm_post_proc_mesh != NULL)
282 delete pcomm_post_proc_mesh;
283}
284
285template <typename E> int PostProcBrokenMeshInMoabBase<E>::getRule(int order) {
286 return -1;
287};
288
289template <typename E>
292
293 auto type = type_from_handle(this->getFEEntityHandle());
294
296
297 try {
298 ref_ele = refElementsMap.at(type);
299 } catch (const out_of_range &e) {
300 SETERRQ1(
301 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
302 "Generation of reference elements for type <%s> is not implemented",
303 moab::CN::EntityTypeName(type));
304 }
305
306 auto set_gauss_pts = [&](auto &level_gauss_pts_on_ref_mesh, auto &level_ref,
307 auto &level_shape_functions,
308
309 auto start_vert_handle, auto start_ele_handle,
310 auto &verts_array, auto &conn, auto &ver_count,
311 auto &ele_count
312
313 ) {
315
316 size_t level = getMaxLevel();
317 level = std::min(level, level_gauss_pts_on_ref_mesh.size() - 1);
318
319 auto &level_ref_gauss_pts = level_gauss_pts_on_ref_mesh[level];
320 auto &level_ref_ele = level_ref[level];
321 auto &shape_functions = level_shape_functions[level];
322 E::gaussPts.resize(level_ref_gauss_pts.size1(), level_ref_gauss_pts.size2(),
323 false);
324 noalias(E::gaussPts) = level_ref_gauss_pts;
325
326 const auto fe_ent = E::numeredEntFiniteElementPtr->getEnt();
327 auto get_fe_coords = [&]() {
328 const EntityHandle *conn;
329 int num_nodes;
331 E::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes, true),
332 "error get connectivity");
333 VectorDouble coords(num_nodes * 3);
335 E::mField.get_moab().get_coords(conn, num_nodes, &*coords.begin()),
336 "error get coordinates");
337 return coords;
338 };
339
340 auto coords = get_fe_coords();
341
342 const int num_nodes = level_ref_gauss_pts.size2();
343 mapGaussPts.resize(level_ref_gauss_pts.size2());
344
347 &*shape_functions.data().begin());
349 &verts_array[0][ver_count], &verts_array[1][ver_count],
350 &verts_array[2][ver_count]);
351 for (int gg = 0; gg != num_nodes; ++gg, ++ver_count) {
352
353 mapGaussPts[gg] = start_vert_handle + ver_count;
354
355 auto set_float_precision = [](const double x) {
356 if (std::abs(x) < std::numeric_limits<float>::epsilon())
357 return 0.;
358 else
359 return x;
360 };
361
362 t_coords(i) = 0;
363 auto t_ele_coords = getFTensor1FromArray<3, 3>(coords);
364 for (int nn = 0; nn != CN::VerticesPerEntity(type); ++nn) {
365 t_coords(i) += t_n * t_ele_coords(i);
366 ++t_ele_coords;
367 ++t_n;
368 }
369
370 for (auto ii : {0, 1, 2})
371 t_coords(ii) = set_float_precision(t_coords(ii));
372
373 ++t_coords;
374 }
375
376 Tag th;
377 int def_in_the_loop = -1;
378 CHKERR getPostProcMesh().tag_get_handle(
379 "NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER, th, MB_TAG_CREAT | MB_TAG_SPARSE,
380 &def_in_the_loop);
381
382 postProcElements.clear();
383 const int num_el = level_ref_ele.size1();
384 const int num_nodes_on_ele = level_ref_ele.size2();
385 auto start_e = start_ele_handle + ele_count;
386 postProcElements = Range(start_e, start_e + num_el - 1);
387 for (auto tt = 0; tt != level_ref_ele.size1(); ++tt, ++ele_count) {
388 for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
389 conn[num_nodes_on_ele * ele_count + nn] =
390 mapGaussPts[level_ref_ele(tt, nn)];
391 }
392 }
393
394 const int n_in_the_loop = E::nInTheLoop;
395 CHKERR getPostProcMesh().tag_clear_data(th, postProcElements,
396 &n_in_the_loop);
397 CHKERR transferTags();
398
400 };
401
402 CHKERR set_gauss_pts(
403
404 ref_ele->levelGaussPtsOnRefMesh, ref_ele->levelRef,
405 ref_ele->levelShapeFunctions,
406
407 ref_ele->startingVertEleHandle, ref_ele->startingEleHandle,
408 ref_ele->verticesOnEleArrays, ref_ele->eleConn, ref_ele->countVertEle,
409 ref_ele->countEle
410
411 );
412
414};
415
416template <typename E>
419
420 auto get_ref_ele = [&](const EntityType type) {
421 PostProcGenerateRefMeshPtr ref_ele_ptr;
422
423 try {
424
425 ref_ele_ptr = refElementsMap.at(type);
426
427 } catch (const out_of_range &e) {
428
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->generateReferenceElementMesh(),
453 "Error when generating reference element");
454 }
455
456 refElementsMap[type] = ref_ele_ptr;
457
458 return ref_ele_ptr;
459 };
460
461 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
462
463 auto miit =
464 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
465 boost::make_tuple(this->getFEName(), this->getLoFERank()));
466 auto hi_miit =
467 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
468 boost::make_tuple(this->getFEName(), this->getHiFERank()));
469
470 const int number_of_ents_in_the_loop = this->getLoopSize();
471 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
472 SETERRQ(E::mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
473 "Wrong size of indicies. Inconsistent size number of iterated "
474 "elements iterated by problem and from range.");
475 }
476
477 for (auto &m : refElementsMap) {
478 m.second->nbVertices = 0;
479 m.second->nbEles = 0;
480 m.second->countEle = 0;
481 m.second->countVertEle = 0;
482 }
483
484 for (; miit != hi_miit; ++miit) {
485 auto type = (*miit)->getEntType();
486 auto ref_ele = get_ref_ele(type);
487
488 // Set pointer to element. So that getDataVectorDofsPtr in getMaxLevel
489 // can work
490 E::numeredEntFiniteElementPtr = *miit;
491 bool add = true;
492 if (E::exeTestHook) {
493 add = E::exeTestHook(this);
494 }
495
496 if (add) {
497 size_t level = getMaxLevel();
498 level = std::min(level, ref_ele->levelGaussPtsOnRefMesh.size() - 1);
499 ref_ele->nbVertices += ref_ele->levelGaussPtsOnRefMesh[level].size2();
500 ref_ele->nbEles += ref_ele->levelRef[level].size1();
501 }
502 }
503
504 auto alloc_vertices_and_elements_on_post_proc_mesh = [&]() {
506
507 ReadUtilIface *iface;
508 CHKERR getPostProcMesh().query_interface(iface);
509
510 for (auto &m : refElementsMap) {
511 if (m.second->nbEles) {
512 CHKERR iface->get_node_coords(3, m.second->nbVertices, 0,
513 m.second->startingVertEleHandle,
514 m.second->verticesOnEleArrays);
515 CHKERR iface->get_element_connect(
516 m.second->nbEles, m.second->levelRef[0].size2(), m.first, 0,
517 m.second->startingEleHandle, m.second->eleConn);
518
519 m.second->countEle = 0;
520 m.second->countVertEle = 0;
521 }
522 }
523
525 };
526
527 CHKERR alloc_vertices_and_elements_on_post_proc_mesh();
528
530}
531
532template <typename E>
535
536 auto update_elements = [&]() {
537 ReadUtilIface *iface;
538 CHKERR getPostProcMesh().query_interface(iface);
540
541 Range ents;
542 for (auto &m : refElementsMap) {
543 if (m.second->nbEles) {
544 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
545 << "Update < " << moab::CN::EntityTypeName(m.first)
546 << " number of processed " << m.second->countEle;
547 CHKERR iface->update_adjacencies(
548 m.second->startingEleHandle, m.second->countEle,
549 m.second->levelRef[0].size2(), m.second->eleConn);
550 ents.merge(Range(m.second->startingEleHandle,
551 m.second->startingEleHandle + m.second->countEle - 1));
552 }
553 }
554
556 };
557
558 auto remove_obsolete_entities = [&]() {
560 Range ents, adj;
561 for (auto &m : refElementsMap) {
562 if (m.second->nbEles) {
563 ents.merge(Range(m.second->startingEleHandle,
564 m.second->startingEleHandle + m.second->countEle - 1));
565 const int dim = moab::CN::Dimension(m.first);
566 for (auto d = 1; d != dim; ++d) {
567 CHKERR getPostProcMesh().get_adjacencies(ents, d, false, adj,
568 moab::Interface::UNION);
569 }
570 }
571 }
572 CHKERR getPostProcMesh().delete_entities(adj);
574 };
575
576 auto set_proc_tags = [&]() {
578 ParallelComm *pcomm_post_proc_mesh =
579 ParallelComm::get_pcomm(&(getPostProcMesh()), MYPCOMM_INDEX);
580 if (pcomm_post_proc_mesh) {
581 Range ents;
582 for (auto &m : refElementsMap) {
583 if (m.second->nbEles) {
584 ents.merge(
585 Range(m.second->startingEleHandle,
586 m.second->startingEleHandle + m.second->countEle - 1));
587 }
588 }
589
590 int rank = E::mField.get_comm_rank();
591 CHKERR getPostProcMesh().tag_clear_data(pcomm_post_proc_mesh->part_tag(),
592 ents, &rank);
593 }
595 };
596
597 CHKERR update_elements();
598 CHKERR remove_obsolete_entities();
599 CHKERR set_proc_tags();
600
602}
603
604template <typename E>
607
608 ParallelComm *pcomm_post_proc_mesh =
609 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
610 if (pcomm_post_proc_mesh != NULL)
611 delete pcomm_post_proc_mesh;
612 CHKERR getPostProcMesh().delete_mesh();
613 pcomm_post_proc_mesh =
614 new ParallelComm(&(getPostProcMesh()), PETSC_COMM_WORLD);
615
616 CHKERR preProcPostProc();
617
619}
620
621template <typename E>
624
625 ParallelComm *pcomm_post_proc_mesh =
626 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
627 if(pcomm_post_proc_mesh == nullptr)
628 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "PComm not allocated");
629
630 CHKERR postProcPostProc();
631 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
632
634}
635
637 return mapGaussPts;
638}
639
641 if (!coreMeshPtr)
642 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Core mesh not set");
643 moab::Interface &post_proc_mesh_interface = *coreMeshPtr;
644 return post_proc_mesh_interface;
645}
646
647template <typename E>
649 return postProcElements;
650}
651
652template <typename E>
654PostProcBrokenMeshInMoabBase<E>::writeFile(const std::string file_name) {
656 ParallelComm *pcomm_post_proc_mesh =
657 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
658 if (pcomm_post_proc_mesh == NULL)
659 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
660 "ParallelComm not allocated");
661 CHKERR getPostProcMesh().write_file(file_name.c_str(), "MOAB",
662 "PARALLEL=WRITE_PART");
664};
665
666template <typename E> struct PostProcBrokenMeshInMoab;
667
668template <>
670 : public PostProcBrokenMeshInMoabBase<VolumeElementForcesAndSourcesCore> {
672 VolumeElementForcesAndSourcesCore>::PostProcBrokenMeshInMoabBase;
673};
674
675template <>
677 : public PostProcBrokenMeshInMoabBase<FaceElementForcesAndSourcesCore> {
679 FaceElementForcesAndSourcesCore>::PostProcBrokenMeshInMoabBase;
680};
681
682template <>
684 : public PostProcBrokenMeshInMoabBase<EdgeElementForcesAndSourcesCore> {
686 EdgeElementForcesAndSourcesCore>::PostProcBrokenMeshInMoabBase;
687};
688
689/**
690 * @brief Post post-proc data at points from hash maps
691 *
692 * @tparam DIM1 dimension of vector in data_map_vec and first column of
693 * data_map_may
694 * @tparam DIM2 dimension of second column in data_map_mat
695 */
696template <int DIM1, int DIM2>
698
699 using DataMapVec = std::map<std::string, boost::shared_ptr<VectorDouble>>;
700 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
701
702 /**
703 * @brief Construct a new OpPostProcMapInMoab object
704 *
705 * @param post_proc_mesh postprocessing mesh
706 * @param map_gauss_pts map of gauss points to nodes of postprocessing mesh
707 * @param data_map_scalar hash map of scalar values (string is name of the
708 * tag)
709 * @param data_map_vec hash map of vector values
710 * @param data_map_mat hash map of second order tensor values
711 * @param data_symm_map_mat hash map of symmetric second order tensor values
712 */
713 OpPostProcMapInMoab(moab::Interface &post_proc_mesh,
714 std::vector<EntityHandle> &map_gauss_pts,
715 DataMapVec data_map_scalar, DataMapMat data_map_vec,
716 DataMapMat data_map_mat, DataMapMat data_symm_map_mat)
719 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
720 dataMapScalar(data_map_scalar), dataMapVec(data_map_vec),
721 dataMapMat(data_map_mat), dataMapSymmMat(data_symm_map_mat) {
722 // Operator is only executed for vertices
723 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
724 }
725 MoFEMErrorCode doWork(int side, EntityType type,
727
728private:
729 moab::Interface &postProcMesh;
730 std::vector<EntityHandle> &mapGaussPts;
735};
736
737template <int DIM1, int DIM2>
742
743 std::array<double, 9> def;
744 std::fill(def.begin(), def.end(), 0);
745
746 auto get_tag = [&](const std::string name, size_t size) {
747 Tag th;
748 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
749 MB_TAG_CREAT | MB_TAG_SPARSE,
750 def.data());
751 return th;
752 };
753
754 MatrixDouble3by3 mat(3, 3);
755
756 auto set_vector_3d = [&](auto &t) -> MatrixDouble3by3 & {
757 mat.clear();
758 for (size_t r = 0; r != DIM1; ++r)
759 mat(0, r) = t(r);
760 return mat;
761 };
762
763 auto set_matrix_3d = [&](auto &t) -> MatrixDouble3by3 & {
764 mat.clear();
765 for (size_t r = 0; r != DIM1; ++r)
766 for (size_t c = 0; c != DIM2; ++c)
767 mat(r, c) = t(r, c);
768 return mat;
769 };
770
771 auto set_matrix_symm_3d = [&](auto &t) -> MatrixDouble3by3 & {
772 mat.clear();
773 for (size_t r = 0; r != DIM1; ++r)
774 for (size_t c = 0; c != DIM1; ++c)
775 mat(r, c) = t(r, c);
776 return mat;
777 };
778
779 auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
780 mat.clear();
781 mat(0, 0) = t;
782 return mat;
783 };
784
785 auto set_float_precision = [](const double x) {
786 if (std::abs(x) < std::numeric_limits<float>::epsilon())
787 return 0.;
788 else
789 return x;
790 };
791
792 auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
793 for (auto &v : mat.data())
794 v = set_float_precision(v);
795 return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
796 &*mat.data().begin());
797 };
798
799 for (auto &m : dataMapScalar) {
800 if (m.second) {
801 auto th = get_tag(m.first, 1);
802 auto t_scl = getFTensor0FromVec(*m.second);
803 auto nb_integration_pts = getGaussPts().size2();
804 size_t gg = 0;
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 size_t gg = 0;
818 for (int gg = 0; gg != nb_integration_pts; ++gg) {
819 CHKERR set_tag(th, gg, set_vector_3d(t_vec));
820 ++t_vec;
821 }
822 }
823 }
824
825 for (auto &m : dataMapMat) {
826 if (m.second) {
827 auto th = get_tag(m.first, 9);
828 auto t_mat = getFTensor2FromMat<DIM1, DIM2>(*m.second);
829 auto nb_integration_pts = getGaussPts().size2();
830 size_t gg = 0;
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 size_t gg = 0;
844 for (int gg = 0; gg != nb_integration_pts; ++gg) {
845 CHKERR set_tag(th, gg, set_matrix_symm_3d(t_mat));
846 ++t_mat;
847 }
848 }
849 }
850
852}
853
854template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseBeginImpl;
855
856template <typename E>
858 : protected PostProcBrokenMeshInMoabBase<E> {
859
861 boost::shared_ptr<moab::Core> core_mesh_ptr)
862 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr) {}
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 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr) {}
891
894
895protected:
897};
898
899template <typename PostProcEle> struct PostProcBrokenMeshInMoabBaseEndImpl;
900
901template <typename E>
903 : protected PostProcBrokenMeshInMoabBase<E> {
904
906 boost::shared_ptr<moab::Core> core_mesh_ptr)
907 : PostProcBrokenMeshInMoabBase<E>(m_field, core_mesh_ptr) {}
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:359
#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
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:1634
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)
PostProcBrokenMeshInMoabBaseContImpl(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr)
MoFEMErrorCode postProcess()
MoFEMErrorCode preProcess()
PostProcBrokenMeshInMoabBaseEndImpl(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr)
MoFEMErrorCode operator()()
FEMethod * getFEMethod()
MoFEMErrorCode preProcess()
MoFEMErrorCode postProcess()
auto & getPostProcElements()
Get postprocessing elements.
MoFEMErrorCode preProcess()
Generate vertices and elements.
PostProcBrokenMeshInMoabBase(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr)
virtual MoFEMErrorCode postProcPostProc()
MoFEMErrorCode setGaussPts(int order)
std::vector< Tag > tagsToTransfer
virtual MoFEMErrorCode transferTags()
PostProcBrokenMeshInMoabBase(MoFEM::Interface &m_field)
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
std::vector< MatrixDouble > levelShapeFunctions
PostProcGenerateRefMeshBase()
int countEle
int nbVertices
std::vector< MatrixDouble > levelGaussPtsOnRefMesh
int nbEles
MoFEMErrorCode getOptions()
std::string optPrefix
EntityHandle startingVertEleHandle
Element for postprocessing. Uses MoAB to generate post-processing mesh.