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 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,
279 :
E(m_field), optionsPrefix(opts_prefix) {}
284 std::string opts_prefix)
291 ParallelComm *pcomm_post_proc_mesh =
293 if (pcomm_post_proc_mesh != NULL)
294 delete pcomm_post_proc_mesh;
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())
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,
432auto 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();
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();
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);
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);
654 moab::Interface &post_proc_mesh_interface = *coreMeshPtr;
655 return post_proc_mesh_interface;
660 return postProcElements;
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,
734 std::fill(&O::doEntities[MBEDGE], &O::doEntities[MBMAXTYPE],
false);
748template <
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);
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);
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);
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));
869 std::string opts_prefix =
"")
874 ParallelComm *pcomm_post_proc_mesh =
876 if (pcomm_post_proc_mesh != NULL)
877 delete pcomm_post_proc_mesh;
879 pcomm_post_proc_mesh =
890template <
typename PostProcEle>
struct PostProcBrokenMeshInMoabBaseContImpl;
898 std::string opts_prefix =
"")
908template <
typename PostProcEle>
struct PostProcBrokenMeshInMoabBaseEndImpl;
916 std::string opts_prefix =
"")
923 ParallelComm *pcomm_post_proc_mesh =
925 if (pcomm_post_proc_mesh ==
nullptr)
927 "PComm not allocated");
928 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
#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 ...
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
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
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.
DataMapMat dataMapSymmMat
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
moab::Interface & postProcMesh
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="")
virtual ~PostProcBrokenMeshInMoabBase()
PostProcBrokenMeshInMoabBase(MoFEM::Interface &m_field, std::string opts_prefix="")
std::vector< Tag > tagsToTransfer
std::string optionsPrefix
Prefix for options.
virtual MoFEMErrorCode transferTags()
virtual MoFEMErrorCode preProcPostProc()
virtual MoFEMErrorCode postProcPostProc()
MoFEMErrorCode setGaussPts(int order)
auto & getMapGaussPts()
Get vector of vectors associated to integration points.
boost::shared_ptr< moab::Core > coreMeshPtr
std::vector< EntityHandle > mapGaussPts
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="")
MoFEMErrorCode postProcess()
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.
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.