v0.14.0
Loading...
Searching...
No Matches
Classes | Public Member Functions | Public Attributes | Private Types | Private Member Functions | List of all members
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 More...
 
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 More...
 
MoFEMErrorCode refineTets (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TET in the meshset More...
 
MoFEMErrorCode refineTets (const Range &tets, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TET in the meshset More...
 
MoFEMErrorCode refineTetsHangingNodes (const Range &tets, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TET in the meshset More...
 
MoFEMErrorCode refineTetsHangingNodes (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TET in the meshset More...
 
MoFEMErrorCode refinePrisms (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET)
 refine PRISM in the meshset More...
 
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 More...
 
MoFEMErrorCode refineTris (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine triangles in the meshset More...
 
MoFEMErrorCode refineTris (const Range &tris, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TRI in the meshset More...
 
MoFEMErrorCode refineTrisHangingNodes (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TRI in the meshset More...
 
MoFEMErrorCode refineTrisHangingNodes (const Range &tris, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
 refine TRI in the meshset More...
 
- Public Member Functions inherited from MoFEM::UnknownInterface
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

Public Attributes

MoFEM::CorecOre
 

Private Types

using SetEdgeBitsFun = boost::function< MoFEMErrorCode(moab::Interface &moab, RefEntity_multiIndex_view_by_ordered_parent_entity &ref_parent_ents_view, EntityHandle tet, BitRefEdges &parent_edges_bit, EntityHandle *edge_new_nodes, int *split_edges)>
 Functions setting edges for refinemnt on enetity level. More...
 

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 More...
 
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. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

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
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

using MoFEM::MeshRefinement::SetEdgeBitsFun = boost::function< MoFEMErrorCode(moab::Interface &moab, RefEntity_multiIndex_view_by_ordered_parent_entity &ref_parent_ents_view, EntityHandle tet, BitRefEdges &parent_edges_bit, EntityHandle *edge_new_nodes, int *split_edges )>
private

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)) {}
CoreTmp< 0 > Core
Definition: Core.hpp:1094

◆ ~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.

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 ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
auto bit
set bit
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
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
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 SETERRQ1(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.
Definition: LogManager.hpp:362
@ VERBOSE
Definition: definitions.h:209
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
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
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 904 of file MeshRefinement.cpp.

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

◆ 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(
738 get_id_for_min_type<MBPRISM>()),
739 refined_finite_elements_ptr->get<Ent_mi_tag>().upper_bound(
740 get_id_for_max_type<MBPRISM>()));
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(), 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, ss.str().c_str());
830 }
831 // find that prism
832 std::bitset<4> ref_prism_bit(0);
833 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.lower_bound(
834 boost::make_tuple(*pit, split_edges.to_ulong()));
835 auto hi_it_by_ref_edges = ref_ele_by_parent_and_ref_edges.upper_bound(
836 boost::make_tuple(*pit, split_edges.to_ulong()));
837 auto it_by_ref_edges2 = it_by_ref_edges;
838 for (int pp = 0; it_by_ref_edges2 != hi_it_by_ref_edges;
839 it_by_ref_edges2++, pp++) {
840 // Add this tet to this ref
841 *(const_cast<RefElement *>(it_by_ref_edges2->get())
842 ->getBitRefLevelPtr()) |= bit;
843 ref_prism_bit.set(pp, 1);
844 if (verb > 2) {
845 std::ostringstream ss;
846 ss << "is refined " << *it_by_ref_edges2 << std::endl;
847 PetscPrintf(m_field.get_comm(), ss.str().c_str());
848 }
849 }
850 if (it_by_ref_edges != hi_it_by_ref_edges) {
851 if (ref_prism_bit.count() != (unsigned int)nb_new_prisms)
852 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
853 "data inconsistency");
854 } else {
855 EntityHandle ref_prisms[4];
856 // create prism
857 for (int pp = 0; pp < nb_new_prisms; pp++) {
858 if (verb > 3) {
859 std::ostringstream ss;
860 ss << "ref prism " << ref_prism_bit << std::endl;
861 PetscPrintf(m_field.get_comm(), ss.str().c_str());
862 }
863 if (!ref_prism_bit.test(pp)) {
864 CHKERR moab.create_element(MBPRISM, &new_prism_conn[6 * pp], 6,
865 ref_prisms[pp]);
866 CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(),
867 &ref_prisms[pp], 1, &*pit);
868 CHKERR moab.tag_set_data(cOre.get_th_RefBitLevel(), &ref_prisms[pp],
869 1, &bit);
870 CHKERR moab.tag_set_data(cOre.get_th_RefBitEdge(), &ref_prisms[pp], 1,
871 &split_edges);
872 auto p_ent =
873 const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
874 ->insert(boost::shared_ptr<RefEntity>(new RefEntity(
875 m_field.get_basic_entity_data_ptr(), ref_prisms[pp])));
876 std::pair<RefElement_multiIndex::iterator, bool> p_fe;
877 try {
878 p_fe =
879 const_cast<RefElement_multiIndex *>(refined_finite_elements_ptr)
880 ->insert(boost::shared_ptr<RefElement>(
881 new RefElement_PRISM(*p_ent.first)));
882 } catch (MoFEMException const &e) {
883 SETERRQ(PETSC_COMM_SELF, e.errorCode, e.errorMessage);
884 }
885 ref_prism_bit.set(pp);
886 CHKERR cOre.addPrismToDatabase(ref_prisms[pp]);
887 if (verb > 2) {
888 std::ostringstream ss;
889 ss << "add prism: " << **p_fe.first << std::endl;
890 if (verb > 7) {
891 for (int nn = 0; nn < 6; nn++) {
892 ss << new_prism_conn[nn] << " ";
893 }
894 ss << std::endl;
895 }
896 PetscPrintf(m_field.get_comm(), ss.str().c_str());
897 }
898 }
899 }
900 }
901 }
903}
@ NOISY
Definition: definitions.h:211
@ VERY_NOISY
Definition: definitions.h:212
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< 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< RefElement >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getEnt > > > > RefElement_multiIndex
std::bitset< BITREFEDGES_SIZE > BitRefEdges
Definition: Types.hpp:34
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:493
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

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()
Definition: definitions.h:447
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440

◆ 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}
static Index< 'p', 3 > p
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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.
Definition: Templates.hpp:1916
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 921 of file MeshRefinement.cpp.

923 {
924 Interface &m_field = cOre;
925 moab::Interface &moab = m_field.get_moab();
927 Range tris;
928 CHKERR moab.get_entities_by_type(meshset, MBTRI, tris, false);
929 CHKERR refineTris(tris, bit, verb, debug);
931}
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 933 of file MeshRefinement.cpp.

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

◆ 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 1014 of file MeshRefinement.cpp.

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

◆ 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 1002 of file MeshRefinement.cpp.

1004 {
1005 Interface &m_field = cOre;
1006 moab::Interface &moab = m_field.get_moab();
1008 Range tris;
1009 CHKERR moab.get_entities_by_type(meshset, MBTRI, tris, false);
1010 CHKERR refineTrisHangingNodes(tris, bit, verb, debug);
1012}
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 964 of file MeshRefinement.cpp.

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

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: