12#ifndef __POSTPROCBROKENMESHINMOABBASE_HPP
13#define __POSTPROCBROKENMESHINMOABBASE_HPP
25 std::vector<MatrixDouble>
28 std::vector<MatrixDouble>
55 boost::shared_ptr<PostProcGenerateRefMeshBase>;
162 boost::shared_ptr<moab::Core>
coreMeshPtr = boost::make_shared<moab::Core>();
164 boost::shared_ptr<moab::Core> core_mesh_ptr);
170 std::map<EntityType, PostProcGenerateRefMeshPtr>
186 std::vector<Tag> tags_to_transfer) {
188 tagsToTransfer.swap(tags_to_transfer);
193 auto get_element_max_dofs_order = [&]() {
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;
202 const auto dof_max_order = get_element_max_dofs_order();
203 return (dof_max_order > 0) ? (dof_max_order - 1) / 2 : 0;
210 auto &calc_mesh = this->mField.get_moab();
212 auto name = [&](
auto tag) {
218 auto data_type = [&](
auto tag) {
219 moab::DataType data_type;
225 auto type = [&](
auto tag) {
227 CHK_MOAB_THROW(calc_mesh.tag_get_type(tag, type),
"get tag type");
231 auto length = [&](
auto tag) {
233 CHK_MOAB_THROW(calc_mesh.tag_get_length(tag, length),
"get length ");
237 auto default_value = [&](
auto tag) {
240 CHK_MOAB_THROW(calc_mesh.tag_get_default_value(tag, def_val, size),
241 "get default tag value");
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),
253 for (
auto tag : tagsToTransfer) {
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,
279 ParallelComm *pcomm_post_proc_mesh =
281 if (pcomm_post_proc_mesh != NULL)
282 delete pcomm_post_proc_mesh;
298 ref_ele = refElementsMap.at(type);
299 }
catch (
const out_of_range &e) {
302 "Generation of reference elements for type <%s> is not implemented",
303 moab::CN::EntityTypeName(type));
306 auto set_gauss_pts = [&](
auto &level_gauss_pts_on_ref_mesh,
auto &level_ref,
307 auto &level_shape_functions,
309 auto start_vert_handle,
auto start_ele_handle,
310 auto &verts_array,
auto &conn,
auto &ver_count,
316 size_t level = getMaxLevel();
317 level = std::min(level, level_gauss_pts_on_ref_mesh.size() - 1);
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(),
324 noalias(E::gaussPts) = level_ref_gauss_pts;
326 const auto fe_ent = E::numeredEntFiniteElementPtr->getEnt();
327 auto get_fe_coords = [&]() {
331 E::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes,
true),
332 "error get connectivity");
335 E::mField.get_moab().get_coords(conn, num_nodes, &*coords.begin()),
336 "error get coordinates");
340 auto coords = get_fe_coords();
342 const int num_nodes = level_ref_gauss_pts.size2();
343 mapGaussPts.resize(level_ref_gauss_pts.size2());
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) {
353 mapGaussPts[gg] = start_vert_handle + ver_count;
355 auto set_float_precision = [](
const double x) {
356 if (std::abs(x) < std::numeric_limits<float>::epsilon())
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);
370 for (
auto ii : {0, 1, 2})
371 t_coords(ii) = set_float_precision(t_coords(ii));
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,
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)];
394 const int n_in_the_loop = E::nInTheLoop;
395 CHKERR getPostProcMesh().tag_clear_data(
th, postProcElements,
404 ref_ele->levelGaussPtsOnRefMesh, ref_ele->levelRef,
405 ref_ele->levelShapeFunctions,
407 ref_ele->startingVertEleHandle, ref_ele->startingEleHandle,
408 ref_ele->verticesOnEleArrays, ref_ele->eleConn, ref_ele->countVertEle,
420 auto get_ref_ele = [&](
const EntityType type) {
425 ref_ele_ptr = refElementsMap.at(type);
427 }
catch (
const out_of_range &e) {
431 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTET>>();
434 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBHEX>>();
437 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
440 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBQUAD>>();
443 ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBEDGE>>();
447 <<
"Generation of reference elements for type < "
448 << moab::CN::EntityTypeName(type) <<
" > is not implemented";
453 "Error when generating reference element");
456 refElementsMap[type] = ref_ele_ptr;
461 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
464 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
465 boost::make_tuple(this->getFEName(), this->getLoFERank()));
467 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
468 boost::make_tuple(this->getFEName(), this->getHiFERank()));
470 const int number_of_ents_in_the_loop = this->getLoopSize();
471 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
473 "Wrong size of indicies. Inconsistent size number of iterated "
474 "elements iterated by problem and from range.");
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;
484 for (; miit != hi_miit; ++miit) {
485 auto type = (*miit)->getEntType();
486 auto ref_ele = get_ref_ele(type);
490 E::numeredEntFiniteElementPtr = *miit;
492 if (E::exeTestHook) {
493 add = E::exeTestHook(
this);
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();
504 auto alloc_vertices_and_elements_on_post_proc_mesh = [&]() {
507 ReadUtilIface *iface;
508 CHKERR getPostProcMesh().query_interface(iface);
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);
519 m.second->countEle = 0;
520 m.second->countVertEle = 0;
527 CHKERR alloc_vertices_and_elements_on_post_proc_mesh();
536 auto update_elements = [&]() {
537 ReadUtilIface *iface;
538 CHKERR getPostProcMesh().query_interface(iface);
542 for (
auto &
m : refElementsMap) {
543 if (
m.second->nbEles) {
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));
558 auto remove_obsolete_entities = [&]() {
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);
572 CHKERR getPostProcMesh().delete_entities(adj);
576 auto set_proc_tags = [&]() {
578 ParallelComm *pcomm_post_proc_mesh =
579 ParallelComm::get_pcomm(&(getPostProcMesh()),
MYPCOMM_INDEX);
580 if (pcomm_post_proc_mesh) {
582 for (
auto &
m : refElementsMap) {
583 if (
m.second->nbEles) {
585 Range(
m.second->startingEleHandle,
586 m.second->startingEleHandle +
m.second->countEle - 1));
590 int rank = E::mField.get_comm_rank();
591 CHKERR getPostProcMesh().tag_clear_data(pcomm_post_proc_mesh->part_tag(),
598 CHKERR remove_obsolete_entities();
608 ParallelComm *pcomm_post_proc_mesh =
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);
625 ParallelComm *pcomm_post_proc_mesh =
627 if(pcomm_post_proc_mesh ==
nullptr)
630 CHKERR postProcPostProc();
631 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
643 moab::Interface &post_proc_mesh_interface = *coreMeshPtr;
644 return post_proc_mesh_interface;
649 return postProcElements;
656 ParallelComm *pcomm_post_proc_mesh =
658 if (pcomm_post_proc_mesh == NULL)
660 "ParallelComm not allocated");
661 CHKERR getPostProcMesh().write_file(file_name.c_str(),
"MOAB",
662 "PARALLEL=WRITE_PART");
696template <
int DIM1,
int DIM2>
699 using DataMapVec = std::map<std::string, boost::shared_ptr<VectorDouble>>;
700 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
714 std::vector<EntityHandle> &map_gauss_pts,
737template <
int DIM1,
int DIM2>
743 std::array<double, 9> def;
744 std::fill(def.begin(), def.end(), 0);
746 auto get_tag = [&](
const std::string name,
size_t size) {
748 CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE,
th,
749 MB_TAG_CREAT | MB_TAG_SPARSE,
758 for (
size_t r = 0; r != DIM1; ++r)
765 for (
size_t r = 0; r != DIM1; ++r)
766 for (
size_t c = 0;
c != DIM2; ++
c)
773 for (
size_t r = 0; r != DIM1; ++r)
774 for (
size_t c = 0;
c != DIM1; ++
c)
785 auto set_float_precision = [](
const double x) {
786 if (std::abs(x) < std::numeric_limits<float>::epsilon())
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());
799 for (
auto &
m : dataMapScalar) {
801 auto th = get_tag(
m.first, 1);
803 auto nb_integration_pts = getGaussPts().size2();
805 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
806 CHKERR set_tag(
th, gg, set_scalar(t_scl));
812 for (
auto &
m : dataMapVec) {
814 auto th = get_tag(
m.first, 3);
815 auto t_vec = getFTensor1FromMat<DIM1>(*
m.second);
816 auto nb_integration_pts = getGaussPts().size2();
818 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
819 CHKERR set_tag(
th, gg, set_vector_3d(t_vec));
825 for (
auto &
m : dataMapMat) {
827 auto th = get_tag(
m.first, 9);
828 auto t_mat = getFTensor2FromMat<DIM1, DIM2>(*
m.second);
829 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();
844 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
845 CHKERR set_tag(
th, gg, set_matrix_symm_3d(t_mat));
861 boost::shared_ptr<moab::Core> core_mesh_ptr)
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;
899template <
typename PostProcEle>
struct PostProcBrokenMeshInMoabBaseEndImpl;
906 boost::shared_ptr<moab::Core> core_mesh_ptr)
913 ParallelComm *pcomm_post_proc_mesh =
915 if (pcomm_post_proc_mesh ==
nullptr)
917 "PComm not allocated");
918 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)
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()()
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.
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
std::vector< MatrixDouble > levelShapeFunctions
PostProcGenerateRefMeshBase()
std::vector< MatrixDouble > levelGaussPtsOnRefMesh
MoFEMErrorCode getOptions()
EntityHandle startingVertEleHandle
Element for postprocessing. Uses MoAB to generate post-processing mesh.
Volume finite element base.