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>
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;
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) {
228 CHK_MOAB_THROW(calc_mesh.tag_get_type(tag, type),
"get tag type");
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 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),
254 for (
auto tag : tagsToTransfer) {
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,
269 :
E(m_field), optionsPrefix(opts_prefix) {}
274 std::string opts_prefix)
281 ParallelComm *pcomm_post_proc_mesh =
283 if (pcomm_post_proc_mesh != NULL)
284 delete pcomm_post_proc_mesh;
300 ref_ele = refElementsMap.at(type);
301 }
catch (
const out_of_range &e) {
304 "Generation of reference elements for type <%s> is not implemented",
305 moab::CN::EntityTypeName(type));
308 auto set_gauss_pts = [&](
auto &level_gauss_pts_on_ref_mesh,
auto &level_ref,
309 auto &level_shape_functions,
311 auto start_vert_handle,
auto start_ele_handle,
312 auto &verts_array,
auto &conn,
auto &ver_count,
318 size_t level = getMaxLevel();
319 level = std::min(level, level_gauss_pts_on_ref_mesh.size() - 1);
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(),
326 noalias(E::gaussPts) = level_ref_gauss_pts;
328 const auto fe_ent = E::numeredEntFiniteElementPtr->getEnt();
329 auto get_fe_coords = [&]() {
333 E::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes,
true),
334 "error get connectivity");
337 E::mField.get_moab().get_coords(conn, num_nodes, &*coords.begin()),
338 "error get coordinates");
342 auto coords = get_fe_coords();
344 const int num_nodes = level_ref_gauss_pts.size2();
345 mapGaussPts.resize(level_ref_gauss_pts.size2());
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) {
355 mapGaussPts[gg] = start_vert_handle + ver_count;
357 auto set_float_precision = [](
const double x) {
358 if (std::abs(x) < std::numeric_limits<float>::epsilon())
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);
372 for (
auto ii : {0, 1, 2})
373 t_coords(ii) = set_float_precision(t_coords(ii));
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,
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)];
396 const int n_in_the_loop = E::nInTheLoop;
397 CHKERR getPostProcMesh().tag_clear_data(
th, postProcElements,
406 ref_ele->levelGaussPtsOnRefMesh, ref_ele->levelRef,
407 ref_ele->levelShapeFunctions,
409 ref_ele->startingVertEleHandle, ref_ele->startingEleHandle,
410 ref_ele->verticesOnEleArrays, ref_ele->eleConn, ref_ele->countVertEle,
422 auto get_ref_ele = [&](
const EntityType type) {
427 ref_ele_ptr = refElementsMap.at(type);
429 }
catch (
const out_of_range &e) {
433 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTET>>();
436 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBHEX>>();
439 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
442 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBQUAD>>();
445 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBEDGE>>();
449 <<
"Generation of reference elements for type < "
450 << moab::CN::EntityTypeName(type) <<
" > is not implemented";
456 "Error when generating reference element");
459 refElementsMap[type] = ref_ele_ptr;
464 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
467 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
468 boost::make_tuple(this->getFEName(), this->getLoFERank()));
470 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
471 boost::make_tuple(this->getFEName(), this->getHiFERank()));
473 const int number_of_ents_in_the_loop = this->getLoopSize();
474 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
476 "Wrong size of indicies. Inconsistent size number of iterated "
477 "elements iterated by problem and from range.");
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;
487 for (; miit != hi_miit; ++miit) {
488 auto type = (*miit)->getEntType();
489 auto ref_ele = get_ref_ele(type);
493 E::numeredEntFiniteElementPtr = *miit;
495 if (E::exeTestHook) {
496 add = E::exeTestHook(
this);
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();
507 auto alloc_vertices_and_elements_on_post_proc_mesh = [&]() {
510 ReadUtilIface *iface;
511 CHKERR getPostProcMesh().query_interface(iface);
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);
522 m.second->countEle = 0;
523 m.second->countVertEle = 0;
530 CHKERR alloc_vertices_and_elements_on_post_proc_mesh();
539 auto update_elements = [&]() {
540 ReadUtilIface *iface;
541 CHKERR getPostProcMesh().query_interface(iface);
545 for (
auto &
m : refElementsMap) {
546 if (
m.second->nbEles) {
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));
561 auto remove_obsolete_entities = [&]() {
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);
575 CHKERR getPostProcMesh().delete_entities(adj);
579 auto set_proc_tags = [&]() {
581 ParallelComm *pcomm_post_proc_mesh =
582 ParallelComm::get_pcomm(&(getPostProcMesh()),
MYPCOMM_INDEX);
583 if (pcomm_post_proc_mesh) {
585 for (
auto &
m : refElementsMap) {
586 if (
m.second->nbEles) {
588 Range(
m.second->startingEleHandle,
589 m.second->startingEleHandle +
m.second->countEle - 1));
593 int rank = E::mField.get_comm_rank();
594 CHKERR getPostProcMesh().tag_clear_data(pcomm_post_proc_mesh->part_tag(),
601 CHKERR remove_obsolete_entities();
611 ParallelComm *pcomm_post_proc_mesh =
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);
628 ParallelComm *pcomm_post_proc_mesh =
630 if(pcomm_post_proc_mesh ==
nullptr)
633 CHKERR postProcPostProc();
634 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
646 moab::Interface &post_proc_mesh_interface = *coreMeshPtr;
647 return post_proc_mesh_interface;
652 return postProcElements;
659 ParallelComm *pcomm_post_proc_mesh =
661 if (pcomm_post_proc_mesh == NULL)
663 "ParallelComm not allocated");
664 CHKERR getPostProcMesh().write_file(file_name.c_str(),
"MOAB",
665 "PARALLEL=WRITE_PART");
699template <
int DIM1,
int DIM2>
702 using DataMapVec = std::map<std::string, boost::shared_ptr<VectorDouble>>;
703 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
717 std::vector<EntityHandle> &map_gauss_pts,
740template <
int DIM1,
int DIM2>
746 std::array<double, 9> def;
747 std::fill(def.begin(), def.end(), 0);
749 auto get_tag = [&](
const std::string name,
size_t size) {
751 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE,
th,
752 MB_TAG_CREAT | MB_TAG_SPARSE,
761 for (
size_t r = 0; r != DIM1; ++r)
768 for (
size_t r = 0; r != DIM1; ++r)
769 for (
size_t c = 0;
c != DIM2; ++
c)
776 for (
size_t r = 0; r != DIM1; ++r)
777 for (
size_t c = 0;
c != DIM1; ++
c)
788 auto set_float_precision = [](
const double x) {
789 if (std::abs(x) < std::numeric_limits<float>::epsilon())
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());
802 for (
auto &
m : dataMapScalar) {
804 auto th = get_tag(
m.first, 1);
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));
814 for (
auto &
m : dataMapVec) {
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));
826 for (
auto &
m : dataMapMat) {
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));
838 for (
auto &
m : dataMapSymmMat) {
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));
861 std::string opts_prefix =
"")
866 ParallelComm *pcomm_post_proc_mesh =
868 if (pcomm_post_proc_mesh != NULL)
869 delete pcomm_post_proc_mesh;
871 pcomm_post_proc_mesh =
882template <
typename PostProcEle>
struct PostProcBrokenMeshInMoabBaseContImpl;
890 std::string opts_prefix =
"")
900template <
typename PostProcEle>
struct PostProcBrokenMeshInMoabBaseEndImpl;
908 std::string opts_prefix =
"")
915 ParallelComm *pcomm_post_proc_mesh =
917 if (pcomm_post_proc_mesh ==
nullptr)
919 "PComm not allocated");
920 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#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()
#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
#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 ...
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.
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.
implementation of Data Operators for Forces and Sources
auto type_from_handle(const EntityHandle h)
get type from entity handle
boost::shared_ptr< PostProcGenerateRefMeshBase > PostProcGenerateRefMeshPtr
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
constexpr double t
plate stiffness
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 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
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()
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()()
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.
virtual ~PostProcBrokenMeshInMoabBase()
MoFEMErrorCode setTagsToTransfer(std::vector< Tag > tags_to_transfer)
Set tags to be transferred to post-processing mesh.
std::map< EntityType, PostProcGenerateRefMeshPtr > refElementsMap
virtual MoFEMErrorCode preProcPostProc()
EntityHandle startingEleHandle
virtual MoFEMErrorCode generateReferenceElementMesh()=0
virtual ~PostProcGenerateRefMeshBase()=default
std::vector< ublas::matrix< int > > levelRef
std::vector< double * > verticesOnEleArrays
virtual MoFEMErrorCode getOptions(std::string prefix)
std::vector< MatrixDouble > levelShapeFunctions
PostProcGenerateRefMeshBase()
std::vector< MatrixDouble > levelGaussPtsOnRefMesh
EntityHandle startingVertEleHandle
Element for postprocessing. Uses MoAB to generate post-processing mesh.
Volume finite element base.