v0.13.2
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, 6> coords;
126 &coords[0], &coords[1], &coords[2]};
128 &coords[3], &coords[4], &coords[5]};
130
131 Range add_bit;
132 for (auto p_eit = edges.pair_begin(); p_eit != edges.pair_end(); ++p_eit) {
133 auto miit_view = ref_parent_ents_view.lower_bound(p_eit->first);
134
135 Range edge_having_parent_vertex;
136 if (miit_view != ref_parent_ents_view.end()) {
137 for (auto hi_miit_view = ref_parent_ents_view.upper_bound(p_eit->second);
138 miit_view != hi_miit_view; ++miit_view) {
139 edge_having_parent_vertex.insert(edge_having_parent_vertex.end(),
140 miit_view->get()->getParentEnt());
141 add_bit.insert(add_bit.end(), miit_view->get()->getEnt());
142 }
143 }
144
145 Range add_vertex_edges =
146 subtract(Range(p_eit->first, p_eit->second), edge_having_parent_vertex);
147
148 for (auto e : add_vertex_edges)
149 parent_edge.emplace_back(e);
150
151 for (auto e : add_vertex_edges) {
152
153 const EntityHandle *conn;
154 int num_nodes;
155 CHKERR moab.get_connectivity(e, conn, num_nodes, true);
156 if (PetscUnlikely(num_nodes != 2)) {
157 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
158 "edge should have 2 edges");
159 }
160 CHKERR moab.get_coords(conn, num_nodes, coords.data());
161 t_0(i) += t_1(i);
162 t_0(i) *= 0.5;
163
164 for (auto j : {0, 1, 2})
165 vert_coords[j].emplace_back(t_0(j));
166 }
167 }
168
169 CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit, bit);
170
171 if (!vert_coords[0].empty()) {
172 ReadUtilIface *read_util;
173 CHKERR moab.query_interface(read_util);
174 int num_nodes = vert_coords[0].size();
175 vector<double *> arrays_coord;
176 CHKERR read_util->get_node_coords(3, num_nodes, 0, start_v, arrays_coord);
177 Range verts(start_v, start_v + num_nodes - 1);
178 for (auto dd : {0, 1, 2}) {
179 std::copy(vert_coords[dd].begin(), vert_coords[dd].end(),
180 arrays_coord[dd]);
181 }
182 CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), verts,
183 &*parent_edge.begin());
184 CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
185 verts, bit, verb);
186 }
188}
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:345
@ 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 898 of file MeshRefinement.cpp.

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

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

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

192 {
193 Interface &m_field = cOre;
194 moab::Interface &moab = m_field.get_moab();
196 Range tets;
197 CHKERR moab.get_entities_by_type(meshset, MBTET, tets, false);
198 CHKERR refineTets(tets, bit, verb, debug);
200}
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 237 of file MeshRefinement.cpp.

239 {
240
242
243 auto set_edge_bits = [](moab::Interface &moab,
245 &ref_parent_ents_view,
246 EntityHandle tet, BitRefEdges &parent_edges_bit,
247 EntityHandle *edge_new_nodes, int *split_edges) {
249
250 for (int ee = 0; ee != 6; ++ee) {
251 EntityHandle edge;
252 CHKERR moab.side_element(tet, 1, ee, edge);
253 auto eit = ref_parent_ents_view.find(edge);
254 if (eit != ref_parent_ents_view.end()) {
255 edge_new_nodes[ee] = (*eit)->getEnt();
256 split_edges[parent_edges_bit.count()] = ee;
257 parent_edges_bit.set(ee, 1);
258 }
259 }
260
262 };
263
264 CHKERR refineTets(_tets, bit, set_edge_bits, verb, debug);
265
267}
#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 318 of file MeshRefinement.cpp.

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

308 {
309 Interface &m_field = cOre;
310 moab::Interface &moab = m_field.get_moab();
312 Range tets;
313 CHKERR moab.get_entities_by_type(meshset, MBTET, tets, false);
316}
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 269 of file MeshRefinement.cpp.

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

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

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

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

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

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

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

998 {
999 Interface &m_field = cOre;
1000 moab::Interface &moab = m_field.get_moab();
1002 Range tris;
1003 CHKERR moab.get_entities_by_type(meshset, MBTRI, tris, false);
1004 CHKERR refineTrisHangingNodes(tris, bit, verb, debug);
1006}
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 958 of file MeshRefinement.cpp.

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

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: