v0.14.0
Loading...
Searching...
No Matches
Classes | Public Member Functions | Public Attributes | Private Types | Private Member Functions | List of all members
MoFEM::MeshRefinement Struct Reference

Mesh refinement interface. More...

#include <src/interfaces/MeshRefinement.hpp>

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

Classes

struct  SetParent
 

Public Member Functions

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

Public Attributes

MoFEM::CorecOre
 

Private Types

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

Private Member Functions

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

Additional Inherited Members

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

Detailed Description

Mesh refinement interface.

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

Bug:

Not working on partitioned meshes

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

If outsourced, class member functions should follow name convention

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

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

Definition at line 26 of file MeshRefinement.hpp.

Member Typedef Documentation

◆ SetEdgeBitsFun

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

Functions setting edges for refinemnt on enetity level.

Definition at line 195 of file MeshRefinement.hpp.

Constructor & Destructor Documentation

◆ MeshRefinement()

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

Definition at line 39 of file MeshRefinement.cpp.

40 : cOre(const_cast<Core &>(core)) {}
CoreTmp< 0 > Core
Definition: Core.hpp:1094

◆ ~MeshRefinement()

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

Member Function Documentation

◆ addVerticesInTheMiddleOfEdges() [1/2]

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

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

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

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

Definition at line 42 of file MeshRefinement.cpp.

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

◆ addVerticesInTheMiddleOfEdges() [2/2]

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

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

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

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

Definition at line 96 of file MeshRefinement.cpp.

97 {
98 Interface &m_field = cOre;
99 moab::Interface &moab = m_field.get_moab();
100 auto refined_ents_ptr = m_field.get_ref_ents();
102 auto miit =
103 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
104 .lower_bound(boost::make_tuple(MBVERTEX, MBEDGE));
105 auto hi_miit =
106 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
107 .upper_bound(boost::make_tuple(MBVERTEX, MBEDGE));
109 ref_parent_ents_view.insert(miit, hi_miit);
110
111 Range edges = ents.subset_by_type(MBEDGE);
112 if (verb >= VERBOSE) {
113 MOFEM_TAG_AND_LOG("WORLD", Sev::verbose, "MeshRefinemnt")
114 << "ref level " << bit << " nb. edges to refine " << edges.size();
115 }
116
117 std::array<std::vector<double>, 3> vert_coords;
118 for (auto &vc : vert_coords)
119 vc.reserve(edges.size());
120
121 std::vector<EntityHandle> parent_edge;
122 parent_edge.reserve(edges.size());
123
124 std::array<double, 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:362
@ VERBOSE
Definition: definitions.h:209
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_non_unique< const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< boost::shared_ptr< RefEntity >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_ordered_parent_entity
Tag get_th_RefParentHandle() const
Definition: Core.hpp:197

◆ query_interface()

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

Implements MoFEM::UnknownInterface.

Definition at line 33 of file MeshRefinement.cpp.

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

◆ refineMeshset()

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

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

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

Definition at line 897 of file MeshRefinement.cpp.

899 {
900 Interface &m_field = cOre;
901 auto refined_ents_ptr = m_field.get_ref_ents();
903 auto miit = refined_ents_ptr->find(meshset);
904 if (miit == refined_ents_ptr->end()) {
905 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
906 "this meshset is not in ref database");
907 }
908 CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
909 meshset, bit, meshset, MBMAXTYPE, recursive, verb);
910 *(const_cast<RefEntity *>(miit->get())->getBitRefLevelPtr()) |= bit;
912}
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 715 of file MeshRefinement.cpp.

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

◆ refineTets() [1/3]

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

refine TET in the meshset

Parameters
EntityHandlemeshset
BitRefLevelbitLevel
verbverbosity level

Definition at line 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
553 std::array<EntityHandle, 8> ref_tets;
554 for (int tt = 0; tt != nb_new_tets; ++tt, ++start_e)
555 ref_tets[tt] = start_e;
556
557 if (nb_new_tets) {
558
559 CHKERR moab.tag_clear_data(cOre.get_th_RefBitEdge(), &ref_tets[0],
560 nb_new_tets, &parent_edges_bit);
561 CHKERR moab.tag_clear_data(cOre.get_th_RefParentHandle(), &ref_tets[0],
562 nb_new_tets, &tit);
563
564 // hash map of nodes (RefEntity) by edges (EntityHandle)
565 std::map<EntityHandle /*edge*/, EntityHandle /*node*/>
566 map_ref_nodes_by_edges;
567
568 Range tit_edges;
569 CHKERR moab.get_adjacencies(&tit, 1, 1, false, tit_edges);
570 for (auto edge : tit_edges) {
571 auto eit = ref_parent_ents_view.find(edge);
572 if (eit != ref_parent_ents_view.end()) {
573 map_ref_nodes_by_edges[(*eit)->getParentEnt()] = eit->get()->getEnt();
574 }
575 }
576
577 // find parents for new edges and faces
578 // get tet edges and faces
579 Range tit_faces;
580 CHKERR moab.get_adjacencies(&tit, 1, 2, false, tit_faces);
581 Range edges_nodes[6], faces_nodes[4];
582 // for edges - add ref nodes
583 // edges_nodes[ee] - contains all nodes on edge ee including mid nodes if
584 // exist
585 Range::iterator eit = tit_edges.begin();
586 for (int ee = 0; eit != tit_edges.end(); eit++, ee++) {
587 CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee], true);
588 auto map_miit = map_ref_nodes_by_edges.find(*eit);
589 if (map_miit != map_ref_nodes_by_edges.end()) {
590 edges_nodes[ee].insert(map_miit->second);
591 }
592 }
593 // for faces - add ref nodes
594 // faces_nodes[ff] - contains all nodes on face ff including mid nodes if
595 // exist
596 Range::iterator fit = tit_faces.begin();
597 for (int ff = 0; fit != tit_faces.end(); fit++, ff++) {
598 CHKERR moab.get_connectivity(&*fit, 1, faces_nodes[ff], true);
599 // Get edges on face and loop over those edges to add mid-nodes to range
600 Range fit_edges;
601 CHKERR moab.get_adjacencies(&*fit, 1, 1, false, fit_edges);
602 for (Range::iterator eit2 = fit_edges.begin(); eit2 != fit_edges.end();
603 eit2++) {
604 auto map_miit = map_ref_nodes_by_edges.find(*eit2);
605 if (map_miit != map_ref_nodes_by_edges.end()) {
606 faces_nodes[ff].insert(map_miit->second);
607 }
608 }
609 }
610 // add ref nodes to tet
611 // tet_nodes contains all nodes on tet including mid edge nodes
612 Range tet_nodes;
613 CHKERR moab.get_connectivity(&tit, 1, tet_nodes, true);
614 for (std::map<EntityHandle, EntityHandle>::iterator map_miit =
615 map_ref_nodes_by_edges.begin();
616 map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
617 tet_nodes.insert(map_miit->second);
618 }
619 Range ref_edges;
620 // Get all all edges of refined tets
621 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 1, false,
622 ref_edges, moab::Interface::UNION);
623 // Check for all ref edge and set parents
624 for (Range::iterator reit = ref_edges.begin(); reit != ref_edges.end();
625 reit++) {
626 Range ref_edges_nodes;
627 CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes, true);
628
629 // Check if ref edge is an coarse edge (loop over coarse tet edges)
630 int ee = 0;
631 for (; ee < 6; ee++) {
632 // Two nodes are common (node[0],node[1],ref_node (if exist))
633 // this tests if given edge is contained by edge of refined
634 // tetrahedral
635 if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
636 EntityHandle edge = tit_edges[ee];
637 CHKERR set_parent(*reit, edge, refined_ents_ptr, cOre);
638 break;
639 }
640 }
641 if (ee < 6)
642 continue; // this refined edge is contained by edge of tetrahedral
643 // check if ref edge is in coarse face
644 int ff = 0;
645 for (; ff < 4; ff++) {
646 // two nodes are common (node[0],node[1],ref_node (if exist))
647 // this tests if given face is contained by face of tetrahedral
648 if (intersect(faces_nodes[ff], ref_edges_nodes).size() == 2) {
649 EntityHandle face = tit_faces[ff];
650 CHKERR set_parent(*reit, face, refined_ents_ptr, cOre);
651 break;
652 }
653 }
654 if (ff < 4)
655 continue; // this refined edge is contained by face of tetrahedral
656
657 // check if ref edge is in coarse tetrahedral (i.e. that is internal
658 // edge of refined tetrahedral)
659 if (intersect(tet_nodes, ref_edges_nodes).size() == 2) {
660 CHKERR set_parent(*reit, tit, refined_ents_ptr, cOre);
661 continue;
662 }
663
664 // Refined edge is not child of any edge, face or tetrahedral, this is
665 // Impossible edge
666 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
667 "impossible refined edge");
668 }
669
670 Range ref_faces;
671 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 2, false,
672 ref_faces, moab::Interface::UNION);
673 Tag th_interface_side;
674 const int def_side[] = {0};
675 CHKERR moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
676 th_interface_side,
677 MB_TAG_CREAT | MB_TAG_SPARSE, def_side);
678 // Check for all ref faces
679 for (auto rfit : ref_faces) {
680 Range ref_faces_nodes;
681 CHKERR moab.get_connectivity(&rfit, 1, ref_faces_nodes, true);
682 // Check if ref face is in coarse face
683 int ff = 0;
684 for (; ff < 4; ff++) {
685 // Check if refined triangle is contained by face of tetrahedral
686 if (intersect(faces_nodes[ff], ref_faces_nodes).size() == 3) {
687 EntityHandle face = tit_faces[ff];
688 CHKERR set_parent(rfit, face, refined_ents_ptr, cOre);
689 int side = 0;
690 // Set face side if it is on interface
691 CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
692 CHKERR moab.tag_set_data(th_interface_side, &rfit, 1, &side);
693 break;
694 }
695 }
696 if (ff < 4)
697 continue; // this face is contained by one of tetrahedrons
698 // check if ref face is in coarse tetrahedral
699 // this is ref face which is contained by tetrahedral volume
700 if (intersect(tet_nodes, ref_faces_nodes).size() == 3) {
701 CHKERR set_parent(rfit, tit, refined_ents_ptr, cOre);
702 continue;
703 }
704 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
705 "impossible refined face");
706 }
707 }
708 }
709
710 CHKERR set_parent(refined_ents_ptr);
711 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
712 ents_to_set_bit, bit, true, verb, &ents);
714}
static Index< 'p', 3 > p
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
int tet_type_3(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
void tet_type_6(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
int tet_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
int tet_type_5(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
int tet_type_4(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
void tet_type_1(const EntityHandle *conn, const int split_edge, const EntityHandle edge_new_node, EntityHandle *new_tets_conn)
auto rangeInserter(const I f, const I s, boost::function< bool(I it)> tester, boost::function< MoFEMErrorCode(I f, I s)> inserter)
Insert ranges.
Definition: Templates.hpp:1826
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 914 of file MeshRefinement.cpp.

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

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

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

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

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

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

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

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: