12 #ifndef __POSTPROCBROKENMESHINMOABBASE_HPP
13 #define __POSTPROCBROKENMESHINMOABBASE_HPP
25 std::vector<MatrixDouble>
28 std::vector<MatrixDouble>
54 boost::shared_ptr<PostProcGenerateRefMeshBase>;
97 std::string opts_prefix =
"");
99 boost::shared_ptr<moab::Core> core_mesh_ptr,
100 std::string opts_prefix =
"");
165 boost::shared_ptr<moab::Core>
coreMeshPtr = boost::make_shared<moab::Core>();
170 std::map<EntityType, PostProcGenerateRefMeshPtr>
185 template <
typename E>
187 std::vector<Tag> tags_to_transfer) {
189 tagsToTransfer.swap(tags_to_transfer);
194 auto get_element_max_dofs_order = [&]() {
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;
203 const auto dof_max_order = get_element_max_dofs_order();
204 return (dof_max_order > 0) ? (dof_max_order - 1) / 2 : 0;
207 template <
typename E>
211 auto &calc_mesh = this->mField.get_moab();
213 auto name = [&](
auto tag) {
219 auto data_type = [&](
auto tag) {
220 moab::DataType data_type;
226 auto type = [&](
auto tag) {
232 auto length = [&](
auto tag) {
234 CHK_MOAB_THROW(calc_mesh.tag_get_length(tag, length),
"get length ");
238 auto default_value = [&](
auto tag) {
241 CHK_MOAB_THROW(calc_mesh.tag_get_default_value(tag, def_val, size),
242 "get default tag value");
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),
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());
264 for (
auto tag : tagsToTransfer) {
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,
276 template <
typename E>
279 :
E(m_field), optionsPrefix(opts_prefix) {}
281 template <
typename E>
284 std::string opts_prefix)
286 coreMeshPtr = core_mesh_ptr;
289 template <
typename E>
291 ParallelComm *pcomm_post_proc_mesh =
293 if (pcomm_post_proc_mesh != NULL)
294 delete pcomm_post_proc_mesh;
301 template <
typename E>
310 ref_ele = refElementsMap.at(
type);
311 }
catch (
const out_of_range &e) {
314 "Generation of reference elements for type <%s> is not implemented",
315 moab::CN::EntityTypeName(
type));
318 auto set_gauss_pts = [&](
auto &level_gauss_pts_on_ref_mesh,
auto &level_ref,
319 auto &level_shape_functions,
321 auto start_vert_handle,
auto start_ele_handle,
322 auto &verts_array,
auto &conn,
auto &ver_count,
328 size_t level = getMaxLevel();
329 level = std::min(level, level_gauss_pts_on_ref_mesh.size() - 1);
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(),
336 noalias(E::gaussPts) = level_ref_gauss_pts;
338 const auto fe_ent = E::numeredEntFiniteElementPtr->getEnt();
339 auto get_fe_coords = [&]() {
343 E::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes,
true),
344 "error get connectivity");
347 E::mField.get_moab().get_coords(conn, num_nodes, &*coords.begin()),
348 "error get coordinates");
352 auto coords = get_fe_coords();
354 const int num_nodes = level_ref_gauss_pts.size2();
355 mapGaussPts.resize(level_ref_gauss_pts.size2());
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) {
365 mapGaussPts[gg] = start_vert_handle + ver_count;
367 auto set_float_precision = [](
const double x) {
368 if (std::abs(x) < std::numeric_limits<float>::epsilon())
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);
382 for (
auto ii : {0, 1, 2})
383 t_coords(ii) = set_float_precision(t_coords(ii));
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,
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)];
406 const int n_in_the_loop = E::nInTheLoop;
407 CHKERR getPostProcMesh().tag_clear_data(
th, postProcElements,
416 ref_ele->levelGaussPtsOnRefMesh, ref_ele->levelRef,
417 ref_ele->levelShapeFunctions,
419 ref_ele->startingVertEleHandle, ref_ele->startingEleHandle,
420 ref_ele->verticesOnEleArrays, ref_ele->eleConn, ref_ele->countVertEle,
428 template <
typename E>
432 auto get_ref_ele = [&](
const EntityType
type) {
435 auto it = refElementsMap.find(
type);
436 if (it != refElementsMap.end()) {
437 ref_ele_ptr = it->second;
441 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTET>>();
444 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBHEX>>();
447 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
450 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBQUAD>>();
453 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBEDGE>>();
457 <<
"Generation of reference elements for type < "
458 << moab::CN::EntityTypeName(
type) <<
" > is not implemented";
464 "Error when generating reference element");
466 refElementsMap[
type] = ref_ele_ptr;
472 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
475 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
476 boost::make_tuple(this->getFEName(), this->getLoFERank()));
478 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
479 boost::make_tuple(this->getFEName(), this->getHiFERank()));
481 const int number_of_ents_in_the_loop = this->getLoopSize();
482 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
484 "Wrong size of indicies. Inconsistent size number of iterated "
485 "elements iterated by problem and from range.");
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;
495 for (; miit != hi_miit; ++miit) {
496 auto type = (*miit)->getEntType();
497 auto ref_ele = get_ref_ele(
type);
501 E::numeredEntFiniteElementPtr = *miit;
503 if (E::exeTestHook) {
504 add = E::exeTestHook(
this);
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();
515 auto alloc_vertices_and_elements_on_post_proc_mesh = [&]() {
518 ReadUtilIface *iface;
519 CHKERR getPostProcMesh().query_interface(iface);
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);
530 m.second->countEle = 0;
531 m.second->countVertEle = 0;
538 CHKERR alloc_vertices_and_elements_on_post_proc_mesh();
543 template <
typename E>
547 auto update_elements = [&]() {
548 ReadUtilIface *iface;
549 CHKERR getPostProcMesh().query_interface(iface);
553 for (
auto &
m : refElementsMap) {
554 if (
m.second->nbEles) {
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));
569 auto remove_obsolete_entities = [&]() {
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);
583 CHKERR getPostProcMesh().delete_entities(adj);
587 auto set_proc_tags = [&]() {
589 ParallelComm *pcomm_post_proc_mesh =
590 ParallelComm::get_pcomm(&(getPostProcMesh()),
MYPCOMM_INDEX);
591 if (pcomm_post_proc_mesh) {
593 for (
auto &
m : refElementsMap) {
594 if (
m.second->nbEles) {
596 Range(
m.second->startingEleHandle,
597 m.second->startingEleHandle +
m.second->countEle - 1));
601 int rank = E::mField.get_comm_rank();
602 CHKERR getPostProcMesh().tag_clear_data(pcomm_post_proc_mesh->part_tag(),
609 CHKERR remove_obsolete_entities();
615 template <
typename E>
619 ParallelComm *pcomm_post_proc_mesh =
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);
632 template <
typename E>
636 ParallelComm *pcomm_post_proc_mesh =
638 if(pcomm_post_proc_mesh ==
nullptr)
641 CHKERR postProcPostProc();
642 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
655 return post_proc_mesh_interface;
658 template <
typename E>
660 return postProcElements;
663 template <
typename E>
667 ParallelComm *pcomm_post_proc_mesh =
669 if (pcomm_post_proc_mesh == NULL)
671 "ParallelComm not allocated");
672 CHKERR getPostProcMesh().write_file(file_name.c_str(),
"MOAB",
673 "PARALLEL=WRITE_PART");
711 using DataMapVec = std::map<std::string, boost::shared_ptr<VectorDouble>>;
712 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
726 std::vector<EntityHandle> &map_gauss_pts,
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) {
734 std::fill(&O::doEntities[MBEDGE], &O::doEntities[MBMAXTYPE],
false);
748 template <
int DIM1,
int DIM2,
typename O>
754 std::array<double, 9> def;
755 std::fill(def.begin(), def.end(), 0);
757 auto get_tag = [&](
const std::string name,
size_t size) {
759 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE,
th,
760 MB_TAG_CREAT | MB_TAG_SPARSE,
769 for (
size_t r = 0;
r !=
DIM1; ++
r)
776 for (
size_t r = 0;
r !=
DIM1; ++
r)
777 for (
size_t c = 0;
c !=
DIM2; ++
c)
784 for (
size_t r = 0;
r !=
DIM1; ++
r)
785 for (
size_t c = 0;
c !=
DIM1; ++
c)
796 auto set_float_precision = [](
const double x) {
797 if (std::abs(x) < std::numeric_limits<float>::epsilon())
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());
810 for (
auto &
m : dataMapScalar) {
812 auto th = get_tag(
m.first, 1);
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));
822 for (
auto &
m : dataMapVec) {
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));
834 for (
auto &
m : dataMapMat) {
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));
846 for (
auto &
m : dataMapSymmMat) {
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));
863 template <
typename E>
869 std::string opts_prefix =
"")
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);
892 template <
typename E>
898 std::string opts_prefix =
"")
910 template <
typename E>
916 std::string opts_prefix =
"")
923 ParallelComm *pcomm_post_proc_mesh =
924 ParallelComm::get_pcomm(&this->getPostProcMesh(),
MYPCOMM_INDEX);
925 if (pcomm_post_proc_mesh ==
nullptr)
927 "PComm not allocated");
928 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
953 template <
typename E>
967 #endif //__POSTPROCBROKENMESHINMOABBASE_HPP