v0.10.0
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 ()=default
 
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_generate_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 moab::Interface & get_moab()=0
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
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415

◆ 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  auto refined_ents_ptr = m_field.get_ref_ents();
97  auto miit =
98  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
99  .lower_bound(boost::make_tuple(MBVERTEX, MBEDGE));
100  auto hi_miit =
101  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
102  .upper_bound(boost::make_tuple(MBVERTEX, MBEDGE));
104  ref_parent_ents_view.insert(miit, hi_miit);
105 
106  Range edges = ents.subset_by_type(MBEDGE);
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  for (auto hi_miit_view = ref_parent_ents_view.upper_bound(p_eit->second);
135  miit_view != hi_miit_view; ++miit_view) {
136  edge_having_parent_vertex.insert(edge_having_parent_vertex.end(),
137  miit_view->get()->getParentEnt());
138  add_bit.insert(add_bit.end(), miit_view->get()->getEnt());
139  }
140  }
141 
142  Range add_vertex_edges =
143  subtract(Range(p_eit->first, p_eit->second), edge_having_parent_vertex);
144 
145  for (auto e : add_vertex_edges)
146  parent_edge.emplace_back(e);
147 
148  for (auto e : add_vertex_edges) {
149 
150  const EntityHandle *conn;
151  int num_nodes;
152  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
153  if (PetscUnlikely(num_nodes != 2)) {
154  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
155  "edge should have 2 edges");
156  }
157  CHKERR moab.get_coords(conn, num_nodes, coords.data());
158  t_0(i) += t_1(i);
159  t_0(i) *= 0.5;
160 
161  for (auto j : {0, 1, 2})
162  vert_coords[j].emplace_back(t_0(j));
163  }
164  }
165 
166  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit, bit);
167 
168  if (!vert_coords[0].empty()) {
169  ReadUtilIface *read_util;
170  CHKERR moab.query_interface(read_util);
171  int num_nodes = vert_coords[0].size();
172  vector<double *> arrays_coord;
173  CHKERR read_util->get_node_coords(3, num_nodes, 0, start_v, arrays_coord);
174  Range verts(start_v, start_v + num_nodes - 1);
175  for (auto dd : {0, 1, 2}) {
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 }
virtual moab::Interface & get_moab()=0
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_non_unique< const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< boost::shared_ptr< RefEntity >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >> > > RefEntity_multiIndex_view_by_ordered_parent_entity
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
Tag get_th_RefParentHandle() const
Definition: Core.hpp:192
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
static Index< 'i', 3 > i
static Index< 'j', 3 > j
#define CHKERR
Inline error check.
Definition: definitions.h:604
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415

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

956  {
957  Interface &m_field = cOre;
958  auto refined_ents_ptr = m_field.get_ref_ents();
960  typedef const RefEntity_multiIndex::index<Ent_mi_tag>::type RefEntsByEnt;
961  RefEntsByEnt::iterator miit = refined_ents_ptr->find(meshset);
962  if (miit == refined_ents_ptr->end()) {
963  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
964  "this meshset is not in ref database");
965  }
966  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
967  meshset, bit, meshset, MBEDGE, recursive, verb);
968  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
969  meshset, bit, meshset, MBTRI, recursive, verb);
970  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
971  meshset, bit, meshset, MBTET, recursive, verb);
972  *(const_cast<RefEntity *>(miit->get())->getBitRefLevelPtr()) |= bit;
974 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
MoFEMErrorCode get_ref_ents(const RefEntity_multiIndex **refined_entities_ptr) const
Get ref entities multi-index from database.
Definition: Core.cpp:802
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415

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

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

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

191  {
192  Interface &m_field = cOre;
193  moab::Interface &moab = m_field.get_moab();
195  Range tets;
196  CHKERR moab.get_entities_by_type(meshset, MBTET, tets, false);
197  CHKERR refine_TET(tets, bit, respect_interface, verb, ref_edges_ptr, debug);
199 }
virtual moab::Interface & get_moab()=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
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:604
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415

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

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

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: