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

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_middel_of_edges(meshset, bit, recursive, verb,
77  start_v);
78  }
virtual MoFEMErrorCode add_vertices_in_the_middel_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_middel_of_edges(edges, bit, verb, start_v);
105  }
virtual MoFEMErrorCode add_vertices_in_the_middel_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_middel_of_edges() [1/2]

MoFEMErrorCode MoFEM::MeshRefinement::add_vertices_in_the_middel_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_middel_of_edges(edges, bit, verb, start_v);
90 }
virtual MoFEMErrorCode add_vertices_in_the_middel_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:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ add_vertices_in_the_middel_of_edges() [2/2]

MoFEMErrorCode MoFEM::MeshRefinement::add_vertices_in_the_middel_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  Range edges = ents.subset_by_type(MBEDGE);
99  auto miit =
100  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
101  .lower_bound(boost::make_tuple(MBVERTEX, MBEDGE));
102  auto hi_miit =
103  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
104  .upper_bound(boost::make_tuple(MBVERTEX, MBEDGE));
106  ref_parent_ents_view.insert(miit, hi_miit);
107  if (verb >= VERBOSE) {
108  std::ostringstream ss;
109  ss << "ref level " << bit << " nb. edges to refine " << edges.size()
110  << std::endl;
111  PetscPrintf(m_field.get_comm(), ss.str().c_str());
112  }
113 
114  std::array<std::vector<double>, 3> vert_coords;
115  for (auto &vc : vert_coords)
116  vc.reserve(edges.size());
117 
118  std::vector<EntityHandle> parent_edge;
119  parent_edge.reserve(edges.size());
120 
121  std::array<double, 6> coords;
123  &coords[0], &coords[1], &coords[2]};
125  &coords[3], &coords[4], &coords[5]};
127 
128  Range add_bit;
129  for (auto p_eit = edges.pair_begin(); p_eit != edges.pair_end(); ++p_eit) {
130  auto miit_view = ref_parent_ents_view.lower_bound(p_eit->first);
131 
132  Range edge_having_parent_vertex;
133  if (miit_view != ref_parent_ents_view.end()) {
134  auto hi_miit_view = ref_parent_ents_view.upper_bound(p_eit->second);
135  for (; miit_view != hi_miit_view; ++miit_view)
136  edge_having_parent_vertex.insert(miit_view->get()->getParentEnt());
137  }
138 
139  Range add_vertex_edges =
140  subtract(Range(p_eit->first, p_eit->second), edge_having_parent_vertex);
141 
142  for (auto e : add_vertex_edges)
143  parent_edge.emplace_back(e);
144 
145  for (auto e : add_vertex_edges) {
146 
147  const EntityHandle *conn;
148  int num_nodes;
149  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
150  if (PetscUnlikely(num_nodes != 2)) {
151  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
152  "edge should have 2 edges");
153  }
154  CHKERR moab.get_coords(conn, num_nodes, coords.data());
155  t_0(i) += t_1(i);
156  t_0(i) *= 0.5;
157 
158  for (auto j : {0, 1, 2})
159  vert_coords[j].emplace_back(t_0(j));
160  }
161  }
162 
163  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit, bit);
164  if (!vert_coords[0].empty()) {
165  int num_nodes = vert_coords[0].size();
166  vector<double *> arrays_coord;
167  EntityHandle startv;
168  {
169  ReadUtilIface *iface;
170  CHKERR moab.query_interface(iface);
171  CHKERR iface->get_node_coords(3, num_nodes, start_v, startv,
172  arrays_coord);
173  }
174  Range verts(startv, startv + num_nodes - 1);
175  for (int dd = 0; dd != 3; ++dd) {
176  std::copy(vert_coords[dd].begin(), vert_coords[dd].end(),
177  arrays_coord[dd]);
178  }
179  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), verts,
180  &*parent_edge.begin());
181  CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
182  verts, bit, verb);
183 }
185 } // namespace MoFEM
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:476
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:595
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:406

◆ 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:500
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
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 896 of file MeshRefinement.cpp.

898  {
899  Interface &m_field = cOre;
900  const RefEntity_multiIndex *refined_ents_ptr;
902  CHKERR m_field.get_ref_ents(&refined_ents_ptr);
903  typedef const RefEntity_multiIndex::index<Ent_mi_tag>::type RefEntsByEnt;
904  RefEntsByEnt::iterator miit = refined_ents_ptr->find(meshset);
905  if (miit == refined_ents_ptr->end()) {
906  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
907  "this meshset is not in ref database");
908  }
909  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
910  meshset, bit, meshset, MBEDGE, recursive, verb);
911  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
912  meshset, bit, meshset, MBTRI, recursive, verb);
913  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
914  meshset, bit, meshset, MBTET, recursive, verb);
915  *(const_cast<RefEntity *>(miit->get())->getBitRefLevelPtr()) |= bit;
917 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
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:406

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

691  {
692 
693  Interface &m_field = cOre;
694  moab::Interface &moab = m_field.get_moab();
695  const RefEntity_multiIndex *refined_ents_ptr;
696  const RefElement_multiIndex *refined_finite_elements_ptr;
697 
698  // FIXME: refinement is based on entity handlers, should work on global ids of
699  // nodes, this will allow parallelise algorithm in the future
700 
702 
703  ierr = m_field.get_ref_ents(&refined_ents_ptr);
704  CHKERRG(ierr);
705  ierr = m_field.get_ref_finite_elements(&refined_finite_elements_ptr);
706 
707  typedef const RefEntity_multiIndex::index<Ent_mi_tag>::type RefEntsByEnt;
708 
709  typedef const RefElement_multiIndex_parents_view::index<
710  Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>::type
711  RefEntByParentAndRefEdges;
712  RefElement_multiIndex_parents_view ref_ele_parent_view;
713  ref_ele_parent_view.insert(
714  refined_finite_elements_ptr->get<EntType_mi_tag>().lower_bound(MBPRISM),
715  refined_finite_elements_ptr->get<EntType_mi_tag>().upper_bound(MBPRISM));
716  RefEntByParentAndRefEdges &ref_ele_by_parent_and_ref_edges =
717  ref_ele_parent_view
718  .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
719  // find all vertices which parent is edge
720  typedef const RefEntity_multiIndex::index<
721  Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
722  RefEntsByComposite &ref_ents_by_comp =
723  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
725  ref_parent_ents_view.insert(
726  ref_ents_by_comp.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE)),
727  ref_ents_by_comp.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE)));
728  Range prisms;
729  rval = moab.get_entities_by_type(meshset, MBPRISM, prisms, false);
731  Range::iterator pit = prisms.begin();
732  for (; pit != prisms.end(); pit++) {
733  RefEntsByEnt::iterator miit_prism =
734  refined_ents_ptr->get<Ent_mi_tag>().find(*pit);
735  if (miit_prism == refined_ents_ptr->end()) {
736  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
737  "this prism is not in ref database");
738  }
739  if (verb >= NOISY) {
740  std::ostringstream ss;
741  ss << "ref prism " << **miit_prism << std::endl;
742  PetscPrintf(m_field.get_comm(), ss.str().c_str());
743  }
744  // prism connectivity
745  int num_nodes;
746  const EntityHandle *conn;
747  rval = moab.get_connectivity(*pit, conn, num_nodes, true);
749  assert(num_nodes == 6);
750  // edges connectivity
751  EntityHandle edges[6];
752  for (int ee = 0; ee < 3; ee++) {
753  rval = moab.side_element(*pit, 1, ee, edges[ee]);
755  }
756  for (int ee = 6; ee < 9; ee++) {
757  rval = moab.side_element(*pit, 1, ee, edges[ee - 3]);
759  }
760  // detect split edges
761  BitRefEdges split_edges(0);
762  EntityHandle edge_nodes[6];
763  std::fill(&edge_nodes[0], &edge_nodes[6], no_handle);
764  for (int ee = 0; ee < 6; ee++) {
765  RefEntity_multiIndex_view_by_hashed_parent_entity::iterator miit_view =
766  ref_parent_ents_view.find(edges[ee]);
767  if (miit_view != ref_parent_ents_view.end()) {
768  if (((*miit_view)->getBitRefLevel() & bit).any()) {
769  edge_nodes[ee] = (*miit_view)->getRefEnt();
770  split_edges.set(ee);
771  }
772  }
773  }
774  if (split_edges.count() == 0) {
775  *(const_cast<RefEntity *>(miit_prism->get())->getBitRefLevelPtr()) |= bit;
776  if (verb >= VERY_NOISY)
777  PetscPrintf(m_field.get_comm(), "no refinement");
778  continue;
779  }
780  // check consistency
781  if (verb >= NOISY) {
782  std::ostringstream ss;
783  ss << "prism split edges " << split_edges << " count "
784  << split_edges.count() << std::endl;
785  PetscPrintf(m_field.get_comm(), ss.str().c_str());
786  }
787  // prism ref
788  EntityHandle new_prism_conn[4 * 6];
789  std::fill(&new_prism_conn[0], &new_prism_conn[4 * 6], no_handle);
790  int nb_new_prisms = 0;
791  switch (split_edges.count()) {
792  case 0:
793  break;
794  case 2:
795  ierr = prism_type_1(conn, split_edges, edge_nodes, new_prism_conn);
796  CHKERRG(ierr);
797  nb_new_prisms = 2;
798  break;
799  case 4:
800  ierr = prism_type_2(conn, split_edges, edge_nodes, new_prism_conn);
801  CHKERRG(ierr);
802  nb_new_prisms = 3;
803  break;
804  case 6:
805  ierr = prism_type_3(conn, split_edges, edge_nodes, new_prism_conn);
806  CHKERRG(ierr);
807  nb_new_prisms = 4;
808  break;
809  default:
810  std::ostringstream ss;
811  ss << split_edges << " : [ " << conn[0] << " " << conn[1] << " "
812  << conn[2] << " " << conn[3] << " " << conn[4] << " " << conn[5]
813  << " ]";
814  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, ss.str().c_str());
815  }
816  // find that prism
817  std::bitset<4> ref_prism_bit(0);
818  RefEntByParentAndRefEdges::iterator it_by_ref_edges =
819  ref_ele_by_parent_and_ref_edges.lower_bound(
820  boost::make_tuple(*pit, split_edges.to_ulong()));
821  RefEntByParentAndRefEdges::iterator hi_it_by_ref_edges =
822  ref_ele_by_parent_and_ref_edges.upper_bound(
823  boost::make_tuple(*pit, split_edges.to_ulong()));
824  RefEntByParentAndRefEdges::iterator it_by_ref_edges2 = it_by_ref_edges;
825  for (int pp = 0; it_by_ref_edges2 != hi_it_by_ref_edges;
826  it_by_ref_edges2++, pp++) {
827  // Add this tet to this ref
828  *(const_cast<RefElement *>(it_by_ref_edges2->get())
829  ->getBitRefLevelPtr()) |= bit;
830  ref_prism_bit.set(pp, 1);
831  if (verb > 2) {
832  std::ostringstream ss;
833  ss << "is refined " << *it_by_ref_edges2 << std::endl;
834  PetscPrintf(m_field.get_comm(), ss.str().c_str());
835  }
836  }
837  if (it_by_ref_edges != hi_it_by_ref_edges) {
838  if (ref_prism_bit.count() != (unsigned int)nb_new_prisms)
839  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
840  "data inconsistency");
841  } else {
842  EntityHandle ref_prisms[4];
843  // create prism
844  for (int pp = 0; pp < nb_new_prisms; pp++) {
845  if (verb > 3) {
846  std::ostringstream ss;
847  ss << "ref prism " << ref_prism_bit << std::endl;
848  PetscPrintf(m_field.get_comm(), ss.str().c_str());
849  }
850  if (!ref_prism_bit.test(pp)) {
851  rval = moab.create_element(MBPRISM, &new_prism_conn[6 * pp], 6,
852  ref_prisms[pp]);
854  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(),
855  &ref_prisms[pp], 1, &*pit);
857  rval = moab.tag_set_data(cOre.get_th_RefBitLevel(), &ref_prisms[pp],
858  1, &bit);
860  rval = moab.tag_set_data(cOre.get_th_RefBitEdge(), &ref_prisms[pp], 1,
861  &split_edges);
863  std::pair<RefEntity_multiIndex::iterator, bool> p_ent =
864  const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
865  ->insert(boost::shared_ptr<RefEntity>(new RefEntity(
866  m_field.get_basic_entity_data_ptr(), ref_prisms[pp])));
867  std::pair<RefElement_multiIndex::iterator, bool> p_fe;
868  try {
869  p_fe =
870  const_cast<RefElement_multiIndex *>(refined_finite_elements_ptr)
871  ->insert(boost::shared_ptr<RefElement>(
872  new RefElement_PRISM(*p_ent.first)));
873  } catch (MoFEMException const &e) {
874  SETERRQ(PETSC_COMM_SELF, e.errorCode, e.errorMessage);
875  }
876  ref_prism_bit.set(pp);
877  ierr = cOre.addPrismToDatabase(ref_prisms[pp]);
878  CHKERRG(ierr);
879  if (verb > 2) {
880  std::ostringstream ss;
881  ss << "add prism: " << **p_fe.first << std::endl;
882  if (verb > 7) {
883  for (int nn = 0; nn < 6; nn++) {
884  ss << new_prism_conn[nn] << " ";
885  }
886  ss << std::endl;
887  }
888  PetscPrintf(m_field.get_comm(), ss.str().c_str());
889  }
890  }
891  }
892  }
893  }
895 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:514
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 MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:500
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition: Common.hpp:27
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:543
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
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
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:266
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
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
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

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

190  {
191  Interface &m_field = cOre;
192  moab::Interface &moab = m_field.get_moab();
194  Range tets;
195  CHKERR moab.get_entities_by_type(meshset, MBTET, tets, false);
196  CHKERR refine_TET(tets, bit, respect_interface, verb, ref_edges_ptr);
198 }
virtual MoFEMErrorCode refine_TET(const EntityHandle meshset, const BitRefLevel &bit, const bool respect_interface=false, int verb=QUIET, Range *ref_edges=NULL)
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:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

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

203  {
204 
205  Interface &m_field = cOre;
206  moab::Interface &moab = m_field.get_moab();
207  const RefEntity_multiIndex *refined_ents_ptr;
208  const RefElement_multiIndex *refined_finite_elements_ptr;
210 
211  // Check if refinement is correct
212  struct Check {
213  map<EntityHandle, EntityHandle> entParentMap;
214  MoFEMErrorCode operator()(Range *ref_edges,
215  const RefEntity_multiIndex *ref_ents_ptr,
216  MoFEM::Core &core) {
218  if (!ref_edges)
220  for (Range::iterator eit = ref_edges->begin(); eit != ref_edges->end();
221  ++eit) {
222  RefEntity_multiIndex::index<
223  Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vit =
224  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
225  boost::make_tuple(MBVERTEX, *eit));
226  if (vit ==
227  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
228  RefEntity_multiIndex::iterator e_eit = ref_ents_ptr->find(*eit);
229  if (e_eit == ref_ents_ptr->end()) {
230  SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
231  "Edge not found %ld", *eit);
232  }
233  cerr << "Parent edge" << endl << **e_eit << endl;
234  if (entParentMap.find(*eit) != entParentMap.end()) {
235  RefEntity_multiIndex::iterator v_eit =
236  ref_ents_ptr->find(entParentMap[*eit]);
237  if (v_eit == ref_ents_ptr->end()) {
238  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
239  "Vertex not found");
240  }
241  cerr << "Vertex " << **v_eit << endl;
242  }
243  RefEntity_multiIndex::index<Ent_Ent_mi_tag>::type::iterator ee_it,
244  ee_hi_it;
245  ee_it = ref_ents_ptr->get<Ent_Ent_mi_tag>().lower_bound(*eit);
246  ee_hi_it = ref_ents_ptr->get<Ent_Ent_mi_tag>().upper_bound(*eit);
247  for (; ee_it != ee_hi_it; ++ee_it) {
248  cerr << "Ent having edge parent by parent " << **ee_it << endl;
249  }
250  RefEntity_multiIndex tmp_index;
251  tmp_index.insert(ref_ents_ptr->begin(), ref_ents_ptr->end());
252  RefEntity_multiIndex::index<
253  Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vvit =
254  tmp_index.get<Composite_ParentEnt_And_EntType_mi_tag>().find(
255  boost::make_tuple(MBVERTEX, *eit));
256  if (vvit !=
257  tmp_index.get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
258  cerr << "Tmp idx Vertex " << **vvit << endl;
259  }
260  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
261  "No vertex on trim edges, that make no sense");
262  } else {
263  entParentMap[vit->get()->getParentEnt()] = vit->get()->getRefEnt();
264  }
265  }
267  }
268  };
269 
270  struct SetParent {
271  map<EntityHandle, EntityHandle> parentsToChange;
272  MoFEMErrorCode operator()(const EntityHandle ent, const EntityHandle parent,
273  const RefEntity_multiIndex *ref_ents_ptr,
274  MoFEM::Core &cOre) {
275  MoFEM::Interface &m_field = cOre;
277  RefEntity_multiIndex::iterator it = ref_ents_ptr->find(ent);
278  if (it != ref_ents_ptr->end()) {
279  if (it->get()->getParentEnt() != parent && ent != parent) {
280  parentsToChange[ent] = parent;
281  }
282  } else {
283  if (ent != parent) {
284  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
285  &ent, 1, &parent);
286  }
287  }
289  }
290  MoFEMErrorCode operator()(const RefEntity_multiIndex *ref_ents_ptr) {
292  for (map<EntityHandle, EntityHandle>::iterator mit =
293  parentsToChange.begin();
294  mit != parentsToChange.end(); ++mit) {
295  RefEntity_multiIndex::iterator it = ref_ents_ptr->find(mit->first);
296  bool success = const_cast<RefEntity_multiIndex *>(ref_ents_ptr)
297  ->modify(it, RefEntity_change_parent(mit->second));
298  if (!success) {
299  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
300  "impossible to set parent");
301  }
302  }
304  }
305  };
306  SetParent set_parent;
307 
308  Range ents_to_set_bit;
309 
310  CHKERR m_field.get_ref_ents(&refined_ents_ptr);
311  CHKERR m_field.get_ref_finite_elements(&refined_finite_elements_ptr);
312 
313  Check check;
314  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
315 
316  // FIXME: refinement is based on entity handlers, should work on global ids of
317  // nodes, this will allow parallelise algorithm in the future
318 
319  // Find all vertices which parent is edge
320  typedef const RefEntity_multiIndex::index<
321  Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
322  RefEntsByComposite &ref_ents =
323  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
325  ref_parent_ents_view.insert(
326  ref_ents.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE)),
327  ref_ents.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE)));
328  typedef const RefElement_multiIndex::index<Ent_mi_tag>::type RefElementByEnt;
329  RefElementByEnt &ref_finite_element =
330  refined_finite_elements_ptr->get<Ent_mi_tag>();
331  typedef const RefElement_multiIndex_parents_view::index<
332  Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>::type
333  RefEntByParentAndRefEdges;
334  RefElement_multiIndex_parents_view ref_ele_parent_view;
335  ref_ele_parent_view.insert(
336  refined_finite_elements_ptr->get<EntType_mi_tag>().lower_bound(MBTET),
337  refined_finite_elements_ptr->get<EntType_mi_tag>().upper_bound(MBTET));
338  RefEntByParentAndRefEdges &ref_ele_by_parent_and_ref_edges =
339  ref_ele_parent_view
340  .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
341 
342  if (respect_interface) {
343  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
344  "not implemented, set last parameter in refine_TET to false");
345  }
346 
347  // make loop over all tets which going to be refined
348  Range tets = _tets.subset_by_type(MBTET);
349  Range::iterator tit = tets.begin();
350  for (; tit != tets.end(); ++tit) {
351 
352  if (ref_finite_element.find(*tit) == ref_finite_element.end()) {
353  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
354  "this tet is not in refinedFiniteElements");
355  }
356 
357  // get tet connectivity
358  const EntityHandle *conn;
359  int num_nodes;
360  CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
361 
362  // get edges
363  BitRefEdges parent_edges_bit(0);
364  EntityHandle edge_new_nodes[] = {0, 0, 0, 0, 0, 0};
365  int split_edges[] = {-1, -1, -1, -1, -1, -1};
366  // hash map of nodes (RefEntity) by edges (EntityHandle)
367  std::map<EntityHandle /*edge*/, EntityHandle /*node*/>
368  map_ref_nodes_by_edges;
369  for (int ee = 0; ee != 6; ++ee) {
370  EntityHandle edge;
371  CHKERR moab.side_element(*tit, 1, ee, edge);
372  RefEntity_multiIndex_view_by_hashed_parent_entity::iterator eit =
373  ref_parent_ents_view.find(edge);
374  if (eit != ref_parent_ents_view.end()) {
375  if (((*eit)->getBitRefLevel() & bit).any()) {
376  edge_new_nodes[ee] = (*eit)->getRefEnt();
377  map_ref_nodes_by_edges[(*eit)->getParentEnt()] =
378  eit->get()->getRefEnt();
379  {
380  const EntityHandle *conn_edge;
381  int num_nodes;
382  moab.get_connectivity(edge, conn_edge, num_nodes, true);
383  if (conn_edge[0] == edge_new_nodes[ee]) {
384  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
385  "node 0 on the edges is mid node, that make no sense");
386  }
387  if (conn_edge[1] == edge_new_nodes[ee]) {
388  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
389  "node 1 on the edges is mid node, that make no sense");
390  }
391  }
392  split_edges[parent_edges_bit.count()] = ee;
393  parent_edges_bit.set(ee, 1);
394  }
395  }
396  }
397 
398  // test if nodes used to refine are not part of tet
399  for (int ee = 0; ee != 6; ee++) {
400  if (edge_new_nodes[ee] == no_handle)
401  continue;
402  for (int nn = 0; nn != 4; nn++) {
403  if (conn[nn] == edge_new_nodes[ee]) {
404  std::cerr << "problem on edge: " << ee << endl;
405  std::cerr << "tet conn : " << conn[0] << " " << conn[1] << " "
406  << conn[2] << " " << conn[3] << " : " << edge_new_nodes[ee]
407  << std::endl;
408  for (int eee = 0; eee != 6; ++eee) {
409  EntityHandle edge;
410  CHKERR moab.side_element(*tit, 1, eee, edge);
411  const EntityHandle *conn_edge;
412  int num_nodes;
413  CHKERR moab.get_connectivity(edge, conn_edge, num_nodes, true);
414  std::cerr << eee << " : " << conn_edge[0] << " " << conn_edge[1]
415  << std::endl;
416  }
417  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
418  "nodes used to refine are not part of tet");
419  }
420  }
421  }
422 
423  // swap nodes forward
424  EntityHandle _conn_[4];
425  std::copy(&conn[0], &conn[4], &_conn_[0]);
426  // build connectivity for rf tets
427  EntityHandle new_tets_conns[8 * 4];
428  std::fill(&new_tets_conns[0], &new_tets_conns[8 * 4], no_handle);
429  int sub_type = -1, nb_new_tets = 0;
430  switch (parent_edges_bit.count()) {
431  case 0: {
432  RefEntity_multiIndex::iterator tit_miit;
433  tit_miit = refined_ents_ptr->find(*tit);
434  if (tit_miit == refined_ents_ptr->end())
435  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
436  "can not find this tet");
437  ents_to_set_bit.insert(*tit);
438  continue;
439  } break;
440  case 1:
441  sub_type = 0;
442  tet_type_1(_conn_, split_edges[0], edge_new_nodes[split_edges[0]],
443  new_tets_conns);
444  nb_new_tets = 2;
445  break;
446  case 2:
447  sub_type =
448  tet_type_2(_conn_, split_edges, edge_new_nodes, new_tets_conns);
449  if (sub_type & (4 | 8 | 16)) {
450  nb_new_tets = 3;
451  break;
452  } else if (sub_type == 1) {
453  nb_new_tets = 4;
454  break;
455  };
456  assert(0);
457  break;
458  case 3:
459  sub_type =
460  tet_type_3(_conn_, split_edges, edge_new_nodes, new_tets_conns);
461  if (sub_type <= 4) {
462  nb_new_tets = 4;
463  break;
464  } else if (sub_type <= 7) {
465  nb_new_tets = 5;
466  break;
467  }
468  assert(0);
469  case 4:
470  sub_type =
471  tet_type_4(_conn_, split_edges, edge_new_nodes, new_tets_conns);
472  if (sub_type == 0) {
473  nb_new_tets = 5;
474  break;
475  } else if (sub_type <= 7) {
476  nb_new_tets = 6;
477  break;
478  }
479  assert(0);
480  case 5:
481  sub_type = tet_type_5(moab, _conn_, edge_new_nodes, new_tets_conns);
482  nb_new_tets = 7;
483  break;
484  case 6:
485  sub_type = 0;
486  tet_type_6(moab, _conn_, edge_new_nodes, new_tets_conns);
487  nb_new_tets = 8;
488  break;
489  default:
490  assert(0);
491  }
492  // find that tets
493  EntityHandle ref_tets[8];
494  std::bitset<8> ref_tets_bit(0);
495  RefEntByParentAndRefEdges::iterator it_by_ref_edges =
496  ref_ele_by_parent_and_ref_edges.lower_bound(
497  boost::make_tuple(*tit, parent_edges_bit.to_ulong()));
498  RefEntByParentAndRefEdges::iterator hi_it_by_ref_edges =
499  ref_ele_by_parent_and_ref_edges.upper_bound(
500  boost::make_tuple(*tit, parent_edges_bit.to_ulong()));
501  // check if tet with this refinement shame already exits
502  if (std::distance(it_by_ref_edges, hi_it_by_ref_edges) ==
503  (unsigned int)nb_new_tets) {
504  for (int tt = 0; it_by_ref_edges != hi_it_by_ref_edges;
505  it_by_ref_edges++, tt++) {
506  EntityHandle tet = it_by_ref_edges->get()->getRefEnt();
507  // set ref tets entities
508  ref_tets[tt] = tet;
509  ref_tets_bit.set(tt, 1);
510  // add this tet if exist to this ref level
511  RefEntity_multiIndex::iterator ref_tet_it;
512  ref_tet_it = refined_ents_ptr->find(tet);
513  if (ref_tet_it == refined_ents_ptr->end()) {
514  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
515  "data inconsistency");
516  }
517  ents_to_set_bit.insert(tet);
518  }
519  } else {
520  // if this element was not refined or was refined with different patterns
521  // of split edges create new elements
522  for (int tt = 0; tt < nb_new_tets; tt++) {
523  if (!ref_tets_bit.test(tt)) {
524  CHKERR moab.create_element(MBTET, &new_tets_conns[4 * tt], 4,
525  ref_tets[tt]);
526  ents_to_set_bit.insert(ref_tets[tt]);
527  Range ref_tet_nodes;
528  CHKERR moab.get_connectivity(&ref_tets[tt], 1, ref_tet_nodes, true);
529  int ref_type[2];
530  ref_type[0] = parent_edges_bit.count();
531  ref_type[1] = sub_type;
532  CHKERR moab.tag_set_data(cOre.get_th_RefType(), &ref_tets[tt], 1,
533  ref_type);
534  CHKERR moab.tag_set_data(cOre.get_th_RefBitEdge(), &ref_tets[tt], 1,
535  &parent_edges_bit);
536  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &ref_tets[tt],
537  1, &*tit);
538  // set bit that this element is now in database
539  ref_tets_bit.set(tt);
540  }
541  }
542  }
543  // find parents for new edges and faces
544  // get tet edges and faces
545  Range tit_edges, tit_faces;
546  CHKERR moab.get_adjacencies(&*tit, 1, 1, false, tit_edges);
547  CHKERR moab.get_adjacencies(&*tit, 1, 2, false, tit_faces);
548  Range edges_nodes[6], faces_nodes[4];
549  // for edges - add ref nodes
550  // edges_nodes[ee] - contains all nodes on edge ee including mid nodes if
551  // exist
552  Range::iterator eit = tit_edges.begin();
553  for (int ee = 0; eit != tit_edges.end(); eit++, ee++) {
554  CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee], true);
555  std::map<EntityHandle, EntityHandle>::iterator map_miit =
556  map_ref_nodes_by_edges.find(*eit);
557  if (map_miit != map_ref_nodes_by_edges.end()) {
558  edges_nodes[ee].insert(map_miit->second);
559  }
560  }
561  // for faces - add ref nodes
562  // faces_nodes[ff] - contains all nodes on face ff including mid nodes if
563  // exist
564  Range::iterator fit = tit_faces.begin();
565  for (int ff = 0; fit != tit_faces.end(); fit++, ff++) {
566  CHKERR moab.get_connectivity(&*fit, 1, faces_nodes[ff], true);
567  // Get edges on face and loop over those edges to add mid-nodes to range
568  Range fit_edges;
569  CHKERR moab.get_adjacencies(&*fit, 1, 1, false, fit_edges);
570  for (Range::iterator eit2 = fit_edges.begin(); eit2 != fit_edges.end();
571  eit2++) {
572  std::map<EntityHandle, EntityHandle>::iterator map_miit =
573  map_ref_nodes_by_edges.find(*eit2);
574  if (map_miit != map_ref_nodes_by_edges.end()) {
575  faces_nodes[ff].insert(map_miit->second);
576  }
577  }
578  }
579  // add ref nodes to tet
580  // tet_nodes contains all nodes on tet including mid edge nodes
581  Range tet_nodes;
582  CHKERR moab.get_connectivity(&*tit, 1, tet_nodes, true);
583  for (std::map<EntityHandle, EntityHandle>::iterator map_miit =
584  map_ref_nodes_by_edges.begin();
585  map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
586  tet_nodes.insert(map_miit->second);
587  }
588  Range ref_edges;
589  // Get all all edges of refined tets
590  CHKERR moab.get_adjacencies(ref_tets, nb_new_tets, 1, true, ref_edges,
591  moab::Interface::UNION);
592  // Check for all ref edge and set parents
593  for (Range::iterator reit = ref_edges.begin(); reit != ref_edges.end();
594  reit++) {
595  Range ref_edges_nodes;
596  CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes, true);
597  if (ref_edges_nodes.size() != 2) {
598  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
599  "data inconsistency, edge should have 2 nodes");
600  }
601  // Check if ref edge is an coarse edge (loop over coarse tet edges)
602  int ee = 0;
603  for (; ee < 6; ee++) {
604  // Two nodes are common (node[0],node[1],ref_node (if exist))
605  // this tests if given edge is contained by edge of refined tetrahedral
606  if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
607  EntityHandle edge = tit_edges[ee];
608  CHKERR set_parent(*reit, edge, refined_ents_ptr, cOre);
609  break;
610  }
611  }
612  if (ee < 6)
613  continue; // this refined edge is contained by edge of tetrahedral
614  // check if ref edge is in coarse face
615  int ff = 0;
616  for (; ff < 4; ff++) {
617  // two nodes are common (node[0],node[1],ref_node (if exist))
618  // this tests if given face is contained by face of tetrahedral
619  if (intersect(faces_nodes[ff], ref_edges_nodes).size() == 2) {
620  EntityHandle face = tit_faces[ff];
621  CHKERR set_parent(*reit, face, refined_ents_ptr, cOre);
622  break;
623  }
624  }
625  if (ff < 4)
626  continue; // this refined edge is contained by face of tetrahedral
627 
628  // check if ref edge is in coarse tetrahedral (i.e. that is internal edge
629  // of refined tetrahedral)
630  if (intersect(tet_nodes, ref_edges_nodes).size() == 2) {
631  CHKERR set_parent(*reit, *tit, refined_ents_ptr, cOre);
632  continue;
633  }
634 
635  // Refined edge is not child of any edge, face or tetrahedral, this is
636  // imposible edge
637  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
638  "impossible refined edge");
639  }
640 
641  Range ref_faces;
642  CHKERR moab.get_adjacencies(ref_tets, nb_new_tets, 2, true, ref_faces,
643  moab::Interface::UNION);
644  Tag th_interface_side;
645  const int def_side[] = {0};
646  CHKERR moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
647  th_interface_side, MB_TAG_CREAT | MB_TAG_SPARSE,
648  def_side);
649  // Check for all ref faces
650  for (Range::iterator rfit = ref_faces.begin(); rfit != ref_faces.end();
651  rfit++) {
652  Range ref_faces_nodes;
653  CHKERR moab.get_connectivity(&*rfit, 1, ref_faces_nodes, true);
654  // Check if ref face is in coarse face
655  int ff = 0;
656  for (; ff < 4; ff++) {
657  // Check if refined triangle is contained by face of tetrahedral
658  if (intersect(faces_nodes[ff], ref_faces_nodes).size() == 3) {
659  EntityHandle face = tit_faces[ff];
660  CHKERR set_parent(*rfit, face, refined_ents_ptr, cOre);
661  int side = 0;
662  // Set face side if it is on interface
663  CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
664  CHKERR moab.tag_set_data(th_interface_side, &*rfit, 1, &side);
665  break;
666  }
667  }
668  if (ff < 4)
669  continue; // this face is contained by one of tetrahedrons
670  // check if ref face is in coarse tetrahedral
671  // this is ref face which is contained by tetrahedral volume
672  if (intersect(tet_nodes, ref_faces_nodes).size() == 3) {
673  CHKERR set_parent(*rfit, *tit, refined_ents_ptr, cOre);
674  continue;
675  }
676  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
677  "impossible refined face");
678  }
679  }
680 
681  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
682  CHKERR set_parent(refined_ents_ptr);
683  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
684  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(ents_to_set_bit,
685  bit, true, verb);
686  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
687 
689 }
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:476
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:507
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.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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:595
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:406
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: