v0.15.0
Loading...
Searching...
No Matches
MoFEM::MeshRefinement Struct Reference

Mesh refinement interface. More...

#include "src/interfaces/MeshRefinement.hpp"

Inheritance diagram for MoFEM::MeshRefinement:
[legend]
Collaboration diagram for MoFEM::MeshRefinement:
[legend]

Classes

struct  SetParent
 

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 
 MeshRefinement (const MoFEM::Core &core)
 
virtual ~MeshRefinement ()=default
 
MoFEMErrorCode addVerticesInTheMiddleOfEdges (const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
 make vertices in the middle of edges in meshset and add them to refinement levels defined by bit
 
MoFEMErrorCode addVerticesInTheMiddleOfEdges (const Range &edges, const BitRefLevel &bit, int verb=QUIET, EntityHandle start_v=0)
 make vertices in the middle of edges in meshset and add them to Refinement levels defined by bit
 
MoFEMErrorCode refineTets (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TET in the meshset
 
MoFEMErrorCode refineTets (const Range &tets, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TET in the meshset
 
MoFEMErrorCode refineTetsHangingNodes (const Range &tets, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TET in the meshset
 
MoFEMErrorCode refineTetsHangingNodes (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TET in the meshset
 
MoFEMErrorCode refinePrisms (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET)
 refine PRISM in the meshset
 
MoFEMErrorCode refineMeshset (const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET)
 refine meshset, i.e. add child of refined entities to meshset
 
MoFEMErrorCode refineTris (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine triangles in the meshset
 
MoFEMErrorCode refineTris (const Range &tris, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TRI in the meshset
 
MoFEMErrorCode refineTrisHangingNodes (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TRI in the meshset
 
MoFEMErrorCode refineTrisHangingNodes (const Range &tris, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TRI in the meshset
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register interface.
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface reference to pointer of interface.
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface.
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface.
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface.
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface.
 
virtual ~UnknownInterface ()=default
 

Public Attributes

MoFEM::CorecOre
 

Private Types

using SetEdgeBitsFun
 Functions setting edges for refinemnt on enetity level.
 

Private Member Functions

MoFEMErrorCode refineTets (const Range &tets, const BitRefLevel &bit, SetEdgeBitsFun set_edge_bits, int verb, const bool debug)
 refine TET in the meshset
 
MoFEMErrorCode refineTris (const Range &tris, const BitRefLevel &bit, SetEdgeBitsFun set_edge_bits, int verb, const bool debug)
 

Additional Inherited Members

- Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode getLibVersion (Version &version)
 Get library version.
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version.
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version.
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version.
 

Detailed Description

Mesh refinement interface.

Currently this class is abstraction to Core interface. In future should be outsourced as independent interface.

Bug

Not working on partitioned meshes

Need to be implemented as a stand alone interface not as a part of core structure which should be only basic database

If outsourced, class member functions should follow name convention

Spelling mistakes will be corrected with names fix to follow name convetion

Examples
EshelbianPlasticity.cpp, bernstein_bezier_generate_base.cpp, child_and_parent.cpp, hanging_node_approx.cpp, and split_sideset.cpp.

Definition at line 26 of file MeshRefinement.hpp.

Member Typedef Documentation

◆ SetEdgeBitsFun

Initial value:
boost::function<
MoFEMErrorCode(moab::Interface &moab,
&ref_parent_ents_view,
EntityHandle tet, BitRefEdges &parent_edges_bit,
EntityHandle *edge_new_nodes, int *split_edges
)>
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFEDGES_SIZE > BitRefEdges
Definition Types.hpp:34
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_non_unique< const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< boost::shared_ptr< RefEntity >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_ordered_parent_entity

Functions setting edges for refinemnt on enetity level.

Definition at line 195 of file MeshRefinement.hpp.

Constructor & Destructor Documentation

◆ MeshRefinement()

MoFEM::MeshRefinement::MeshRefinement ( const MoFEM::Core & core)

Definition at line 39 of file MeshRefinement.cpp.

40 : cOre(const_cast<Core &>(core)) {}

◆ ~MeshRefinement()

virtual MoFEM::MeshRefinement::~MeshRefinement ( )
virtualdefault

Member Function Documentation

◆ addVerticesInTheMiddleOfEdges() [1/2]

MoFEMErrorCode MoFEM::MeshRefinement::addVerticesInTheMiddleOfEdges ( const EntityHandle meshset,
const BitRefLevel & bit,
const bool recursive = false,
int verb = QUIET,
EntityHandle start_v = 0 )

make vertices in the middle of edges in meshset and add them to refinement levels defined by bit

Takes entities fromm meshsets and queried recursively (get entities from meshsets in meshsets, usually have to be used for CUBIT meshset). If meshset does not contain any edges, get entities in dimension 3 and get edge adjacencies.

Parameters
EntityHandlemeshset
BitRefLevelbitLevel
recursiveIf true, meshsets containing meshsets are queried recursively. Returns the contents of meshsets, but not the meshsets themselves if true.
Examples
bernstein_bezier_generate_base.cpp.

Definition at line 42 of file MeshRefinement.cpp.

44 {
45 Interface &m_field = cOre;
46 moab::Interface &moab = m_field.get_moab();
48 Range edges;
49 CHKERR moab.get_entities_by_type(meshset, MBEDGE, edges, recursive);
50 if (edges.empty()) {
51 Range tets;
52 CHKERR moab.get_entities_by_type(meshset, MBTET, tets, recursive);
53 CHKERR moab.get_adjacencies(tets, 1, true, edges, moab::Interface::UNION);
54 if (tets.empty()) {
55 Range prisms;
56 CHKERR moab.get_entities_by_type(meshset, MBPRISM, prisms, recursive);
57 for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
58 const EntityHandle *conn;
59 int num_nodes;
60 CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
61 assert(num_nodes == 6);
62 //
63 Range edge;
64 CHKERR moab.get_adjacencies(&conn[0], 2, 1, true, edge);
65 assert(edge.size() == 1);
66 edges.insert(edge[0]);
67 edge.clear();
68 CHKERR moab.get_adjacencies(&conn[1], 2, 1, true, edge);
69 assert(edge.size() == 1);
70 edges.insert(edge[0]);
71 EntityHandle conn_edge2[] = {conn[2], conn[0]};
72 edge.clear();
73 CHKERR moab.get_adjacencies(conn_edge2, 2, 1, true, edge);
74 assert(edge.size() == 1);
75 edges.insert(edge[0]);
76 //
77 edge.clear();
78 CHKERR moab.get_adjacencies(&conn[3], 2, 1, true, edge);
79 assert(edge.size() == 1);
80 edges.insert(edge[0]);
81 edge.clear();
82 CHKERR moab.get_adjacencies(&conn[4], 2, 1, true, edge);
83 assert(edge.size() == 1);
84 edges.insert(edge[0]);
85 EntityHandle conn_edge8[] = {conn[5], conn[3]};
86 edge.clear();
87 CHKERR moab.get_adjacencies(conn_edge8, 2, 1, true, edge);
88 assert(edge.size() == 1);
89 edges.insert(edge[0]);
90 }
91 }
92 }
93 CHKERR addVerticesInTheMiddleOfEdges(edges, bit, verb, start_v);
95}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto bit
set bit
DeprecatedCoreInterface Interface
moab::Interface & get_moab()
Definition Core.hpp:322
MoFEMErrorCode addVerticesInTheMiddleOfEdges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit

◆ addVerticesInTheMiddleOfEdges() [2/2]

MoFEMErrorCode MoFEM::MeshRefinement::addVerticesInTheMiddleOfEdges ( const Range & edges,
const BitRefLevel & bit,
int verb = QUIET,
EntityHandle start_v = 0 )

make vertices in the middle of edges in meshset and add them to Refinement levels defined by bit

Takes entities from meshsets and queried recursively (get entities from meshsets in meshsets, usually have to be used for CUBIT meshset). If meshset does not contain any edges, get entities in dimension 3 and get edge adjacencies.

Parameters
Rangeconsisting edges for refine
BitRefLevelbitLevel
recursiveIf true, meshsets containing meshsets are queried recursively. Returns the contents of meshsets, but not the meshsets themselves if true.

Definition at line 96 of file MeshRefinement.cpp.

97 {
98 Interface &m_field = cOre;
99 moab::Interface &moab = m_field.get_moab();
100 auto refined_ents_ptr = m_field.get_ref_ents();
102 auto miit =
103 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
104 .lower_bound(boost::make_tuple(MBVERTEX, MBEDGE));
105 auto hi_miit =
106 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
107 .upper_bound(boost::make_tuple(MBVERTEX, MBEDGE));
109 ref_parent_ents_view.insert(miit, hi_miit);
110
111 Range edges = ents.subset_by_type(MBEDGE);
112 if (verb >= VERBOSE) {
113 MOFEM_TAG_AND_LOG("WORLD", Sev::verbose, "MeshRefinemnt")
114 << "ref level " << bit << " nb. edges to refine " << edges.size();
115 }
116
117 std::array<std::vector<double>, 3> vert_coords;
118 for (auto &vc : vert_coords)
119 vc.reserve(edges.size());
120
121 std::vector<EntityHandle> parent_edge;
122 parent_edge.reserve(edges.size());
123
124 std::array<double, 9> coords;
126 &coords[0], &coords[1], &coords[2]};
128 &coords[3], &coords[4], &coords[5]};
130 &coords[6], &coords[7], &coords[8]};
131
132 FTensor::Index<'i', 3> i;
133
134 Range add_bit;
135 for (auto p_eit = edges.pair_begin(); p_eit != edges.pair_end(); ++p_eit) {
136 auto miit_view = ref_parent_ents_view.lower_bound(p_eit->first);
137
138 Range edge_having_parent_vertex;
139 if (miit_view != ref_parent_ents_view.end()) {
140 for (auto hi_miit_view = ref_parent_ents_view.upper_bound(p_eit->second);
141 miit_view != hi_miit_view; ++miit_view) {
142 edge_having_parent_vertex.insert(edge_having_parent_vertex.end(),
143 miit_view->get()->getParentEnt());
144 add_bit.insert(add_bit.end(), miit_view->get()->getEnt());
145 }
146 }
147
148 Range add_vertex_edges =
149 subtract(Range(p_eit->first, p_eit->second), edge_having_parent_vertex);
150
151 for (auto e : add_vertex_edges)
152 parent_edge.emplace_back(e);
153
154 for (auto e : add_vertex_edges) {
155
156 const EntityHandle *conn;
157 int num_nodes;
158 CHKERR moab.get_connectivity(e, conn, num_nodes, false);
159 if(num_nodes == 2) {
160 CHKERR moab.get_coords(conn, num_nodes, coords.data());
161 t_0(i) += t_1(i);
162 t_0(i) *= 0.5;
163 for (auto j : {0, 1, 2})
164 vert_coords[j].emplace_back(t_0(j));
165
166 } else if (num_nodes == 3) {
167 for (auto j : {0, 1, 2})
168 vert_coords[j].emplace_back(t_2(j));
169 } else {
170 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
171 "edge should have 2 or 3 nodes but have %d", num_nodes);
172 }
173 }
174 }
175
176 CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit, bit);
177
178 if (!vert_coords[0].empty()) {
179 ReadUtilIface *read_util;
180 CHKERR moab.query_interface(read_util);
181 int num_nodes = vert_coords[0].size();
182 vector<double *> arrays_coord;
183 CHKERR read_util->get_node_coords(3, num_nodes, 0, start_v, arrays_coord);
184 Range verts(start_v, start_v + num_nodes - 1);
185 for (auto dd : {0, 1, 2}) {
186 std::copy(vert_coords[dd].begin(), vert_coords[dd].end(),
187 arrays_coord[dd]);
188 }
189 CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), verts,
190 &*parent_edge.begin());
191 CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
192 verts, bit, verb);
193 }
195}
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
@ VERBOSE
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
Tag get_th_RefParentHandle() const
Definition Core.hpp:197

◆ query_interface()

MoFEMErrorCode MoFEM::MeshRefinement::query_interface ( boost::typeindex::type_index type_index,
UnknownInterface ** iface ) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 33 of file MeshRefinement.cpp.

34 {
35 *iface = const_cast<MeshRefinement *>(this);
36 return 0;
37}
MeshRefinement(const MoFEM::Core &core)

◆ refineMeshset()

MoFEMErrorCode MoFEM::MeshRefinement::refineMeshset ( const EntityHandle meshset,
const BitRefLevel & bit,
const bool recursive = false,
int verb = QUIET )

refine meshset, i.e. add child of refined entities to meshset

Parameters
EntityHandlemeshset where to save the child refined entities
BitRefLevelbitLevel
recursiveIf true, meshsets containing meshsets are queried recursively. Returns the contents of meshsets, but not the meshsets themselves if true.

Definition at line 905 of file MeshRefinement.cpp.

907 {
908 Interface &m_field = cOre;
909 auto refined_ents_ptr = m_field.get_ref_ents();
911 auto miit = refined_ents_ptr->find(meshset);
912 if (miit == refined_ents_ptr->end()) {
913 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
914 "this meshset is not in ref database");
915 }
916 CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
917 meshset, bit, meshset, MBMAXTYPE, recursive, verb);
918 *(const_cast<RefEntity *>(miit->get())->getBitRefLevelPtr()) |= bit;
920}
RefEntityTmp< 0 > RefEntity
MoFEMErrorCode get_ref_ents(const RefEntity_multiIndex **refined_entities_ptr) const
Get ref entities multi-index from database.
Definition Core.cpp:810

◆ refinePrisms()

MoFEMErrorCode MoFEM::MeshRefinement::refinePrisms ( const EntityHandle meshset,
const BitRefLevel & bit,
int verb = QUIET )

refine PRISM in the meshset

Parameters
EntityHandlemeshset
BitRefLevelbitLevel

Definition at line 722 of file MeshRefinement.cpp.

723 {
724
725 Interface &m_field = cOre;
726 moab::Interface &moab = m_field.get_moab();
727 auto refined_ents_ptr = m_field.get_ref_ents();
728 auto refined_finite_elements_ptr = m_field.get_ref_finite_elements();
729
730 // FIXME: refinement is based on entity handlers, should work on global ids of
731 // nodes, this will allow parallelise algorithm in the future
732
734
735 RefElement_multiIndex_parents_view ref_ele_parent_view;
736 ref_ele_parent_view.insert(
737 refined_finite_elements_ptr->get<Ent_mi_tag>().lower_bound(
739 refined_finite_elements_ptr->get<Ent_mi_tag>().upper_bound(
741 auto &ref_ele_by_parent_and_ref_edges =
742 ref_ele_parent_view
743 .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
744 // find all vertices which parent is edge
745 auto &ref_ents_by_comp =
746 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
748 ref_parent_ents_view.insert(
749 ref_ents_by_comp.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE)),
750 ref_ents_by_comp.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE)));
751 Range prisms;
752 CHKERR moab.get_entities_by_type(meshset, MBPRISM, prisms, false);
753 Range::iterator pit = prisms.begin();
754 for (; pit != prisms.end(); pit++) {
755 auto miit_prism = refined_ents_ptr->get<Ent_mi_tag>().find(*pit);
756 if (miit_prism == refined_ents_ptr->end()) {
757 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
758 "this prism is not in ref database");
759 }
760 if (verb >= NOISY) {
761 MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "MeshRefinment")
762 << "ref prism " << **miit_prism;
763 }
764 // prism connectivity
765 int num_nodes;
766 const EntityHandle *conn;
767 CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
768 assert(num_nodes == 6);
769 // edges connectivity
770 EntityHandle edges[6];
771 for (int ee = 0; ee < 3; ee++) {
772 CHKERR moab.side_element(*pit, 1, ee, edges[ee]);
773 }
774 for (int ee = 6; ee < 9; ee++) {
775 CHKERR moab.side_element(*pit, 1, ee, edges[ee - 3]);
776 }
777 // detect split edges
778 BitRefEdges split_edges(0);
779 EntityHandle edge_nodes[6];
780 std::fill(&edge_nodes[0], &edge_nodes[6], no_handle);
781 for (int ee = 0; ee < 6; ee++) {
782 RefEntity_multiIndex_view_by_hashed_parent_entity::iterator miit_view =
783 ref_parent_ents_view.find(edges[ee]);
784 if (miit_view != ref_parent_ents_view.end()) {
785 if (((*miit_view)->getBitRefLevel() & bit).any()) {
786 edge_nodes[ee] = (*miit_view)->getEnt();
787 split_edges.set(ee);
788 }
789 }
790 }
791 if (split_edges.count() == 0) {
792 *(const_cast<RefEntity *>(miit_prism->get())->getBitRefLevelPtr()) |= bit;
793 if (verb >= VERY_NOISY)
794 MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "MeshRefinment")
795 << "no refinement";
796 continue;
797 }
798 // check consistency
799 if (verb >= NOISY) {
800 std::ostringstream ss;
801 ss << "prism split edges " << split_edges << " count "
802 << split_edges.count() << std::endl;
803 PetscPrintf(m_field.get_comm(), "%s", ss.str().c_str());
804 }
805 // prism ref
806 EntityHandle new_prism_conn[4 * 6];
807 std::fill(&new_prism_conn[0], &new_prism_conn[4 * 6], no_handle);
808 int nb_new_prisms = 0;
809 switch (split_edges.count()) {
810 case 0:
811 break;
812 case 2:
813 CHKERR prism_type_1(conn, split_edges, edge_nodes, new_prism_conn);
814 nb_new_prisms = 2;
815 break;
816 case 4:
817 CHKERR prism_type_2(conn, split_edges, edge_nodes, new_prism_conn);
818 nb_new_prisms = 3;
819 break;
820 case 6:
821 CHKERR prism_type_3(conn, split_edges, edge_nodes, new_prism_conn);
822 nb_new_prisms = 4;
823 break;
824 default:
825 std::ostringstream ss;
826 ss << split_edges << " : [ " << conn[0] << " " << conn[1] << " "
827 << conn[2] << " " << conn[3] << " " << conn[4] << " " << conn[5]
828 << " ]";
829 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "%s",
830 ss.str().c_str());
831 }
832 // find that prism
833 std::bitset<4> ref_prism_bit(0);
834 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.lower_bound(
835 boost::make_tuple(*pit, split_edges.to_ulong()));
836 auto hi_it_by_ref_edges = ref_ele_by_parent_and_ref_edges.upper_bound(
837 boost::make_tuple(*pit, split_edges.to_ulong()));
838 auto it_by_ref_edges2 = it_by_ref_edges;
839 for (int pp = 0; it_by_ref_edges2 != hi_it_by_ref_edges;
840 it_by_ref_edges2++, pp++) {
841 // Add this tet to this ref
842 *(const_cast<RefElement *>(it_by_ref_edges2->get())
843 ->getBitRefLevelPtr()) |= bit;
844 ref_prism_bit.set(pp, 1);
845 if (verb > 2) {
846 std::ostringstream ss;
847 ss << "is refined " << *it_by_ref_edges2 << std::endl;
848 PetscPrintf(m_field.get_comm(), "%s", ss.str().c_str());
849 }
850 }
851 if (it_by_ref_edges != hi_it_by_ref_edges) {
852 if (ref_prism_bit.count() != (unsigned int)nb_new_prisms)
853 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
854 "data inconsistency");
855 } else {
856 EntityHandle ref_prisms[4];
857 // create prism
858 for (int pp = 0; pp < nb_new_prisms; pp++) {
859 if (verb > 3) {
860 std::ostringstream ss;
861 ss << "ref prism " << ref_prism_bit << std::endl;
862 PetscPrintf(m_field.get_comm(), "%s", ss.str().c_str());
863 }
864 if (!ref_prism_bit.test(pp)) {
865 CHKERR moab.create_element(MBPRISM, &new_prism_conn[6 * pp], 6,
866 ref_prisms[pp]);
867 CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(),
868 &ref_prisms[pp], 1, &*pit);
869 CHKERR moab.tag_set_data(cOre.get_th_RefBitLevel(), &ref_prisms[pp],
870 1, &bit);
871 CHKERR moab.tag_set_data(cOre.get_th_RefBitEdge(), &ref_prisms[pp], 1,
872 &split_edges);
873 auto p_ent =
874 const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
875 ->insert(boost::shared_ptr<RefEntity>(new RefEntity(
876 m_field.get_basic_entity_data_ptr(), ref_prisms[pp])));
877 std::pair<RefElement_multiIndex::iterator, bool> p_fe;
878 try {
879 p_fe =
880 const_cast<RefElement_multiIndex *>(refined_finite_elements_ptr)
881 ->insert(boost::shared_ptr<RefElement>(
882 new RefElement_PRISM(*p_ent.first)));
883 } catch (MoFEMException const &e) {
884 SETERRQ(PETSC_COMM_SELF, e.errorCode, "%s", e.errorMessage);
885 }
886 ref_prism_bit.set(pp);
887 CHKERR cOre.addPrismToDatabase(ref_prisms[pp]);
888 if (verb > 2) {
889 std::ostringstream ss;
890 ss << "add prism: " << **p_fe.first << std::endl;
891 if (verb > 7) {
892 for (int nn = 0; nn < 6; nn++) {
893 ss << new_prism_conn[nn] << " ";
894 }
895 ss << std::endl;
896 }
897 PetscPrintf(m_field.get_comm(), "%s", ss.str().c_str());
898 }
899 }
900 }
901 }
902 }
904}
@ NOISY
@ VERY_NOISY
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< hashed_non_unique< const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< boost::shared_ptr< RefEntity >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_hashed_parent_entity
multi-index view of RefEntity by parent entity
multi_index_container< boost::shared_ptr< RefElement >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getEnt > > > > RefElement_multiIndex
EntityHandle get_id_for_max_type()
EntityHandle get_id_for_min_type()
MoFEMErrorCode prism_type_2(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
MoFEMErrorCode prism_type_3(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
MoFEMErrorCode prism_type_1(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
multi_index_container< boost::shared_ptr< RefElement >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getEnt > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getParentEnt > >, ordered_non_unique< tag< Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag >, composite_key< RefElement, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getParentEnt >, const_mem_fun< RefElement, int, &RefElement::getBitRefEdgesUlong > > > > > RefElement_multiIndex_parents_view
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition Common.hpp:12
MoFEMErrorCode addPrismToDatabase(const EntityHandle prism, int verb=DEFAULT_VERBOSITY)
add prim element
Definition Core.cpp:544
Tag get_th_RefBitLevel() const
Definition Core.hpp:198
Tag get_th_RefBitEdge() const
Definition Core.hpp:199

◆ refineTets() [1/3]

MoFEMErrorCode MoFEM::MeshRefinement::refineTets ( const EntityHandle meshset,
const BitRefLevel & bit,
int verb = QUIET,
const bool debug = false )

refine TET in the meshset

Parameters
EntityHandlemeshset
BitRefLevelbitLevel
verbverbosity level
Examples
bernstein_bezier_generate_base.cpp.

Definition at line 197 of file MeshRefinement.cpp.

199 {
200 Interface &m_field = cOre;
201 moab::Interface &moab = m_field.get_moab();
203 Range tets;
204 CHKERR moab.get_entities_by_type(meshset, MBTET, tets, false);
205 CHKERR refineTets(tets, bit, verb, debug);
207}
static const bool debug
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset

◆ refineTets() [2/3]

MoFEMErrorCode MoFEM::MeshRefinement::refineTets ( const Range & tets,
const BitRefLevel & bit,
int verb = QUIET,
const bool debug = false )

refine TET in the meshset

Parameters
Rangeof tets to refine
BitRefLevelbitLevel
BitRefLevelbitLevel
verbverbosity level

Definition at line 244 of file MeshRefinement.cpp.

246 {
247
249
250 auto set_edge_bits = [](moab::Interface &moab,
252 &ref_parent_ents_view,
253 EntityHandle tet, BitRefEdges &parent_edges_bit,
254 EntityHandle *edge_new_nodes, int *split_edges) {
256
257 for (int ee = 0; ee != 6; ++ee) {
258 EntityHandle edge;
259 CHKERR moab.side_element(tet, 1, ee, edge);
260 auto eit = ref_parent_ents_view.find(edge);
261 if (eit != ref_parent_ents_view.end()) {
262 edge_new_nodes[ee] = (*eit)->getEnt();
263 split_edges[parent_edges_bit.count()] = ee;
264 parent_edges_bit.set(ee, 1);
265 }
266 }
267
269 };
270
271 CHKERR refineTets(_tets, bit, set_edge_bits, verb, debug);
272
274}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...

◆ refineTets() [3/3]

MoFEMErrorCode MoFEM::MeshRefinement::refineTets ( const Range & tets,
const BitRefLevel & bit,
SetEdgeBitsFun set_edge_bits,
int verb,
const bool debug )
private

refine TET in the meshset

Parameters
Rangeof tets to refine
BitRefLevelbitLevel
verbverbosity level

Definition at line 325 of file MeshRefinement.cpp.

328 {
329
330 Interface &m_field = cOre;
331 moab::Interface &moab = m_field.get_moab();
332 auto refined_ents_ptr = m_field.get_ref_ents();
333 auto refined_finite_elements_ptr = m_field.get_ref_finite_elements();
334 ReadUtilIface *read_util;
336
337 CHKERR m_field.get_moab().query_interface(read_util);
338
339 Range ents_to_set_bit;
340
341 // FIXME: refinement is based on entity handlers, should work on global ids of
342 // nodes, this will allow parallelise algorithm in the future
343
344 auto get_parent_ents_view = [&](const auto type_child,
345 const auto type_parent) {
346 auto &ref_ents =
347 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
348 auto range =
349 ref_ents.equal_range(boost::make_tuple(type_child, type_parent));
350
352
353 using I = decltype(range.first);
354
355 boost::function<bool(I it)> tester = [&](I it) {
356 return ((*it)->getBitRefLevel() & bit).any();
357 };
358
359 boost::function<MoFEMErrorCode(I f, I s)> inserter = [&](I f, I s) {
360 ref_parent_ents_view.insert(f, s);
361 return 0;
362 };
363
364 rangeInserter(range.first, range.second, tester, inserter);
365
366 return ref_parent_ents_view;
367 };
368
369 auto ref_parent_ents_view = get_parent_ents_view(MBVERTEX, MBEDGE);
370
371 auto tets = _tets.subset_by_type(MBTET);
372
373 auto get_ele_parent_view = [&]() {
374 auto &ref_ents = refined_finite_elements_ptr->get<Ent_mi_tag>();
375 RefElement_multiIndex_parents_view ref_ele_parent_view;
376 for (auto p = tets.pair_begin(); p != tets.pair_end(); ++p) {
377 auto tmi = ref_ents.lower_bound(p->first);
378 auto tme = ref_ents.upper_bound(p->second);
379 ref_ele_parent_view.insert(tmi, tme);
380 }
381 return ref_ele_parent_view;
382 };
383
384 auto ref_ele_parent_view = get_ele_parent_view();
385
386 std::vector<EntityHandle> parent_tets_refinded;
387 std::vector<EntityHandle> vertices_of_created_tets;
388 std::vector<BitRefEdges> parent_edges_bit_vec;
389 std::vector<int> nb_new_tets_vec;
390 std::vector<int> sub_type_vec;
391
392 parent_tets_refinded.reserve(tets.size());
393 vertices_of_created_tets.reserve(4 * tets.size());
394 parent_edges_bit_vec.reserve(tets.size());
395 nb_new_tets_vec.reserve(tets.size());
396 sub_type_vec.reserve(tets.size());
397
398 // make loop over all tets which going to be refined
399 for (auto tit : tets) {
400
401 // get tet connectivity
402 const EntityHandle *conn;
403 int num_nodes;
404 CHKERR moab.get_connectivity(tit, conn, num_nodes, true);
405
406 // get edges
407 BitRefEdges parent_edges_bit(0);
408 EntityHandle edge_new_nodes[] = {0, 0, 0, 0, 0, 0};
409 int split_edges[] = {-1, -1, -1, -1, -1, -1};
410
411 CHKERR set_edge_bits(moab, ref_parent_ents_view, tit, parent_edges_bit,
412 edge_new_nodes, split_edges);
413
414 if (parent_edges_bit.count()) {
415
416 // swap nodes forward
417 EntityHandle _conn_[4];
418 std::copy(&conn[0], &conn[4], &_conn_[0]);
419
420 // build connectivity for ref tets
421 EntityHandle new_tets_conns[8 * 4];
422 std::fill(&new_tets_conns[0], &new_tets_conns[8 * 4], no_handle);
423
424 int sub_type = -1, nb_new_tets = 0;
425 switch (parent_edges_bit.count()) {
426 case 1:
427 sub_type = 0;
428 tet_type_1(_conn_, split_edges[0], edge_new_nodes[split_edges[0]],
429 new_tets_conns);
430 nb_new_tets = 2;
431 break;
432 case 2:
433 sub_type =
434 tet_type_2(_conn_, split_edges, edge_new_nodes, new_tets_conns);
435 if (sub_type & (4 | 8 | 16)) {
436 nb_new_tets = 3;
437 break;
438 } else if (sub_type == 1) {
439 nb_new_tets = 4;
440 break;
441 };
442 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Impossible case");
443 break;
444 case 3:
445 sub_type =
446 tet_type_3(_conn_, split_edges, edge_new_nodes, new_tets_conns);
447 if (sub_type <= 4) {
448 nb_new_tets = 4;
449 break;
450 } else if (sub_type <= 7) {
451 nb_new_tets = 5;
452 break;
453 }
454 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Impossible case");
455 case 4:
456 sub_type =
457 tet_type_4(_conn_, split_edges, edge_new_nodes, new_tets_conns);
458 if (sub_type == 0) {
459 nb_new_tets = 5;
460 break;
461 } else if (sub_type <= 7) {
462 nb_new_tets = 6;
463 break;
464 }
465 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Impossible case");
466 case 5:
467 sub_type = tet_type_5(moab, _conn_, edge_new_nodes, new_tets_conns);
468 nb_new_tets = 7;
469 break;
470 case 6:
471 sub_type = 0;
472 tet_type_6(moab, _conn_, edge_new_nodes, new_tets_conns);
473 nb_new_tets = 8;
474 break;
475 default:
476 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Impossible case");
477 }
478
479 // find that tets
480 auto &ref_ele_by_parent_and_ref_edges =
481 ref_ele_parent_view
482 .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
483 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.equal_range(
484 boost::make_tuple(tit, parent_edges_bit.to_ulong()));
485 // check if tet with this refinement scheme already exits
486 std::vector<EntityHandle> ents_to_set_bit_vec;
487 if (std::distance(it_by_ref_edges.first, it_by_ref_edges.second) ==
488 static_cast<size_t>(nb_new_tets)) {
489 ents_to_set_bit_vec.reserve(nb_new_tets);
490 for (int tt = 0; it_by_ref_edges.first != it_by_ref_edges.second;
491 it_by_ref_edges.first++, tt++) {
492 auto tet = (*it_by_ref_edges.first)->getEnt();
493 ents_to_set_bit_vec.emplace_back(tet);
494 // set ref tets entities
495 if (debug) {
496 // add this tet if exist to this ref level
497 auto ref_tet_it = refined_ents_ptr->find(tet);
498 if (ref_tet_it == refined_ents_ptr->end()) {
499 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
500 "Tetrahedron should be in database");
501 }
502 }
503 }
504 ents_to_set_bit.insert_list(ents_to_set_bit_vec.begin(),
505 ents_to_set_bit_vec.end());
506
507 } else {
508 // if this element was not refined or was refined with different
509 // patterns of split edges create new elements
510
511 parent_tets_refinded.emplace_back(tit);
512 for (int tt = 0; tt != nb_new_tets; ++tt) {
513 for (auto nn : {0, 1, 2, 3})
514 vertices_of_created_tets.emplace_back(new_tets_conns[4 * tt + nn]);
515 }
516 parent_edges_bit_vec.emplace_back(parent_edges_bit);
517 nb_new_tets_vec.emplace_back(nb_new_tets);
518 sub_type_vec.emplace_back(sub_type);
519 }
520 } else {
521
522 if (debug) {
523 RefEntity_multiIndex::iterator tit_miit;
524 tit_miit = refined_ents_ptr->find(tit);
525 if (tit_miit == refined_ents_ptr->end())
526 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
527 "can not find this tet");
528 }
529
530 ents_to_set_bit.insert(tit);
531 }
532 }
533
534 // Create tets
535 EntityHandle start_e = 0;
536 EntityHandle *conn_moab;
537 const int nb_new_tets = vertices_of_created_tets.size() / 4;
538 read_util->get_element_connect(nb_new_tets, 4, MBTET, 0, start_e, conn_moab);
539 std::copy(vertices_of_created_tets.begin(), vertices_of_created_tets.end(),
540 conn_moab);
541 CHKERR read_util->update_adjacencies(start_e, nb_new_tets, 4, conn_moab);
542 ents_to_set_bit.insert(start_e, start_e + nb_new_tets - 1);
543
544 // Create adj entities
545 Range ents;
546 for (auto d : {1, 2}) {
547 CHKERR moab.get_adjacencies(ents_to_set_bit, d, true, ents,
548 moab::Interface::UNION);
549 }
550
551 SetParent set_parent;
552
553 // Set parents and adjacencies
554 for (int idx = 0; idx != parent_tets_refinded.size(); ++idx) {
555
556 const EntityHandle tit = parent_tets_refinded[idx];
557 const BitRefEdges &parent_edges_bit = parent_edges_bit_vec[idx];
558 const int nb_new_tets = nb_new_tets_vec[idx];
559
560 std::array<EntityHandle, 8> ref_tets;
561 for (int tt = 0; tt != nb_new_tets; ++tt, ++start_e)
562 ref_tets[tt] = start_e;
563
564 if (nb_new_tets) {
565
566 CHKERR moab.tag_clear_data(cOre.get_th_RefBitEdge(), &ref_tets[0],
567 nb_new_tets, &parent_edges_bit);
568 CHKERR moab.tag_clear_data(cOre.get_th_RefParentHandle(), &ref_tets[0],
569 nb_new_tets, &tit);
570
571 // hash map of nodes (RefEntity) by edges (EntityHandle)
572 std::map<EntityHandle /*edge*/, EntityHandle /*node*/>
573 map_ref_nodes_by_edges;
574
575 Range tit_edges;
576 CHKERR moab.get_adjacencies(&tit, 1, 1, false, tit_edges);
577 for (auto edge : tit_edges) {
578 auto eit = ref_parent_ents_view.find(edge);
579 if (eit != ref_parent_ents_view.end()) {
580 map_ref_nodes_by_edges[(*eit)->getParentEnt()] = eit->get()->getEnt();
581 }
582 }
583
584 // find parents for new edges and faces
585 // get tet edges and faces
586 Range tit_faces;
587 CHKERR moab.get_adjacencies(&tit, 1, 2, false, tit_faces);
588 Range edges_nodes[6], faces_nodes[4];
589 // for edges - add ref nodes
590 // edges_nodes[ee] - contains all nodes on edge ee including mid nodes if
591 // exist
592 Range::iterator eit = tit_edges.begin();
593 for (int ee = 0; eit != tit_edges.end(); eit++, ee++) {
594 CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee], true);
595 auto map_miit = map_ref_nodes_by_edges.find(*eit);
596 if (map_miit != map_ref_nodes_by_edges.end()) {
597 edges_nodes[ee].insert(map_miit->second);
598 }
599 }
600 // for faces - add ref nodes
601 // faces_nodes[ff] - contains all nodes on face ff including mid nodes if
602 // exist
603 Range::iterator fit = tit_faces.begin();
604 for (int ff = 0; fit != tit_faces.end(); fit++, ff++) {
605 CHKERR moab.get_connectivity(&*fit, 1, faces_nodes[ff], true);
606 // Get edges on face and loop over those edges to add mid-nodes to range
607 Range fit_edges;
608 CHKERR moab.get_adjacencies(&*fit, 1, 1, false, fit_edges);
609 for (Range::iterator eit2 = fit_edges.begin(); eit2 != fit_edges.end();
610 eit2++) {
611 auto map_miit = map_ref_nodes_by_edges.find(*eit2);
612 if (map_miit != map_ref_nodes_by_edges.end()) {
613 faces_nodes[ff].insert(map_miit->second);
614 }
615 }
616 }
617 // add ref nodes to tet
618 // tet_nodes contains all nodes on tet including mid edge nodes
619 Range tet_nodes;
620 CHKERR moab.get_connectivity(&tit, 1, tet_nodes, true);
621 for (std::map<EntityHandle, EntityHandle>::iterator map_miit =
622 map_ref_nodes_by_edges.begin();
623 map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
624 tet_nodes.insert(map_miit->second);
625 }
626 Range ref_edges;
627 // Get all all edges of refined tets
628 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 1, false,
629 ref_edges, moab::Interface::UNION);
630 // Check for all ref edge and set parents
631 for (Range::iterator reit = ref_edges.begin(); reit != ref_edges.end();
632 reit++) {
633 Range ref_edges_nodes;
634 CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes, true);
635
636 // Check if ref edge is an coarse edge (loop over coarse tet edges)
637 int ee = 0;
638 for (; ee < 6; ee++) {
639 // Two nodes are common (node[0],node[1],ref_node (if exist))
640 // this tests if given edge is contained by edge of refined
641 // tetrahedral
642 if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
643 EntityHandle edge = tit_edges[ee];
644 CHKERR set_parent(*reit, edge, refined_ents_ptr, cOre);
645 break;
646 }
647 }
648 if (ee < 6)
649 continue; // this refined edge is contained by edge of tetrahedral
650 // check if ref edge is in coarse face
651 int ff = 0;
652 for (; ff < 4; ff++) {
653 // two nodes are common (node[0],node[1],ref_node (if exist))
654 // this tests if given face is contained by face of tetrahedral
655 if (intersect(faces_nodes[ff], ref_edges_nodes).size() == 2) {
656 EntityHandle face = tit_faces[ff];
657 CHKERR set_parent(*reit, face, refined_ents_ptr, cOre);
658 break;
659 }
660 }
661 if (ff < 4)
662 continue; // this refined edge is contained by face of tetrahedral
663
664 // check if ref edge is in coarse tetrahedral (i.e. that is internal
665 // edge of refined tetrahedral)
666 if (intersect(tet_nodes, ref_edges_nodes).size() == 2) {
667 CHKERR set_parent(*reit, tit, refined_ents_ptr, cOre);
668 continue;
669 }
670
671 // Refined edge is not child of any edge, face or tetrahedral, this is
672 // Impossible edge
673 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
674 "impossible refined edge");
675 }
676
677 Range ref_faces;
678 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 2, false,
679 ref_faces, moab::Interface::UNION);
680 Tag th_interface_side;
681 const int def_side[] = {0};
682 CHKERR moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
683 th_interface_side,
684 MB_TAG_CREAT | MB_TAG_SPARSE, def_side);
685 // Check for all ref faces
686 for (auto rfit : ref_faces) {
687 Range ref_faces_nodes;
688 CHKERR moab.get_connectivity(&rfit, 1, ref_faces_nodes, true);
689 // Check if ref face is in coarse face
690 int ff = 0;
691 for (; ff < 4; ff++) {
692 // Check if refined triangle is contained by face of tetrahedral
693 if (intersect(faces_nodes[ff], ref_faces_nodes).size() == 3) {
694 EntityHandle face = tit_faces[ff];
695 CHKERR set_parent(rfit, face, refined_ents_ptr, cOre);
696 int side = 0;
697 // Set face side if it is on interface
698 CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
699 CHKERR moab.tag_set_data(th_interface_side, &rfit, 1, &side);
700 break;
701 }
702 }
703 if (ff < 4)
704 continue; // this face is contained by one of tetrahedrons
705 // check if ref face is in coarse tetrahedral
706 // this is ref face which is contained by tetrahedral volume
707 if (intersect(tet_nodes, ref_faces_nodes).size() == 3) {
708 CHKERR set_parent(rfit, tit, refined_ents_ptr, cOre);
709 continue;
710 }
711 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
712 "impossible refined face");
713 }
714 }
715 }
716
717 CHKERR set_parent(refined_ents_ptr);
718 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
719 ents_to_set_bit, bit, true, verb, &ents);
721}
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
int tet_type_3(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
void tet_type_6(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
int tet_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
int tet_type_5(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
int tet_type_4(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
void tet_type_1(const EntityHandle *conn, const int split_edge, const EntityHandle edge_new_node, EntityHandle *new_tets_conn)
auto rangeInserter(const I f, const I s, boost::function< bool(I it)> tester, boost::function< MoFEMErrorCode(I f, I s)> inserter)
Insert ranges.
constexpr IntegrationType I

◆ refineTetsHangingNodes() [1/2]

MoFEMErrorCode MoFEM::MeshRefinement::refineTetsHangingNodes ( const EntityHandle meshset,
const BitRefLevel & bit,
int verb = QUIET,
const bool debug = false )

refine TET in the meshset

Parameters
Rangeof tets to refine
BitRefLevelbitLevel
BitRefLevelbitLevel
verbverbosity level

Definition at line 313 of file MeshRefinement.cpp.

315 {
316 Interface &m_field = cOre;
317 moab::Interface &moab = m_field.get_moab();
319 Range tets;
320 CHKERR moab.get_entities_by_type(meshset, MBTET, tets, false);
323}
MoFEMErrorCode refineTetsHangingNodes(const Range &tets, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset

◆ refineTetsHangingNodes() [2/2]

MoFEMErrorCode MoFEM::MeshRefinement::refineTetsHangingNodes ( const Range & tets,
const BitRefLevel & bit,
int verb = QUIET,
const bool debug = false )

refine TET in the meshset

Parameters
Rangeof tets to refine
BitRefLevelbitLevel
BitRefLevelbitLevel
verbverbosity level

Definition at line 276 of file MeshRefinement.cpp.

279 {
280
282
283 auto set_edge_bits = [](moab::Interface &moab,
285 &ref_parent_ents_view,
286 EntityHandle tet, BitRefEdges &parent_edges_bit,
287 EntityHandle *edge_new_nodes, int *split_edges) {
289
290 for (int ee = 0; ee != 6; ++ee) {
291 EntityHandle edge;
292 CHKERR moab.side_element(tet, 1, ee, edge);
293 auto eit = ref_parent_ents_view.find(edge);
294 if (eit != ref_parent_ents_view.end()) {
295 edge_new_nodes[ee] = (*eit)->getEnt();
296 split_edges[parent_edges_bit.count()] = ee;
297 parent_edges_bit.set(ee, 1);
298 }
299 }
300
301 if(parent_edges_bit.count() < 6)
302 parent_edges_bit.reset();
303
305 };
306
307 CHKERR refineTets(_tets, bit, set_edge_bits, verb, debug);
308
310}

◆ refineTris() [1/3]

MoFEMErrorCode MoFEM::MeshRefinement::refineTris ( const EntityHandle meshset,
const BitRefLevel & bit,
int verb = QUIET,
const bool debug = false )

refine triangles in the meshset

Parameters
EntityHandlemeshset
BitRefLevelbitLevel
verbverbosity level

Definition at line 922 of file MeshRefinement.cpp.

924 {
925 Interface &m_field = cOre;
926 moab::Interface &moab = m_field.get_moab();
928 Range tris;
929 CHKERR moab.get_entities_by_type(meshset, MBTRI, tris, false);
930 CHKERR refineTris(tris, bit, verb, debug);
932}
MoFEMErrorCode refineTris(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine triangles in the meshset

◆ refineTris() [2/3]

MoFEMErrorCode MoFEM::MeshRefinement::refineTris ( const Range & tris,
const BitRefLevel & bit,
int verb = QUIET,
const bool debug = false )

refine TRI in the meshset

Parameters
meshsetof entities to refine
BitRefLevelbit level of created entities
verbverbosity level

Definition at line 934 of file MeshRefinement.cpp.

936 {
938
939 auto set_edge_bits = [](moab::Interface &moab,
941 &ref_parent_ents_view,
942 EntityHandle tet, BitRefEdges &parent_edges_bit,
943 EntityHandle *edge_new_nodes, int *split_edges) {
945
946 for (int ee = 0; ee != 3; ++ee) {
947 EntityHandle edge;
948 CHKERR moab.side_element(tet, 1, ee, edge);
949 auto eit = ref_parent_ents_view.find(edge);
950 if (eit != ref_parent_ents_view.end()) {
951 edge_new_nodes[ee] = (*eit)->getEnt();
952 split_edges[parent_edges_bit.count()] = ee;
953 parent_edges_bit.set(ee, 1);
954 }
955 }
956
958 };
959
960 CHKERR refineTris(_tris, bit, set_edge_bits, verb, debug);
961
963}

◆ refineTris() [3/3]

MoFEMErrorCode MoFEM::MeshRefinement::refineTris ( const Range & tris,
const BitRefLevel & bit,
SetEdgeBitsFun set_edge_bits,
int verb,
const bool debug )
private

Definition at line 1015 of file MeshRefinement.cpp.

1018 {
1019 Interface &m_field = cOre;
1020 moab::Interface &moab = m_field.get_moab();
1021 auto refined_ents_ptr = m_field.get_ref_ents();
1022 auto refined_finite_elements_ptr = m_field.get_ref_finite_elements();
1023 ReadUtilIface *read_util;
1024
1026
1027 auto get_parent_ents_view = [&](const auto type_child,
1028 const auto type_parent) {
1029 auto &ref_ents =
1030 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
1031 auto range =
1032 ref_ents.equal_range(boost::make_tuple(type_child, type_parent));
1033
1035
1036 using I = decltype(range.first);
1037
1038 boost::function<bool(I it)> tester = [&](I it) {
1039 return ((*it)->getBitRefLevel() & bit).any();
1040 };
1041
1042 boost::function<MoFEMErrorCode(I f, I s)> inserter = [&](I f, I s) {
1043 ref_parent_ents_view.insert(f, s);
1044 return 0;
1045 };
1046
1047 rangeInserter(range.first, range.second, tester, inserter);
1048
1049 return ref_parent_ents_view;
1050 };
1051
1052 auto ref_parent_ents_view = get_parent_ents_view(MBVERTEX, MBEDGE);
1053
1054 auto tris = _tris.subset_by_type(MBTRI);
1055
1056 auto get_ele_parent_view = [&]() {
1057 auto &ref_ents = refined_finite_elements_ptr->get<Ent_mi_tag>();
1058 RefElement_multiIndex_parents_view ref_ele_parent_view;
1059 for (auto p = tris.pair_begin(); p != tris.pair_end(); ++p) {
1060 auto tmi = ref_ents.lower_bound(p->first);
1061 auto tme = ref_ents.upper_bound(p->second);
1062 ref_ele_parent_view.insert(tmi, tme);
1063 }
1064 return ref_ele_parent_view;
1065 };
1066
1067 auto ref_ele_parent_view = get_ele_parent_view();
1068
1069 Range ents_to_set_bit;
1070
1071 std::vector<EntityHandle> parent_tris_refined; // entities refined
1072 std::vector<EntityHandle> vertices_of_created_tris;
1073 std::vector<BitRefEdges> parent_edges_bit_vec;
1074 std::vector<int> nb_new_tris_vec;
1075
1076 parent_tris_refined.reserve(tris.size());
1077 vertices_of_created_tris.reserve(3 * tris.size());
1078 parent_edges_bit_vec.reserve(tris.size());
1079 nb_new_tris_vec.reserve(tris.size());
1080
1081 // make loop over all tets which going to be refined
1082 for (auto tit : tris) {
1083
1084 // get tet connectivity
1085 const EntityHandle *conn;
1086 int num_nodes;
1087 CHKERR moab.get_connectivity(tit, conn, num_nodes, true);
1088
1089 // get edges
1090 BitRefEdges parent_edges_bit(0);
1091 EntityHandle edge_new_nodes[] = {0, 0, 0};
1092 int split_edges[] = {-1, -1, -1};
1093
1094 CHKERR set_edge_bits(moab, ref_parent_ents_view, tit, parent_edges_bit,
1095 edge_new_nodes, split_edges);
1096
1097 if (parent_edges_bit.count()) {
1098
1099 // swap nodes forward
1100 EntityHandle _conn_[3];
1101 std::copy(&conn[0], &conn[3], &_conn_[0]);
1102
1103 // build connectivity for ref tets
1104 EntityHandle new_tets_conns[4 * 3];
1105 std::fill(&new_tets_conns[0], &new_tets_conns[4 * 3], no_handle);
1106
1107 int nb_new_tris = 0;
1108 switch (parent_edges_bit.count()) {
1109 case 1:
1110 CHKERR tri_type_1(conn, split_edges[0], edge_new_nodes[split_edges[0]],
1111 new_tets_conns);
1112 nb_new_tris += 2;
1113 break;
1114 case 2:
1115 CHKERR tri_type_2(conn, split_edges, edge_new_nodes, new_tets_conns);
1116 nb_new_tris += 3;
1117 break;
1118 case 3:
1119 CHKERR tri_type_3(conn, edge_new_nodes, new_tets_conns);
1120 nb_new_tris += 4;
1121 break;
1122 default:
1123 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Impossible case");
1124 }
1125
1126 // find that tets
1127 auto &ref_ele_by_parent_and_ref_edges =
1128 ref_ele_parent_view
1129 .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
1130 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.equal_range(
1131 boost::make_tuple(tit, parent_edges_bit.to_ulong()));
1132 // check if tet with this refinement scheme already exits
1133 std::vector<EntityHandle> ents_to_set_bit_vec;
1134 if (std::distance(it_by_ref_edges.first, it_by_ref_edges.second) ==
1135 static_cast<size_t>(nb_new_tris)) {
1136 ents_to_set_bit_vec.reserve(nb_new_tris);
1137 for (int tt = 0; it_by_ref_edges.first != it_by_ref_edges.second;
1138 it_by_ref_edges.first++, tt++) {
1139 auto tet = (*it_by_ref_edges.first)->getEnt();
1140 ents_to_set_bit_vec.emplace_back(tet);
1141 // set ref tets entities
1142 if (debug) {
1143 // add this tet if exist to this ref level
1144 auto ref_tet_it = refined_ents_ptr->find(tet);
1145 if (ref_tet_it == refined_ents_ptr->end()) {
1146 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1147 "Tetrahedron should be in database");
1148 }
1149 }
1150 }
1151 ents_to_set_bit.insert_list(ents_to_set_bit_vec.begin(),
1152 ents_to_set_bit_vec.end());
1153
1154 } else {
1155 // if this element was not refined or was refined with different
1156 // patterns of split edges create new elements
1157
1158 parent_tris_refined.emplace_back(tit);
1159 for (int tt = 0; tt != nb_new_tris; ++tt) {
1160 for (auto nn : {0, 1, 2})
1161 vertices_of_created_tris.emplace_back(new_tets_conns[3 * tt + nn]);
1162 }
1163 parent_edges_bit_vec.emplace_back(parent_edges_bit);
1164 nb_new_tris_vec.emplace_back(nb_new_tris);
1165 }
1166 } else {
1167
1168 if (debug) {
1169 RefEntity_multiIndex::iterator tit_miit;
1170 tit_miit = refined_ents_ptr->find(tit);
1171 if (tit_miit == refined_ents_ptr->end())
1172 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1173 "can not find this tet");
1174 }
1175
1176 ents_to_set_bit.insert(tit);
1177 }
1178 }
1179
1180 const int nb_new_tris = vertices_of_created_tris.size() / 3;
1181 EntityHandle start_e = 0;
1182
1183 if (nb_new_tris) {
1184
1185 // Create tets
1186 EntityHandle *conn_moab;
1187 CHKERR m_field.get_moab().query_interface(read_util);
1188 CHKERR read_util->get_element_connect(nb_new_tris, 3, MBTRI, 0, start_e,
1189 conn_moab);
1190 std::copy(vertices_of_created_tris.begin(), vertices_of_created_tris.end(),
1191 conn_moab);
1192 CHKERR read_util->update_adjacencies(start_e, nb_new_tris, 3, conn_moab);
1193 ents_to_set_bit.insert(start_e, start_e + nb_new_tris - 1);
1194
1195 if (debug) {
1196 auto meshset_ptr = get_temp_meshset_ptr(moab);
1197 CHKERR moab.add_entities(*meshset_ptr, ents_to_set_bit);
1198 CHKERR moab.write_file("out_refinded_debug.vtk", "VTK", "",
1199 meshset_ptr->get_ptr(), 1);
1200 }
1201 }
1202
1203 // Create adj entities
1204 Range ents;
1205 for (auto d : {1}) {
1206 CHKERR moab.get_adjacencies(ents_to_set_bit, d, true, ents,
1207 moab::Interface::UNION);
1208 }
1209
1210 SetParent set_parent;
1211
1212 // Set parrents and adjacencies. Iterate over all parent triangles.
1213 for (int idx = 0; idx != parent_tris_refined.size(); ++idx) {
1214
1215 const auto tit = parent_tris_refined[idx]; // parent triangle
1216 const auto &parent_edges_bit =
1217 parent_edges_bit_vec[idx]; // marks what edges are refined
1218 const auto nb_new_tris =
1219 nb_new_tris_vec[idx]; // number of triangles created by refinement
1220
1221 // Iterate over refined triangles
1222 std::array<EntityHandle, 4> ref_tris;
1223 for (int tt = 0; tt != nb_new_tris; ++tt, ++start_e)
1224 ref_tris[tt] = start_e;
1225
1226 if (nb_new_tris) {
1227
1228 CHKERR moab.tag_clear_data(cOre.get_th_RefBitEdge(), &ref_tris[0],
1229 nb_new_tris, &parent_edges_bit);
1230 CHKERR moab.tag_clear_data(cOre.get_th_RefParentHandle(), &ref_tris[0],
1231 nb_new_tris, &tit);
1232
1233 // Hash map of nodes (RefEntity) by edges (EntityHandle)
1234 std::map<EntityHandle /*edge*/, EntityHandle /*node*/>
1235 map_ref_nodes_by_edges;
1236
1237 // find parents for new edges and faces
1238 // get tet edges and faces
1239 Range tit_edges; // parent edges
1240 CHKERR moab.get_adjacencies(&tit, 1, 1, false, tit_edges);
1241 for (auto edge : tit_edges) {
1242 auto eit = ref_parent_ents_view.find(edge);
1243 if (eit != ref_parent_ents_view.end()) {
1244 map_ref_nodes_by_edges[(*eit)->getParentEnt()] = eit->get()->getEnt();
1245 }
1246 }
1247
1248 Range edges_nodes[3];
1249 // for edges - add ref nodes
1250 // edges_nodes[ee] - contains all nodes on edge ee including mid nodes if
1251 // exist
1252 Range::iterator eit = tit_edges.begin();
1253 for (int ee = 0; eit != tit_edges.end(); eit++, ee++) {
1254 CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee], true);
1255 auto map_miit = map_ref_nodes_by_edges.find(*eit);
1256 if (map_miit != map_ref_nodes_by_edges.end()) {
1257 edges_nodes[ee].insert(map_miit->second);
1258 }
1259 }
1260
1261 // add ref nodes to tri
1262 // tet_nodes contains all nodes on tet including mid edge nodes
1263 Range tri_nodes;
1264 CHKERR moab.get_connectivity(&tit, 1, tri_nodes, true);
1265 for (auto map_miit = map_ref_nodes_by_edges.begin();
1266 map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
1267 tri_nodes.insert(map_miit->second);
1268 }
1269
1270 Range ref_edges; // edges on all new tets
1271 // Get all all edges of refined tets
1272 CHKERR moab.get_adjacencies(ref_tris.data(), nb_new_tris, 1, false,
1273 ref_edges, moab::Interface::UNION);
1274
1275 // Check for all ref edge and set parents
1276 for (auto reit = ref_edges.begin(); reit != ref_edges.end(); ++reit) {
1277 Range ref_edges_nodes;
1278 CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes, true);
1279 // Check if ref edge is an coarse edge (loop over coarse tet edges)
1280 int ee = 0;
1281 for (; ee < 3; ee++) {
1282 // Two nodes are common (node[0],node[1],ref_node (if exist))
1283 // this tests if given edge is contained by edge of refined
1284 // triangle
1285 if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
1286 EntityHandle edge = tit_edges[ee];
1287 CHKERR set_parent(*reit, edge, refined_ents_ptr, cOre);
1288 break;
1289 }
1290 }
1291 if (ee < 3)
1292 continue; // this refined edge is contained by edge of triangle
1293
1294 // check if ref edge is in coarse tetrahedral (i.e. that is internal
1295 // edge of refined tetrahedral)
1296 if (intersect(tri_nodes, ref_edges_nodes).size() == 2) {
1297 CHKERR set_parent(*reit, tit, refined_ents_ptr, cOre);
1298 continue;
1299 }
1300
1301 // Refined edge is not child of any edge, face or tetrahedral, this is
1302 // Impossible edge
1303 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1304 "impossible refined edge");
1305 }
1306 }
1307 }
1308
1309 CHKERR set_parent(refined_ents_ptr);
1310 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
1311 ents_to_set_bit, bit, true, verb, &ents);
1312
1314}
MoFEMErrorCode tri_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tris_conn)
MoFEMErrorCode tri_type_1(const EntityHandle *conn, const int split_edge, const EntityHandle edge_new_node, EntityHandle *new_tris_conn)
MoFEMErrorCode tri_type_3(const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tris_conn)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.

◆ refineTrisHangingNodes() [1/2]

MoFEMErrorCode MoFEM::MeshRefinement::refineTrisHangingNodes ( const EntityHandle meshset,
const BitRefLevel & bit,
int verb = QUIET,
const bool debug = false )

refine TRI in the meshset

Parameters
Rangeof entities to refine
BitRefLevelbit level of created entities
verbverbosity level

Definition at line 1003 of file MeshRefinement.cpp.

1005 {
1006 Interface &m_field = cOre;
1007 moab::Interface &moab = m_field.get_moab();
1009 Range tris;
1010 CHKERR moab.get_entities_by_type(meshset, MBTRI, tris, false);
1011 CHKERR refineTrisHangingNodes(tris, bit, verb, debug);
1013}
MoFEMErrorCode refineTrisHangingNodes(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TRI in the meshset

◆ refineTrisHangingNodes() [2/2]

MoFEMErrorCode MoFEM::MeshRefinement::refineTrisHangingNodes ( const Range & tris,
const BitRefLevel & bit,
int verb = QUIET,
const bool debug = false )

refine TRI in the meshset

Parameters
Rangeof entities to refine
BitRefLevelbit level of created entities
verbverbosity level

Definition at line 965 of file MeshRefinement.cpp.

968 {
969
971
972 auto set_edge_bits = [](moab::Interface &moab,
974 &ref_parent_ents_view,
975 EntityHandle tet, BitRefEdges &parent_edges_bit,
976 EntityHandle *edge_new_nodes, int *split_edges) {
978
979 for (int ee = 0; ee != 3; ++ee) {
980 EntityHandle edge;
981 CHKERR moab.side_element(tet, 1, ee, edge);
982 auto eit = ref_parent_ents_view.find(edge);
983 if (eit != ref_parent_ents_view.end()) {
984 edge_new_nodes[ee] = (*eit)->getEnt();
985 split_edges[parent_edges_bit.count()] = ee;
986 parent_edges_bit.set(ee, 1);
987 }
988 }
989
990 if (parent_edges_bit.count() < 3) {
991 parent_edges_bit.reset();
992 }
993
995 };
996
997 CHKERR refineTris(_tris, bit, set_edge_bits, verb, debug);
998
1000};

Member Data Documentation

◆ cOre

MoFEM::Core& MoFEM::MeshRefinement::cOre

Definition at line 31 of file MeshRefinement.hpp.


The documentation for this struct was generated from the following files: