v0.8.23
Public Member Functions | Public Attributes | 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]

Public Member Functions

MoFEMErrorCode query_interface (const MOFEMuuid &uuid, UnknownInterface **iface) const
 
 MeshRefinement (const MoFEM::Core &core)
 
 ~MeshRefinement ()
 destructor More...
 
virtual MoFEMErrorCode add_vertices_in_the_middle_of_edges (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...
 
DEPRECATED MoFEMErrorCode add_verices_in_the_middel_of_edges (const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
 
virtual MoFEMErrorCode add_vertices_in_the_middle_of_edges (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...
 
DEPRECATED MoFEMErrorCode add_verices_in_the_middel_of_edges (const Range &edges, const BitRefLevel &bit, int verb=QUIET, EntityHandle start_v=0)
 
virtual MoFEMErrorCode refine_TET (const EntityHandle meshset, const BitRefLevel &bit, const bool respect_interface=false, int verb=QUIET, Range *ref_edges=NULL, const bool debug=false)
 refine TET in the meshset More...
 
virtual MoFEMErrorCode refine_TET (const Range &tets, const BitRefLevel &bit, const bool respect_interface=false, int verb=QUIET, Range *ref_edges=NULL, const bool debug=false)
 refine TET in the meshset More...
 
virtual MoFEMErrorCode refine_PRISM (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET)
 refine PRISM in the meshset More...
 
virtual MoFEMErrorCode refine_MESHSET (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...
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (const MOFEMuuid &uuid, bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE , bool VERIFY = false>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, IFACE *&iface) const
 Get interface by uuid and return reference to pointer of 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 ()
 
virtual MoFEMErrorCode getLibVersion (Version &version) const
 Get library version. More...
 
virtual const MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version) const
 Get database major version. More...
 
virtual MoFEMErrorCode getInterfaceVersion (Version &version) const
 Get database major version. More...
 
template<>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, UnknownInterface *&iface) const
 

Public Attributes

MoFEM::CorecOre
 

Additional Inherited Members

- Protected Member Functions inherited from MoFEM::UnknownInterface
boost::typeindex::type_index getClassIdx (const MOFEMuuid &uid) const
 Get type name for interface Id. More...
 
MOFEMuuid getUId (const boost::typeindex::type_index &class_idx) const
 Get interface Id for class name. 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_genrate_base.cpp.

Definition at line 41 of file MeshRefinement.hpp.

Constructor & Destructor Documentation

◆ MeshRefinement()

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

Definition at line 34 of file MeshRefinement.cpp.

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

◆ ~MeshRefinement()

MoFEM::MeshRefinement::~MeshRefinement ( )

destructor

Definition at line 50 of file MeshRefinement.hpp.

50 {}

Member Function Documentation

◆ add_verices_in_the_middel_of_edges() [1/2]

DEPRECATED MoFEMErrorCode MoFEM::MeshRefinement::add_verices_in_the_middel_of_edges ( const EntityHandle  meshset,
const BitRefLevel bit,
const bool  recursive = false,
int  verb = QUIET,
EntityHandle  start_v = 0 
)
Deprecated:
Use function with correct spelling

Definition at line 72 of file MeshRefinement.hpp.

75  {
76  return add_vertices_in_the_middle_of_edges(meshset, bit, recursive, verb,
77  start_v);
78  }
virtual MoFEMErrorCode add_vertices_in_the_middle_of_edges(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

◆ add_verices_in_the_middel_of_edges() [2/2]

DEPRECATED MoFEMErrorCode MoFEM::MeshRefinement::add_verices_in_the_middel_of_edges ( const Range &  edges,
const BitRefLevel bit,
int  verb = QUIET,
EntityHandle  start_v = 0 
)
Deprecated:
Use function with correct spelling

Definition at line 101 of file MeshRefinement.hpp.

103  {
104  return add_vertices_in_the_middle_of_edges(edges, bit, verb, start_v);
105  }
virtual MoFEMErrorCode add_vertices_in_the_middle_of_edges(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

◆ add_vertices_in_the_middle_of_edges() [1/2]

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

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

39  {
40  Interface &m_field = cOre;
41  moab::Interface &moab = m_field.get_moab();
43  Range edges;
44  CHKERR moab.get_entities_by_type(meshset, MBEDGE, edges, recursive);
45  if (edges.empty()) {
46  Range tets;
47  CHKERR moab.get_entities_by_type(meshset, MBTET, tets, recursive);
48  CHKERR moab.get_adjacencies(tets, 1, true, edges, moab::Interface::UNION);
49  if (tets.empty()) {
50  Range prisms;
51  CHKERR moab.get_entities_by_type(meshset, MBPRISM, prisms, recursive);
52  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
53  const EntityHandle *conn;
54  int num_nodes;
55  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
56  assert(num_nodes == 6);
57  //
58  Range edge;
59  CHKERR moab.get_adjacencies(&conn[0], 2, 1, true, edge);
60  assert(edge.size() == 1);
61  edges.insert(edge[0]);
62  edge.clear();
63  CHKERR moab.get_adjacencies(&conn[1], 2, 1, true, edge);
64  assert(edge.size() == 1);
65  edges.insert(edge[0]);
66  EntityHandle conn_edge2[] = {conn[2], conn[0]};
67  edge.clear();
68  CHKERR moab.get_adjacencies(conn_edge2, 2, 1, true, edge);
69  assert(edge.size() == 1);
70  edges.insert(edge[0]);
71  //
72  edge.clear();
73  CHKERR moab.get_adjacencies(&conn[3], 2, 1, true, edge);
74  assert(edge.size() == 1);
75  edges.insert(edge[0]);
76  edge.clear();
77  CHKERR moab.get_adjacencies(&conn[4], 2, 1, true, edge);
78  assert(edge.size() == 1);
79  edges.insert(edge[0]);
80  EntityHandle conn_edge8[] = {conn[5], conn[3]};
81  edge.clear();
82  CHKERR moab.get_adjacencies(conn_edge8, 2, 1, true, edge);
83  assert(edge.size() == 1);
84  edges.insert(edge[0]);
85  }
86  }
87  }
88  CHKERR add_vertices_in_the_middle_of_edges(edges, bit, verb, start_v);
90 }
virtual MoFEMErrorCode add_vertices_in_the_middle_of_edges(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
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ add_vertices_in_the_middle_of_edges() [2/2]

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

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

92  {
93  Interface &m_field = cOre;
94  moab::Interface &moab = m_field.get_moab();
95  const RefEntity_multiIndex *refined_ents_ptr;
97  CHKERR m_field.get_ref_ents(&refined_ents_ptr);
98  auto miit =
99  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
100  .lower_bound(boost::make_tuple(MBVERTEX, MBEDGE));
101  auto hi_miit =
102  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
103  .upper_bound(boost::make_tuple(MBVERTEX, MBEDGE));
105  ref_parent_ents_view.insert(miit, hi_miit);
106 
107  Range edges = ents.subset_by_type(MBEDGE);
108  if (verb >= VERBOSE) {
109  std::ostringstream ss;
110  ss << "ref level " << bit << " nb. edges to refine " << edges.size()
111  << std::endl;
112  PetscPrintf(m_field.get_comm(), ss.str().c_str());
113  }
114 
115  std::array<std::vector<double>, 3> vert_coords;
116  for (auto &vc : vert_coords)
117  vc.reserve(edges.size());
118 
119  std::vector<EntityHandle> parent_edge;
120  parent_edge.reserve(edges.size());
121 
122  std::array<double, 6> coords;
124  &coords[0], &coords[1], &coords[2]};
126  &coords[3], &coords[4], &coords[5]};
128 
129  Range add_bit;
130  for (auto p_eit = edges.pair_begin(); p_eit != edges.pair_end(); ++p_eit) {
131  auto miit_view = ref_parent_ents_view.lower_bound(p_eit->first);
132 
133  Range edge_having_parent_vertex;
134  if (miit_view != ref_parent_ents_view.end()) {
135  for (auto hi_miit_view = ref_parent_ents_view.upper_bound(p_eit->second);
136  miit_view != hi_miit_view; ++miit_view) {
137  edge_having_parent_vertex.insert(edge_having_parent_vertex.end(),
138  miit_view->get()->getParentEnt());
139  add_bit.insert(add_bit.end(), miit_view->get()->getRefEnt());
140  }
141  }
142 
143  Range add_vertex_edges =
144  subtract(Range(p_eit->first, p_eit->second), edge_having_parent_vertex);
145 
146  for (auto e : add_vertex_edges)
147  parent_edge.emplace_back(e);
148 
149  for (auto e : add_vertex_edges) {
150 
151  const EntityHandle *conn;
152  int num_nodes;
153  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
154  if (PetscUnlikely(num_nodes != 2)) {
155  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
156  "edge should have 2 edges");
157  }
158  CHKERR moab.get_coords(conn, num_nodes, coords.data());
159  t_0(i) += t_1(i);
160  t_0(i) *= 0.5;
161 
162  for (auto j : {0, 1, 2})
163  vert_coords[j].emplace_back(t_0(j));
164  }
165  }
166 
167  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit, bit);
168 
169  if (!vert_coords[0].empty()) {
170  ReadUtilIface *read_util;
171  CHKERR moab.query_interface(read_util);
172  int num_nodes = vert_coords[0].size();
173  vector<double *> arrays_coord;
174  CHKERR read_util->get_node_coords(3, num_nodes, 0, start_v, arrays_coord);
175  Range verts(start_v, start_v + num_nodes - 1);
176  for (auto dd : {0, 1, 2}) {
177  std::copy(vert_coords[dd].begin(), vert_coords[dd].end(),
178  arrays_coord[dd]);
179  }
180  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), verts,
181  &*parent_edge.begin());
182  CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
183  verts, bit, verb);
184  }
186 }
Tag get_th_RefParentHandle() const
Definition: Core.hpp:150
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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::getRefEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_ordered_parent_entity
#define CHKERR
Inline error check.
Definition: definitions.h:596
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, 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::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ query_interface()

MoFEMErrorCode MoFEM::MeshRefinement::query_interface ( const MOFEMuuid uuid,
UnknownInterface **  iface 
) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 22 of file MeshRefinement.cpp.

23  {
25  *iface = NULL;
26  if (uuid == IDD_MOFEMMeshRefine) {
27  *iface = const_cast<MeshRefinement *>(this);
29  }
30  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
32 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static const MOFEMuuid IDD_MOFEMMeshRefine

◆ refine_MESHSET()

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

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

959  {
960  Interface &m_field = cOre;
961  const RefEntity_multiIndex *refined_ents_ptr;
963  CHKERR m_field.get_ref_ents(&refined_ents_ptr);
964  typedef const RefEntity_multiIndex::index<Ent_mi_tag>::type RefEntsByEnt;
965  RefEntsByEnt::iterator miit = refined_ents_ptr->find(meshset);
966  if (miit == refined_ents_ptr->end()) {
967  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
968  "this meshset is not in ref database");
969  }
970  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
971  meshset, bit, meshset, MBEDGE, recursive, verb);
972  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
973  meshset, bit, meshset, MBTRI, recursive, verb);
974  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
975  meshset, bit, meshset, MBTET, recursive, verb);
976  *(const_cast<RefEntity *>(miit->get())->getBitRefLevelPtr()) |= bit;
978 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, 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::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ refine_PRISM()

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

refine PRISM in the meshset

Parameters
EntityHandlemeshset
BitRefLevelbitLevel

Definition at line 764 of file MeshRefinement.cpp.

765  {
766 
767  Interface &m_field = cOre;
768  moab::Interface &moab = m_field.get_moab();
769  const RefEntity_multiIndex *refined_ents_ptr;
770  const RefElement_multiIndex *refined_finite_elements_ptr;
771 
772  // FIXME: refinement is based on entity handlers, should work on global ids of
773  // nodes, this will allow parallelise algorithm in the future
774 
776 
777  CHKERR m_field.get_ref_ents(&refined_ents_ptr);
778  CHKERR m_field.get_ref_finite_elements(&refined_finite_elements_ptr);
779 
780  typedef const RefEntity_multiIndex::index<Ent_mi_tag>::type RefEntsByEnt;
781 
782  typedef const RefElement_multiIndex_parents_view::index<
783  Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>::type
784  RefEntByParentAndRefEdges;
785  RefElement_multiIndex_parents_view ref_ele_parent_view;
786  ref_ele_parent_view.insert(
787  refined_finite_elements_ptr->get<EntType_mi_tag>().lower_bound(MBPRISM),
788  refined_finite_elements_ptr->get<EntType_mi_tag>().upper_bound(MBPRISM));
789  RefEntByParentAndRefEdges &ref_ele_by_parent_and_ref_edges =
790  ref_ele_parent_view
791  .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
792  // find all vertices which parent is edge
793  typedef const RefEntity_multiIndex::index<
794  Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
795  RefEntsByComposite &ref_ents_by_comp =
796  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
798  ref_parent_ents_view.insert(
799  ref_ents_by_comp.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE)),
800  ref_ents_by_comp.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE)));
801  Range prisms;
802  CHKERR moab.get_entities_by_type(meshset, MBPRISM, prisms, false);
803  Range::iterator pit = prisms.begin();
804  for (; pit != prisms.end(); pit++) {
805  RefEntsByEnt::iterator miit_prism =
806  refined_ents_ptr->get<Ent_mi_tag>().find(*pit);
807  if (miit_prism == refined_ents_ptr->end()) {
808  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
809  "this prism is not in ref database");
810  }
811  if (verb >= NOISY) {
812  std::ostringstream ss;
813  ss << "ref prism " << **miit_prism << std::endl;
814  PetscPrintf(m_field.get_comm(), ss.str().c_str());
815  }
816  // prism connectivity
817  int num_nodes;
818  const EntityHandle *conn;
819  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
820  assert(num_nodes == 6);
821  // edges connectivity
822  EntityHandle edges[6];
823  for (int ee = 0; ee < 3; ee++) {
824  CHKERR moab.side_element(*pit, 1, ee, edges[ee]);
825  }
826  for (int ee = 6; ee < 9; ee++) {
827  CHKERR moab.side_element(*pit, 1, ee, edges[ee - 3]);
828  }
829  // detect split edges
830  BitRefEdges split_edges(0);
831  EntityHandle edge_nodes[6];
832  std::fill(&edge_nodes[0], &edge_nodes[6], no_handle);
833  for (int ee = 0; ee < 6; ee++) {
834  RefEntity_multiIndex_view_by_hashed_parent_entity::iterator miit_view =
835  ref_parent_ents_view.find(edges[ee]);
836  if (miit_view != ref_parent_ents_view.end()) {
837  if (((*miit_view)->getBitRefLevel() & bit).any()) {
838  edge_nodes[ee] = (*miit_view)->getRefEnt();
839  split_edges.set(ee);
840  }
841  }
842  }
843  if (split_edges.count() == 0) {
844  *(const_cast<RefEntity *>(miit_prism->get())->getBitRefLevelPtr()) |= bit;
845  if (verb >= VERY_NOISY)
846  PetscPrintf(m_field.get_comm(), "no refinement");
847  continue;
848  }
849  // check consistency
850  if (verb >= NOISY) {
851  std::ostringstream ss;
852  ss << "prism split edges " << split_edges << " count "
853  << split_edges.count() << std::endl;
854  PetscPrintf(m_field.get_comm(), ss.str().c_str());
855  }
856  // prism ref
857  EntityHandle new_prism_conn[4 * 6];
858  std::fill(&new_prism_conn[0], &new_prism_conn[4 * 6], no_handle);
859  int nb_new_prisms = 0;
860  switch (split_edges.count()) {
861  case 0:
862  break;
863  case 2:
864  CHKERR prism_type_1(conn, split_edges, edge_nodes, new_prism_conn);
865  nb_new_prisms = 2;
866  break;
867  case 4:
868  CHKERR prism_type_2(conn, split_edges, edge_nodes, new_prism_conn);
869  nb_new_prisms = 3;
870  break;
871  case 6:
872  CHKERR prism_type_3(conn, split_edges, edge_nodes, new_prism_conn);
873  nb_new_prisms = 4;
874  break;
875  default:
876  std::ostringstream ss;
877  ss << split_edges << " : [ " << conn[0] << " " << conn[1] << " "
878  << conn[2] << " " << conn[3] << " " << conn[4] << " " << conn[5]
879  << " ]";
880  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, ss.str().c_str());
881  }
882  // find that prism
883  std::bitset<4> ref_prism_bit(0);
884  RefEntByParentAndRefEdges::iterator it_by_ref_edges =
885  ref_ele_by_parent_and_ref_edges.lower_bound(
886  boost::make_tuple(*pit, split_edges.to_ulong()));
887  RefEntByParentAndRefEdges::iterator hi_it_by_ref_edges =
888  ref_ele_by_parent_and_ref_edges.upper_bound(
889  boost::make_tuple(*pit, split_edges.to_ulong()));
890  RefEntByParentAndRefEdges::iterator it_by_ref_edges2 = it_by_ref_edges;
891  for (int pp = 0; it_by_ref_edges2 != hi_it_by_ref_edges;
892  it_by_ref_edges2++, pp++) {
893  // Add this tet to this ref
894  *(const_cast<RefElement *>(it_by_ref_edges2->get())
895  ->getBitRefLevelPtr()) |= bit;
896  ref_prism_bit.set(pp, 1);
897  if (verb > 2) {
898  std::ostringstream ss;
899  ss << "is refined " << *it_by_ref_edges2 << std::endl;
900  PetscPrintf(m_field.get_comm(), ss.str().c_str());
901  }
902  }
903  if (it_by_ref_edges != hi_it_by_ref_edges) {
904  if (ref_prism_bit.count() != (unsigned int)nb_new_prisms)
905  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
906  "data inconsistency");
907  } else {
908  EntityHandle ref_prisms[4];
909  // create prism
910  for (int pp = 0; pp < nb_new_prisms; pp++) {
911  if (verb > 3) {
912  std::ostringstream ss;
913  ss << "ref prism " << ref_prism_bit << std::endl;
914  PetscPrintf(m_field.get_comm(), ss.str().c_str());
915  }
916  if (!ref_prism_bit.test(pp)) {
917  CHKERR moab.create_element(MBPRISM, &new_prism_conn[6 * pp], 6,
918  ref_prisms[pp]);
919  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(),
920  &ref_prisms[pp], 1, &*pit);
921  CHKERR moab.tag_set_data(cOre.get_th_RefBitLevel(), &ref_prisms[pp],
922  1, &bit);
923  CHKERR moab.tag_set_data(cOre.get_th_RefBitEdge(), &ref_prisms[pp], 1,
924  &split_edges);
925  std::pair<RefEntity_multiIndex::iterator, bool> p_ent =
926  const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
927  ->insert(boost::shared_ptr<RefEntity>(new RefEntity(
928  m_field.get_basic_entity_data_ptr(), ref_prisms[pp])));
929  std::pair<RefElement_multiIndex::iterator, bool> p_fe;
930  try {
931  p_fe =
932  const_cast<RefElement_multiIndex *>(refined_finite_elements_ptr)
933  ->insert(boost::shared_ptr<RefElement>(
934  new RefElement_PRISM(*p_ent.first)));
935  } catch (MoFEMException const &e) {
936  SETERRQ(PETSC_COMM_SELF, e.errorCode, e.errorMessage);
937  }
938  ref_prism_bit.set(pp);
939  CHKERR cOre.addPrismToDatabase(ref_prisms[pp]);
940  if (verb > 2) {
941  std::ostringstream ss;
942  ss << "add prism: " << **p_fe.first << std::endl;
943  if (verb > 7) {
944  for (int nn = 0; nn < 6; nn++) {
945  ss << new_prism_conn[nn] << " ";
946  }
947  ss << std::endl;
948  }
949  PetscPrintf(m_field.get_comm(), ss.str().c_str());
950  }
951  }
952  }
953  }
954  }
956 }
Tag get_th_RefParentHandle() const
Definition: Core.hpp:150
Tag get_th_RefBitLevel() const
Definition: Core.hpp:151
multi_index_container< boost::shared_ptr< RefElement >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getRefEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityType, &RefElement::getEntType > > > > RefElement_multiIndex
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition: Common.hpp:27
MoFEMErrorCode prism_type_2(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
Tag get_th_RefBitEdge() const
Definition: Core.hpp:152
std::bitset< BITREFEDGES_SIZE > BitRefEdges
Definition: Types.hpp:45
MoFEMErrorCode prism_type_1(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
MoFEMErrorCode addPrismToDatabase(const EntityHandle prism, int verb=DEFAULT_VERBOSITY)
add prim element
Definition: Core.cpp:267
MoFEMErrorCode prism_type_3(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::getRefEnt > >, 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
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::getRefEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_hashed_parent_entity
multi-index view of RefEntity by parent entity
#define CHKERR
Inline error check.
Definition: definitions.h:596
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, 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::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ refine_TET() [1/2]

MoFEMErrorCode MoFEM::MeshRefinement::refine_TET ( const EntityHandle  meshset,
const BitRefLevel bit,
const bool  respect_interface = false,
int  verb = QUIET,
Range *  ref_edges = NULL,
const bool  debug = false 
)
virtual

refine TET in the meshset

Parameters
EntityHandlemeshset
BitRefLevelbitLevel
boolrespect_interface If TRUE, interface elements would be refined
verbverbosity level
Range* pointer to range in not null check consistency of refinement

Definition at line 188 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 refine_TET(tets, bit, respect_interface, verb, ref_edges_ptr, debug);
200 }
virtual MoFEMErrorCode refine_TET(const EntityHandle meshset, const BitRefLevel &bit, const bool respect_interface=false, int verb=QUIET, Range *ref_edges=NULL, const bool debug=false)
refine TET in the meshset
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ refine_TET() [2/2]

MoFEMErrorCode MoFEM::MeshRefinement::refine_TET ( const Range &  tets,
const BitRefLevel bit,
const bool  respect_interface = false,
int  verb = QUIET,
Range *  ref_edges = NULL,
const bool  debug = false 
)
virtual

refine TET in the meshset

Parameters
Rangeof tets to refine
BitRefLevelbitLevel
BitRefLevelbitLevel
boolrespect_interface If TRUE, interface elements would be refined
verbverbosity level
Range* pointer to range in not null check consistency of refinement

Definition at line 202 of file MeshRefinement.cpp.

206  {
207 
208  Interface &m_field = cOre;
209  moab::Interface &moab = m_field.get_moab();
210  const RefEntity_multiIndex *refined_ents_ptr;
211  const RefElement_multiIndex *refined_finite_elements_ptr;
212  ReadUtilIface *read_util;
214 
215  CHKERR m_field.get_moab().query_interface(read_util);
216 
217  // Check if refinement is correct
218  struct Check {
219  map<EntityHandle, EntityHandle> entParentMap;
220  MoFEMErrorCode operator()(Range *ref_edges,
221  const RefEntity_multiIndex *ref_ents_ptr,
222  MoFEM::Core &core) {
224  if (!ref_edges)
226  for (Range::iterator eit = ref_edges->begin(); eit != ref_edges->end();
227  ++eit) {
228  RefEntity_multiIndex::index<
229  Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vit =
230  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
231  boost::make_tuple(MBVERTEX, *eit));
232  if (vit ==
233  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
234  RefEntity_multiIndex::iterator e_eit = ref_ents_ptr->find(*eit);
235  if (e_eit == ref_ents_ptr->end()) {
236  SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
237  "Edge not found %ld", *eit);
238  }
239  cerr << "Parent edge" << endl << **e_eit << endl;
240  if (entParentMap.find(*eit) != entParentMap.end()) {
241  RefEntity_multiIndex::iterator v_eit =
242  ref_ents_ptr->find(entParentMap[*eit]);
243  if (v_eit == ref_ents_ptr->end()) {
244  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
245  "Vertex not found");
246  }
247  cerr << "Vertex " << **v_eit << endl;
248  }
249  RefEntity_multiIndex::index<Ent_Ent_mi_tag>::type::iterator ee_it,
250  ee_hi_it;
251  ee_it = ref_ents_ptr->get<Ent_Ent_mi_tag>().lower_bound(*eit);
252  ee_hi_it = ref_ents_ptr->get<Ent_Ent_mi_tag>().upper_bound(*eit);
253  for (; ee_it != ee_hi_it; ++ee_it) {
254  cerr << "Ent having edge parent by parent " << **ee_it << endl;
255  }
256  RefEntity_multiIndex tmp_index;
257  tmp_index.insert(ref_ents_ptr->begin(), ref_ents_ptr->end());
258  RefEntity_multiIndex::index<
259  Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vvit =
260  tmp_index.get<Composite_ParentEnt_And_EntType_mi_tag>().find(
261  boost::make_tuple(MBVERTEX, *eit));
262  if (vvit !=
263  tmp_index.get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
264  cerr << "Tmp idx Vertex " << **vvit << endl;
265  }
266  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
267  "No vertex on trim edges, that make no sense");
268  } else {
269  entParentMap[vit->get()->getParentEnt()] = vit->get()->getRefEnt();
270  }
271  }
273  }
274  };
275 
276  struct SetParent {
277  map<EntityHandle, EntityHandle> parentsToChange;
278  MoFEMErrorCode operator()(const EntityHandle ent, const EntityHandle parent,
279  const RefEntity_multiIndex *ref_ents_ptr,
280  MoFEM::Core &cOre) {
281  MoFEM::Interface &m_field = cOre;
283  RefEntity_multiIndex::iterator it = ref_ents_ptr->find(ent);
284  if (it != ref_ents_ptr->end()) {
285  if (it->get()->getParentEnt() != parent && ent != parent) {
286  parentsToChange[ent] = parent;
287  }
288  } else {
289  if (ent != parent) {
290  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
291  &ent, 1, &parent);
292  }
293  }
295  }
296  MoFEMErrorCode operator()(const RefEntity_multiIndex *ref_ents_ptr) {
298  for (map<EntityHandle, EntityHandle>::iterator mit =
299  parentsToChange.begin();
300  mit != parentsToChange.end(); ++mit) {
301  RefEntity_multiIndex::iterator it = ref_ents_ptr->find(mit->first);
302  bool success = const_cast<RefEntity_multiIndex *>(ref_ents_ptr)
303  ->modify(it, RefEntity_change_parent(mit->second));
304  if (!success) {
305  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
306  "impossible to set parent");
307  }
308  }
310  }
311  };
312  SetParent set_parent;
313 
314  Range ents_to_set_bit;
315 
316  CHKERR m_field.get_ref_ents(&refined_ents_ptr);
317  CHKERR m_field.get_ref_finite_elements(&refined_finite_elements_ptr);
318 
319  Check check;
320  if (debug)
321  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
322 
323  // FIXME: refinement is based on entity handlers, should work on global ids of
324  // nodes, this will allow parallelise algorithm in the future
325 
326  // Find all vertices which parent is edge
327  typedef const RefEntity_multiIndex::index<
328  Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
329  RefEntsByComposite &ref_ents =
330  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
332  ref_parent_ents_view.insert(
333  ref_ents.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE)),
334  ref_ents.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE)));
335  typedef const RefElement_multiIndex::index<Ent_mi_tag>::type RefElementByEnt;
336  RefElementByEnt &ref_finite_element =
337  refined_finite_elements_ptr->get<Ent_mi_tag>();
338  typedef const RefElement_multiIndex_parents_view::index<
339  Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>::type
340  RefEntByParentAndRefEdges;
341  RefElement_multiIndex_parents_view ref_ele_parent_view;
342  ref_ele_parent_view.insert(
343  refined_finite_elements_ptr->get<EntType_mi_tag>().lower_bound(MBTET),
344  refined_finite_elements_ptr->get<EntType_mi_tag>().upper_bound(MBTET));
345  RefEntByParentAndRefEdges &ref_ele_by_parent_and_ref_edges =
346  ref_ele_parent_view
347  .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
348 
349  if (respect_interface) {
350  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
351  "not implemented, set last parameter in refine_TET to false");
352  }
353 
354  Range tets = _tets.subset_by_type(MBTET);
355 
356  std::vector<EntityHandle> parent_ents_refined_and_created;
357  std::vector<EntityHandle> vertices_of_created_tets;
358  std::vector<BitRefEdges> parent_edges_bit_vec;
359  std::vector<int> nb_new_tets_vec;
360  std::vector<int> sub_type_vec;
361 
362  parent_ents_refined_and_created.reserve(tets.size());
363  vertices_of_created_tets.reserve(4 * tets.size());
364  parent_edges_bit_vec.reserve(tets.size());
365  nb_new_tets_vec.reserve(tets.size());
366  sub_type_vec.reserve(tets.size());
367 
368  // make loop over all tets which going to be refined
369  for (auto tit : tets) {
370 
371  if (ref_finite_element.find(tit) == ref_finite_element.end()) {
372  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
373  "this tet is not in refinedFiniteElements");
374  }
375 
376  // get tet connectivity
377  const EntityHandle *conn;
378  int num_nodes;
379  CHKERR moab.get_connectivity(tit, conn, num_nodes, true);
380 
381  // get edges
382  BitRefEdges parent_edges_bit(0);
383  EntityHandle edge_new_nodes[] = {0, 0, 0, 0, 0, 0};
384  int split_edges[] = {-1, -1, -1, -1, -1, -1};
385 
386  for (int ee = 0; ee != 6; ++ee) {
387  EntityHandle edge;
388  CHKERR moab.side_element(tit, 1, ee, edge);
389  RefEntity_multiIndex_view_by_hashed_parent_entity::iterator eit =
390  ref_parent_ents_view.find(edge);
391  if (eit != ref_parent_ents_view.end()) {
392  if (((*eit)->getBitRefLevel() & bit).any()) {
393  edge_new_nodes[ee] = (*eit)->getRefEnt();
394  {
395  const EntityHandle *conn_edge;
396  int num_nodes;
397  moab.get_connectivity(edge, conn_edge, num_nodes, true);
398  if (conn_edge[0] == edge_new_nodes[ee]) {
399  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
400  "node 0 on the edges is mid node, that make no sense");
401  }
402  if (conn_edge[1] == edge_new_nodes[ee]) {
403  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
404  "node 1 on the edges is mid node, that make no sense");
405  }
406  }
407  split_edges[parent_edges_bit.count()] = ee;
408  parent_edges_bit.set(ee, 1);
409  }
410  }
411  }
412 
413  // test if nodes used to refine are not part of tet
414  if (debug) {
415  for (int ee = 0; ee != 6; ee++) {
416  if (edge_new_nodes[ee] == no_handle)
417  continue;
418  for (int nn = 0; nn != 4; nn++) {
419  if (conn[nn] == edge_new_nodes[ee]) {
420  std::cerr << "problem on edge: " << ee << endl;
421  std::cerr << "tet conn : " << conn[0] << " " << conn[1] << " "
422  << conn[2] << " " << conn[3] << " : "
423  << edge_new_nodes[ee] << std::endl;
424  for (int eee = 0; eee != 6; ++eee) {
425  EntityHandle edge;
426  CHKERR moab.side_element(tit, 1, eee, edge);
427  const EntityHandle *conn_edge;
428  int num_nodes;
429  CHKERR moab.get_connectivity(edge, conn_edge, num_nodes, true);
430  std::cerr << eee << " : " << conn_edge[0] << " " << conn_edge[1]
431  << std::endl;
432  }
433  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
434  "nodes used to refine are not part of tet");
435  }
436  }
437  }
438  }
439 
440  // swap nodes forward
441  EntityHandle _conn_[4];
442  std::copy(&conn[0], &conn[4], &_conn_[0]);
443 
444  // build connectivity for ref tets
445  EntityHandle new_tets_conns[8 * 4];
446  std::fill(&new_tets_conns[0], &new_tets_conns[8 * 4], no_handle);
447 
448  int sub_type = -1, nb_new_tets = 0;
449  switch (parent_edges_bit.count()) {
450  case 0: {
451  RefEntity_multiIndex::iterator tit_miit;
452  tit_miit = refined_ents_ptr->find(tit);
453  if (tit_miit == refined_ents_ptr->end())
454  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
455  "can not find this tet");
456  ents_to_set_bit.insert(tit);
457  continue;
458  } break;
459  case 1:
460  sub_type = 0;
461  tet_type_1(_conn_, split_edges[0], edge_new_nodes[split_edges[0]],
462  new_tets_conns);
463  nb_new_tets = 2;
464  break;
465  case 2:
466  sub_type =
467  tet_type_2(_conn_, split_edges, edge_new_nodes, new_tets_conns);
468  if (sub_type & (4 | 8 | 16)) {
469  nb_new_tets = 3;
470  break;
471  } else if (sub_type == 1) {
472  nb_new_tets = 4;
473  break;
474  };
475  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Imposible case");
476  break;
477  case 3:
478  sub_type =
479  tet_type_3(_conn_, split_edges, edge_new_nodes, new_tets_conns);
480  if (sub_type <= 4) {
481  nb_new_tets = 4;
482  break;
483  } else if (sub_type <= 7) {
484  nb_new_tets = 5;
485  break;
486  }
487  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Imposible case");
488  case 4:
489  sub_type =
490  tet_type_4(_conn_, split_edges, edge_new_nodes, new_tets_conns);
491  if (sub_type == 0) {
492  nb_new_tets = 5;
493  break;
494  } else if (sub_type <= 7) {
495  nb_new_tets = 6;
496  break;
497  }
498  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Imposible case");
499  case 5:
500  sub_type = tet_type_5(moab, _conn_, edge_new_nodes, new_tets_conns);
501  nb_new_tets = 7;
502  break;
503  case 6:
504  sub_type = 0;
505  tet_type_6(moab, _conn_, edge_new_nodes, new_tets_conns);
506  nb_new_tets = 8;
507  break;
508  default:
509  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Imposible case");
510  }
511 
512  // find that tets
513  auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.equal_range(
514  boost::make_tuple(tit, parent_edges_bit.to_ulong()));
515  // check if tet with this refinement scheme already exits
516  std::vector<EntityHandle> ents_to_set_bit_vec;
517  if (std::distance(it_by_ref_edges.first, it_by_ref_edges.second) ==
518  static_cast<size_t>(nb_new_tets)) {
519  ents_to_set_bit_vec.reserve(nb_new_tets);
520  for (int tt = 0; it_by_ref_edges.first != it_by_ref_edges.second;
521  it_by_ref_edges.first++, tt++) {
522  auto tet = (*it_by_ref_edges.first)->getRefEnt();
523  ents_to_set_bit_vec.emplace_back(tet);
524  // set ref tets entities
525  if (debug) {
526  // add this tet if exist to this ref level
527  RefEntity_multiIndex::iterator ref_tet_it;
528  ref_tet_it = refined_ents_ptr->find(tet);
529  if (ref_tet_it == refined_ents_ptr->end()) {
530  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
531  "Tetrahedron should be in database");
532  }
533  }
534  }
535  ents_to_set_bit.insert_list(ents_to_set_bit_vec.begin(),
536  ents_to_set_bit_vec.end());
537 
538  } else {
539  // if this element was not refined or was refined with different patterns
540  // of split edges create new elements
541 
542  parent_ents_refined_and_created.emplace_back(tit);
543  for (int tt = 0; tt != nb_new_tets; ++tt) {
544  for (auto nn : {0, 1, 2, 3})
545  vertices_of_created_tets.emplace_back(new_tets_conns[4 * tt + nn]);
546  }
547  parent_edges_bit_vec.emplace_back(parent_edges_bit);
548  nb_new_tets_vec.emplace_back(nb_new_tets);
549  sub_type_vec.emplace_back(sub_type);
550  }
551  }
552 
553  // Create tets
554  EntityHandle start_e = 0;
555  EntityHandle *conn_moab;
556  const int nb_new_tets = vertices_of_created_tets.size() / 4;
557  read_util->get_element_connect(nb_new_tets, 4, MBTET, 0, start_e, conn_moab);
558  std::copy(vertices_of_created_tets.begin(), vertices_of_created_tets.end(),
559  conn_moab);
560  CHKERR read_util->update_adjacencies(start_e, nb_new_tets, 4, conn_moab);
561  ents_to_set_bit.insert(start_e, start_e + nb_new_tets - 1);
562 
563  // Create adj entities
564  for (auto d : {1, 2}) {
565  Range ents;
566  CHKERR moab.get_adjacencies(ents_to_set_bit, d, true, ents,
567  moab::Interface::UNION);
568  }
569 
570  // Set parrents and adjacencies
571  for (int idx = 0; idx != parent_ents_refined_and_created.size(); ++idx) {
572 
573  const EntityHandle tit = parent_ents_refined_and_created[idx];
574  const BitRefEdges &parent_edges_bit = parent_edges_bit_vec[idx];
575  const int nb_new_tets = nb_new_tets_vec[idx];
576  const int sub_type = sub_type_vec[idx];
577 
578  std::array<EntityHandle, 8> ref_tets;
579  for (int tt = 0; tt != nb_new_tets; ++tt, ++start_e)
580  ref_tets[tt] = start_e;
581 
582  if (nb_new_tets) {
583 
584  int ref_type[2];
585  ref_type[0] = parent_edges_bit.count();
586  ref_type[1] = sub_type;
587  for (int tt = 0; tt != nb_new_tets; ++tt) {
588  CHKERR moab.tag_set_data(cOre.get_th_RefType(), &ref_tets[tt], 1,
589  ref_type);
590  CHKERR moab.tag_set_data(cOre.get_th_RefBitEdge(), &ref_tets[tt], 1,
591  &parent_edges_bit);
592  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &ref_tets[tt],
593  1, &tit);
594  }
595 
596  // hash map of nodes (RefEntity) by edges (EntityHandle)
597  std::map<EntityHandle /*edge*/, EntityHandle /*node*/>
598  map_ref_nodes_by_edges;
599 
600  Range tet_edges;
601  CHKERR moab.get_adjacencies(&tit, 1, 1, false, tet_edges);
602  for (auto edge : tet_edges) {
603  RefEntity_multiIndex_view_by_hashed_parent_entity::iterator eit =
604  ref_parent_ents_view.find(edge);
605  if (eit != ref_parent_ents_view.end()) {
606  if (((*eit)->getBitRefLevel() & bit).any()) {
607  map_ref_nodes_by_edges[(*eit)->getParentEnt()] =
608  eit->get()->getRefEnt();
609  }
610  }
611  }
612 
613  // find parents for new edges and faces
614  // get tet edges and faces
615  Range tit_edges, tit_faces;
616  CHKERR moab.get_adjacencies(&tit, 1, 1, false, tit_edges);
617  CHKERR moab.get_adjacencies(&tit, 1, 2, false, tit_faces);
618  Range edges_nodes[6], faces_nodes[4];
619  // for edges - add ref nodes
620  // edges_nodes[ee] - contains all nodes on edge ee including mid nodes if
621  // exist
622  Range::iterator eit = tit_edges.begin();
623  for (int ee = 0; eit != tit_edges.end(); eit++, ee++) {
624  CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee], true);
625  std::map<EntityHandle, EntityHandle>::iterator map_miit =
626  map_ref_nodes_by_edges.find(*eit);
627  if (map_miit != map_ref_nodes_by_edges.end()) {
628  edges_nodes[ee].insert(map_miit->second);
629  }
630  }
631  // for faces - add ref nodes
632  // faces_nodes[ff] - contains all nodes on face ff including mid nodes if
633  // exist
634  Range::iterator fit = tit_faces.begin();
635  for (int ff = 0; fit != tit_faces.end(); fit++, ff++) {
636  CHKERR moab.get_connectivity(&*fit, 1, faces_nodes[ff], true);
637  // Get edges on face and loop over those edges to add mid-nodes to range
638  Range fit_edges;
639  CHKERR moab.get_adjacencies(&*fit, 1, 1, false, fit_edges);
640  for (Range::iterator eit2 = fit_edges.begin(); eit2 != fit_edges.end();
641  eit2++) {
642  std::map<EntityHandle, EntityHandle>::iterator map_miit =
643  map_ref_nodes_by_edges.find(*eit2);
644  if (map_miit != map_ref_nodes_by_edges.end()) {
645  faces_nodes[ff].insert(map_miit->second);
646  }
647  }
648  }
649  // add ref nodes to tet
650  // tet_nodes contains all nodes on tet including mid edge nodes
651  Range tet_nodes;
652  CHKERR moab.get_connectivity(&tit, 1, tet_nodes, true);
653  for (std::map<EntityHandle, EntityHandle>::iterator map_miit =
654  map_ref_nodes_by_edges.begin();
655  map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
656  tet_nodes.insert(map_miit->second);
657  }
658  Range ref_edges;
659  // Get all all edges of refined tets
660  CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 1, false,
661  ref_edges, moab::Interface::UNION);
662  // Check for all ref edge and set parents
663  for (Range::iterator reit = ref_edges.begin(); reit != ref_edges.end();
664  reit++) {
665  Range ref_edges_nodes;
666  CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes, true);
667  if (ref_edges_nodes.size() != 2) {
668  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
669  "data inconsistency, edge should have 2 nodes");
670  }
671  // Check if ref edge is an coarse edge (loop over coarse tet edges)
672  int ee = 0;
673  for (; ee < 6; ee++) {
674  // Two nodes are common (node[0],node[1],ref_node (if exist))
675  // this tests if given edge is contained by edge of refined
676  // tetrahedral
677  if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
678  EntityHandle edge = tit_edges[ee];
679  CHKERR set_parent(*reit, edge, refined_ents_ptr, cOre);
680  break;
681  }
682  }
683  if (ee < 6)
684  continue; // this refined edge is contained by edge of tetrahedral
685  // check if ref edge is in coarse face
686  int ff = 0;
687  for (; ff < 4; ff++) {
688  // two nodes are common (node[0],node[1],ref_node (if exist))
689  // this tests if given face is contained by face of tetrahedral
690  if (intersect(faces_nodes[ff], ref_edges_nodes).size() == 2) {
691  EntityHandle face = tit_faces[ff];
692  CHKERR set_parent(*reit, face, refined_ents_ptr, cOre);
693  break;
694  }
695  }
696  if (ff < 4)
697  continue; // this refined edge is contained by face of tetrahedral
698 
699  // check if ref edge is in coarse tetrahedral (i.e. that is internal
700  // edge of refined tetrahedral)
701  if (intersect(tet_nodes, ref_edges_nodes).size() == 2) {
702  CHKERR set_parent(*reit, tit, refined_ents_ptr, cOre);
703  continue;
704  }
705 
706  // Refined edge is not child of any edge, face or tetrahedral, this is
707  // imposible edge
708  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
709  "impossible refined edge");
710  }
711 
712  Range ref_faces;
713  CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 2, false,
714  ref_faces, moab::Interface::UNION);
715  Tag th_interface_side;
716  const int def_side[] = {0};
717  CHKERR moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
718  th_interface_side,
719  MB_TAG_CREAT | MB_TAG_SPARSE, def_side);
720  // Check for all ref faces
721  for (auto rfit : ref_faces) {
722  Range ref_faces_nodes;
723  CHKERR moab.get_connectivity(&rfit, 1, ref_faces_nodes, true);
724  // Check if ref face is in coarse face
725  int ff = 0;
726  for (; ff < 4; ff++) {
727  // Check if refined triangle is contained by face of tetrahedral
728  if (intersect(faces_nodes[ff], ref_faces_nodes).size() == 3) {
729  EntityHandle face = tit_faces[ff];
730  CHKERR set_parent(rfit, face, refined_ents_ptr, cOre);
731  int side = 0;
732  // Set face side if it is on interface
733  CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
734  CHKERR moab.tag_set_data(th_interface_side, &rfit, 1, &side);
735  break;
736  }
737  }
738  if (ff < 4)
739  continue; // this face is contained by one of tetrahedrons
740  // check if ref face is in coarse tetrahedral
741  // this is ref face which is contained by tetrahedral volume
742  if (intersect(tet_nodes, ref_faces_nodes).size() == 3) {
743  CHKERR set_parent(rfit, tit, refined_ents_ptr, cOre);
744  continue;
745  }
746  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
747  "impossible refined face");
748  }
749  }
750  }
751 
752  if (debug)
753  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
754  CHKERR set_parent(refined_ents_ptr);
755  if (debug)
756  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
757  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(ents_to_set_bit,
758  bit, true, verb);
759  if (debug)
760  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
761 
763 }
virtual MoFEMErrorCode get_ref_ents(const RefEntity_multiIndex **refined_ents_ptr) const =0
Get ref entities multi-index from database.
Tag get_th_RefParentHandle() const
Definition: Core.hpp:150
int tet_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
virtual moab::Interface & get_moab()=0
multi_index_container< boost::shared_ptr< RefElement >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getRefEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityType, &RefElement::getEntType > > > > RefElement_multiIndex
moab::Interface & get_moab()
Definition: Core.hpp:266
int tet_type_5(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition: Common.hpp:27
void tet_type_1(const EntityHandle *conn, const int split_edge, const EntityHandle edge_new_node, EntityHandle *new_tets_conn)
Tag get_th_RefBitEdge() const
Definition: Core.hpp:152
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
std::bitset< BITREFEDGES_SIZE > BitRefEdges
Definition: Types.hpp:45
void tet_type_6(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
static const bool debug
multi_index_container< boost::shared_ptr< RefElement >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getRefEnt > >, 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
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::getRefEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_hashed_parent_entity
multi-index view of RefEntity by parent entity
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode get_ref_finite_elements(const RefElement_multiIndex **refined_finite_elements_ptr) const =0
Get ref finite elements multi-index form database.
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, 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::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
int tet_type_3(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
virtual MPI_Comm & get_comm() const =0
Tag get_th_RefType() const
Definition: Core.hpp:153
int tet_type_4(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)

Member Data Documentation

◆ cOre

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

Definition at line 46 of file MeshRefinement.hpp.


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