v0.13.1
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 TET 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 198 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:1086

◆ ~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:1955
virtual moab::Interface & get_moab()=0
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
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
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:494
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_IMPOSIBLE_CASE, "Imposible 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_IMPOSIBLE_CASE, "Imposible 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_IMPOSIBLE_CASE, "Imposible 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_IMPOSIBLE_CASE, "Imposible 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 int ref_type[2];
561 ref_type[0] = parent_edges_bit.count();
562 ref_type[1] = sub_type;
563
564 CHKERR moab.tag_clear_data(cOre.get_th_RefType(), &ref_tets[0],
565 nb_new_tets, ref_type);
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 // imposible 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
static Index< 'I', 3 > I
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:35
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
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:1504
Tag get_th_RefType() const
Definition: Core.hpp:200

◆ 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 TET 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 TET 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
Rangeof tets to refine
BitRefLevelbitLevel
BitRefLevelbitLevel
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_IMPOSIBLE_CASE, "Imposible 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 int ref_type[2];
1228 ref_type[0] = parent_edges_bit.count();
1229 ref_type[1] = 0;
1230
1231 // Set tags for reference dtriangles
1232 CHKERR moab.tag_clear_data(cOre.get_th_RefType(), &ref_tris[0],
1233 nb_new_tris, ref_type);
1234 CHKERR moab.tag_clear_data(cOre.get_th_RefBitEdge(), &ref_tris[0],
1235 nb_new_tris, &parent_edges_bit);
1236 CHKERR moab.tag_clear_data(cOre.get_th_RefParentHandle(), &ref_tris[0],
1237 nb_new_tris, &tit);
1238
1239 // Hash map of nodes (RefEntity) by edges (EntityHandle)
1240 std::map<EntityHandle /*edge*/, EntityHandle /*node*/>
1241 map_ref_nodes_by_edges;
1242
1243 // find parents for new edges and faces
1244 // get tet edges and faces
1245 Range tit_edges; // parent edges
1246 CHKERR moab.get_adjacencies(&tit, 1, 1, false, tit_edges);
1247 for (auto edge : tit_edges) {
1248 auto eit = ref_parent_ents_view.find(edge);
1249 if (eit != ref_parent_ents_view.end()) {
1250 map_ref_nodes_by_edges[(*eit)->getParentEnt()] = eit->get()->getEnt();
1251 }
1252 }
1253
1254 Range edges_nodes[3];
1255 // for edges - add ref nodes
1256 // edges_nodes[ee] - contains all nodes on edge ee including mid nodes if
1257 // exist
1258 Range::iterator eit = tit_edges.begin();
1259 for (int ee = 0; eit != tit_edges.end(); eit++, ee++) {
1260 CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee], true);
1261 auto map_miit = map_ref_nodes_by_edges.find(*eit);
1262 if (map_miit != map_ref_nodes_by_edges.end()) {
1263 edges_nodes[ee].insert(map_miit->second);
1264 }
1265 }
1266
1267 // add ref nodes to tri
1268 // tet_nodes contains all nodes on tet including mid edge nodes
1269 Range tri_nodes;
1270 CHKERR moab.get_connectivity(&tit, 1, tri_nodes, true);
1271 for (auto map_miit = map_ref_nodes_by_edges.begin();
1272 map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
1273 tri_nodes.insert(map_miit->second);
1274 }
1275
1276 Range ref_edges; // edges on all new tets
1277 // Get all all edges of refined tets
1278 CHKERR moab.get_adjacencies(ref_tris.data(), nb_new_tris, 1, false,
1279 ref_edges, moab::Interface::UNION);
1280
1281 // Check for all ref edge and set parents
1282 for (auto reit = ref_edges.begin(); reit != ref_edges.end(); ++reit) {
1283 Range ref_edges_nodes;
1284 CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes, true);
1285 // Check if ref edge is an coarse edge (loop over coarse tet edges)
1286 int ee = 0;
1287 for (; ee < 3; ee++) {
1288 // Two nodes are common (node[0],node[1],ref_node (if exist))
1289 // this tests if given edge is contained by edge of refined
1290 // triangle
1291 if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
1292 EntityHandle edge = tit_edges[ee];
1293 CHKERR set_parent(*reit, edge, refined_ents_ptr, cOre);
1294 break;
1295 }
1296 }
1297 if (ee < 3)
1298 continue; // this refined edge is contained by edge of triangel
1299
1300 // check if ref edge is in coarse tetrahedral (i.e. that is internal
1301 // edge of refined tetrahedral)
1302 if (intersect(tri_nodes, ref_edges_nodes).size() == 2) {
1303 CHKERR set_parent(*reit, tit, refined_ents_ptr, cOre);
1304 continue;
1305 }
1306
1307 // Refined edge is not child of any edge, face or tetrahedral, this is
1308 // imposible edge
1309 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1310 "impossible refined edge");
1311 }
1312 }
1313 }
1314
1315 CHKERR set_parent(refined_ents_ptr);
1316 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
1317 ents_to_set_bit, bit, true, verb, &ents);
1318
1320}
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 temprary meshset.
Definition: Templates.hpp:1465

◆ 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 tets to refine
BitRefLevelbitLevel
BitRefLevelbitLevel
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 tets to refine
BitRefLevelbitLevel
BitRefLevelbitLevel
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: