v0.9.1
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:483
#define CHKERR
Inline error check.
Definition: definitions.h:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413

◆ 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()->getRefEnt());
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 }
Tag get_th_RefParentHandle() const
Definition: Core.hpp:150
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26

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

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

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

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

◆ 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:483
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413

◆ 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()->getRefEnt();
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<EntType_mi_tag>().lower_bound(MBTET),
340  refined_finite_elements_ptr->get<EntType_mi_tag>().upper_bound(MBTET));
341  RefEntByParentAndRefEdges &ref_ele_by_parent_and_ref_edges =
342  ref_ele_parent_view
343  .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
344 
345  if (respect_interface) {
346  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
347  "not implemented, set last parameter in refine_TET to false");
348  }
349 
350  Range tets = _tets.subset_by_type(MBTET);
351 
352  std::vector<EntityHandle> parent_ents_refined_and_created;
353  std::vector<EntityHandle> vertices_of_created_tets;
354  std::vector<BitRefEdges> parent_edges_bit_vec;
355  std::vector<int> nb_new_tets_vec;
356  std::vector<int> sub_type_vec;
357 
358  parent_ents_refined_and_created.reserve(tets.size());
359  vertices_of_created_tets.reserve(4 * tets.size());
360  parent_edges_bit_vec.reserve(tets.size());
361  nb_new_tets_vec.reserve(tets.size());
362  sub_type_vec.reserve(tets.size());
363 
364  // make loop over all tets which going to be refined
365  for (auto tit : tets) {
366 
367  if (ref_finite_element.find(tit) == ref_finite_element.end()) {
368  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
369  "this tet is not in refinedFiniteElements");
370  }
371 
372  // get tet connectivity
373  const EntityHandle *conn;
374  int num_nodes;
375  CHKERR moab.get_connectivity(tit, conn, num_nodes, true);
376 
377  // get edges
378  BitRefEdges parent_edges_bit(0);
379  EntityHandle edge_new_nodes[] = {0, 0, 0, 0, 0, 0};
380  int split_edges[] = {-1, -1, -1, -1, -1, -1};
381 
382  for (int ee = 0; ee != 6; ++ee) {
383  EntityHandle edge;
384  CHKERR moab.side_element(tit, 1, ee, edge);
385  RefEntity_multiIndex_view_by_hashed_parent_entity::iterator eit =
386  ref_parent_ents_view.find(edge);
387  if (eit != ref_parent_ents_view.end()) {
388  if (((*eit)->getBitRefLevel() & bit).any()) {
389  edge_new_nodes[ee] = (*eit)->getRefEnt();
390  {
391  const EntityHandle *conn_edge;
392  int num_nodes;
393  moab.get_connectivity(edge, conn_edge, num_nodes, true);
394  if (conn_edge[0] == edge_new_nodes[ee]) {
395  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
396  "node 0 on the edges is mid node, that make no sense");
397  }
398  if (conn_edge[1] == edge_new_nodes[ee]) {
399  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
400  "node 1 on the edges is mid node, that make no sense");
401  }
402  }
403  split_edges[parent_edges_bit.count()] = ee;
404  parent_edges_bit.set(ee, 1);
405  }
406  }
407  }
408 
409  // test if nodes used to refine are not part of tet
410  if (debug) {
411  for (int ee = 0; ee != 6; ee++) {
412  if (edge_new_nodes[ee] == no_handle)
413  continue;
414  for (int nn = 0; nn != 4; nn++) {
415  if (conn[nn] == edge_new_nodes[ee]) {
416  std::cerr << "problem on edge: " << ee << endl;
417  std::cerr << "tet conn : " << conn[0] << " " << conn[1] << " "
418  << conn[2] << " " << conn[3] << " : "
419  << edge_new_nodes[ee] << std::endl;
420  for (int eee = 0; eee != 6; ++eee) {
421  EntityHandle edge;
422  CHKERR moab.side_element(tit, 1, eee, edge);
423  const EntityHandle *conn_edge;
424  int num_nodes;
425  CHKERR moab.get_connectivity(edge, conn_edge, num_nodes, true);
426  std::cerr << eee << " : " << conn_edge[0] << " " << conn_edge[1]
427  << std::endl;
428  }
429  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
430  "nodes used to refine are not part of tet");
431  }
432  }
433  }
434  }
435 
436  // swap nodes forward
437  EntityHandle _conn_[4];
438  std::copy(&conn[0], &conn[4], &_conn_[0]);
439 
440  // build connectivity for ref tets
441  EntityHandle new_tets_conns[8 * 4];
442  std::fill(&new_tets_conns[0], &new_tets_conns[8 * 4], no_handle);
443 
444  int sub_type = -1, nb_new_tets = 0;
445  switch (parent_edges_bit.count()) {
446  case 0: {
447  RefEntity_multiIndex::iterator tit_miit;
448  tit_miit = refined_ents_ptr->find(tit);
449  if (tit_miit == refined_ents_ptr->end())
450  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
451  "can not find this tet");
452  ents_to_set_bit.insert(tit);
453  continue;
454  } break;
455  case 1:
456  sub_type = 0;
457  tet_type_1(_conn_, split_edges[0], edge_new_nodes[split_edges[0]],
458  new_tets_conns);
459  nb_new_tets = 2;
460  break;
461  case 2:
462  sub_type =
463  tet_type_2(_conn_, split_edges, edge_new_nodes, new_tets_conns);
464  if (sub_type & (4 | 8 | 16)) {
465  nb_new_tets = 3;
466  break;
467  } else if (sub_type == 1) {
468  nb_new_tets = 4;
469  break;
470  };
471  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Imposible case");
472  break;
473  case 3:
474  sub_type =
475  tet_type_3(_conn_, split_edges, edge_new_nodes, new_tets_conns);
476  if (sub_type <= 4) {
477  nb_new_tets = 4;
478  break;
479  } else if (sub_type <= 7) {
480  nb_new_tets = 5;
481  break;
482  }
483  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Imposible case");
484  case 4:
485  sub_type =
486  tet_type_4(_conn_, split_edges, edge_new_nodes, new_tets_conns);
487  if (sub_type == 0) {
488  nb_new_tets = 5;
489  break;
490  } else if (sub_type <= 7) {
491  nb_new_tets = 6;
492  break;
493  }
494  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Imposible case");
495  case 5:
496  sub_type = tet_type_5(moab, _conn_, edge_new_nodes, new_tets_conns);
497  nb_new_tets = 7;
498  break;
499  case 6:
500  sub_type = 0;
501  tet_type_6(moab, _conn_, edge_new_nodes, new_tets_conns);
502  nb_new_tets = 8;
503  break;
504  default:
505  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Imposible case");
506  }
507 
508  // find that tets
509  auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.equal_range(
510  boost::make_tuple(tit, parent_edges_bit.to_ulong()));
511  // check if tet with this refinement scheme already exits
512  std::vector<EntityHandle> ents_to_set_bit_vec;
513  if (std::distance(it_by_ref_edges.first, it_by_ref_edges.second) ==
514  static_cast<size_t>(nb_new_tets)) {
515  ents_to_set_bit_vec.reserve(nb_new_tets);
516  for (int tt = 0; it_by_ref_edges.first != it_by_ref_edges.second;
517  it_by_ref_edges.first++, tt++) {
518  auto tet = (*it_by_ref_edges.first)->getRefEnt();
519  ents_to_set_bit_vec.emplace_back(tet);
520  // set ref tets entities
521  if (debug) {
522  // add this tet if exist to this ref level
523  RefEntity_multiIndex::iterator ref_tet_it;
524  ref_tet_it = refined_ents_ptr->find(tet);
525  if (ref_tet_it == refined_ents_ptr->end()) {
526  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
527  "Tetrahedron should be in database");
528  }
529  }
530  }
531  ents_to_set_bit.insert_list(ents_to_set_bit_vec.begin(),
532  ents_to_set_bit_vec.end());
533 
534  } else {
535  // if this element was not refined or was refined with different patterns
536  // of split edges create new elements
537 
538  parent_ents_refined_and_created.emplace_back(tit);
539  for (int tt = 0; tt != nb_new_tets; ++tt) {
540  for (auto nn : {0, 1, 2, 3})
541  vertices_of_created_tets.emplace_back(new_tets_conns[4 * tt + nn]);
542  }
543  parent_edges_bit_vec.emplace_back(parent_edges_bit);
544  nb_new_tets_vec.emplace_back(nb_new_tets);
545  sub_type_vec.emplace_back(sub_type);
546  }
547  }
548 
549  // Create tets
550  EntityHandle start_e = 0;
551  EntityHandle *conn_moab;
552  const int nb_new_tets = vertices_of_created_tets.size() / 4;
553  read_util->get_element_connect(nb_new_tets, 4, MBTET, 0, start_e, conn_moab);
554  std::copy(vertices_of_created_tets.begin(), vertices_of_created_tets.end(),
555  conn_moab);
556  CHKERR read_util->update_adjacencies(start_e, nb_new_tets, 4, conn_moab);
557  ents_to_set_bit.insert(start_e, start_e + nb_new_tets - 1);
558 
559  // Create adj entities
560  for (auto d : {1, 2}) {
561  Range ents;
562  CHKERR moab.get_adjacencies(ents_to_set_bit, d, true, ents,
563  moab::Interface::UNION);
564  }
565 
566  // Set parrents and adjacencies
567  for (int idx = 0; idx != parent_ents_refined_and_created.size(); ++idx) {
568 
569  const EntityHandle tit = parent_ents_refined_and_created[idx];
570  const BitRefEdges &parent_edges_bit = parent_edges_bit_vec[idx];
571  const int nb_new_tets = nb_new_tets_vec[idx];
572  const int sub_type = sub_type_vec[idx];
573 
574  std::array<EntityHandle, 8> ref_tets;
575  for (int tt = 0; tt != nb_new_tets; ++tt, ++start_e)
576  ref_tets[tt] = start_e;
577 
578  if (nb_new_tets) {
579 
580  int ref_type[2];
581  ref_type[0] = parent_edges_bit.count();
582  ref_type[1] = sub_type;
583  for (int tt = 0; tt != nb_new_tets; ++tt) {
584  CHKERR moab.tag_set_data(cOre.get_th_RefType(), &ref_tets[tt], 1,
585  ref_type);
586  CHKERR moab.tag_set_data(cOre.get_th_RefBitEdge(), &ref_tets[tt], 1,
587  &parent_edges_bit);
588  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &ref_tets[tt],
589  1, &tit);
590  }
591 
592  // hash map of nodes (RefEntity) by edges (EntityHandle)
593  std::map<EntityHandle /*edge*/, EntityHandle /*node*/>
594  map_ref_nodes_by_edges;
595 
596  Range tet_edges;
597  CHKERR moab.get_adjacencies(&tit, 1, 1, false, tet_edges);
598  for (auto edge : tet_edges) {
599  RefEntity_multiIndex_view_by_hashed_parent_entity::iterator eit =
600  ref_parent_ents_view.find(edge);
601  if (eit != ref_parent_ents_view.end()) {
602  if (((*eit)->getBitRefLevel() & bit).any()) {
603  map_ref_nodes_by_edges[(*eit)->getParentEnt()] =
604  eit->get()->getRefEnt();
605  }
606  }
607  }
608 
609  // find parents for new edges and faces
610  // get tet edges and faces
611  Range tit_edges, tit_faces;
612  CHKERR moab.get_adjacencies(&tit, 1, 1, false, tit_edges);
613  CHKERR moab.get_adjacencies(&tit, 1, 2, false, tit_faces);
614  Range edges_nodes[6], faces_nodes[4];
615  // for edges - add ref nodes
616  // edges_nodes[ee] - contains all nodes on edge ee including mid nodes if
617  // exist
618  Range::iterator eit = tit_edges.begin();
619  for (int ee = 0; eit != tit_edges.end(); eit++, ee++) {
620  CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee], true);
621  std::map<EntityHandle, EntityHandle>::iterator map_miit =
622  map_ref_nodes_by_edges.find(*eit);
623  if (map_miit != map_ref_nodes_by_edges.end()) {
624  edges_nodes[ee].insert(map_miit->second);
625  }
626  }
627  // for faces - add ref nodes
628  // faces_nodes[ff] - contains all nodes on face ff including mid nodes if
629  // exist
630  Range::iterator fit = tit_faces.begin();
631  for (int ff = 0; fit != tit_faces.end(); fit++, ff++) {
632  CHKERR moab.get_connectivity(&*fit, 1, faces_nodes[ff], true);
633  // Get edges on face and loop over those edges to add mid-nodes to range
634  Range fit_edges;
635  CHKERR moab.get_adjacencies(&*fit, 1, 1, false, fit_edges);
636  for (Range::iterator eit2 = fit_edges.begin(); eit2 != fit_edges.end();
637  eit2++) {
638  std::map<EntityHandle, EntityHandle>::iterator map_miit =
639  map_ref_nodes_by_edges.find(*eit2);
640  if (map_miit != map_ref_nodes_by_edges.end()) {
641  faces_nodes[ff].insert(map_miit->second);
642  }
643  }
644  }
645  // add ref nodes to tet
646  // tet_nodes contains all nodes on tet including mid edge nodes
647  Range tet_nodes;
648  CHKERR moab.get_connectivity(&tit, 1, tet_nodes, true);
649  for (std::map<EntityHandle, EntityHandle>::iterator map_miit =
650  map_ref_nodes_by_edges.begin();
651  map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
652  tet_nodes.insert(map_miit->second);
653  }
654  Range ref_edges;
655  // Get all all edges of refined tets
656  CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 1, false,
657  ref_edges, moab::Interface::UNION);
658  // Check for all ref edge and set parents
659  for (Range::iterator reit = ref_edges.begin(); reit != ref_edges.end();
660  reit++) {
661  Range ref_edges_nodes;
662  CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes, true);
663  if (ref_edges_nodes.size() != 2) {
664  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
665  "data inconsistency, edge should have 2 nodes");
666  }
667  // Check if ref edge is an coarse edge (loop over coarse tet edges)
668  int ee = 0;
669  for (; ee < 6; ee++) {
670  // Two nodes are common (node[0],node[1],ref_node (if exist))
671  // this tests if given edge is contained by edge of refined
672  // tetrahedral
673  if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
674  EntityHandle edge = tit_edges[ee];
675  CHKERR set_parent(*reit, edge, refined_ents_ptr, cOre);
676  break;
677  }
678  }
679  if (ee < 6)
680  continue; // this refined edge is contained by edge of tetrahedral
681  // check if ref edge is in coarse face
682  int ff = 0;
683  for (; ff < 4; ff++) {
684  // two nodes are common (node[0],node[1],ref_node (if exist))
685  // this tests if given face is contained by face of tetrahedral
686  if (intersect(faces_nodes[ff], ref_edges_nodes).size() == 2) {
687  EntityHandle face = tit_faces[ff];
688  CHKERR set_parent(*reit, face, refined_ents_ptr, cOre);
689  break;
690  }
691  }
692  if (ff < 4)
693  continue; // this refined edge is contained by face of tetrahedral
694 
695  // check if ref edge is in coarse tetrahedral (i.e. that is internal
696  // edge of refined tetrahedral)
697  if (intersect(tet_nodes, ref_edges_nodes).size() == 2) {
698  CHKERR set_parent(*reit, tit, refined_ents_ptr, cOre);
699  continue;
700  }
701 
702  // Refined edge is not child of any edge, face or tetrahedral, this is
703  // imposible edge
704  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
705  "impossible refined edge");
706  }
707 
708  Range ref_faces;
709  CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 2, false,
710  ref_faces, moab::Interface::UNION);
711  Tag th_interface_side;
712  const int def_side[] = {0};
713  CHKERR moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
714  th_interface_side,
715  MB_TAG_CREAT | MB_TAG_SPARSE, def_side);
716  // Check for all ref faces
717  for (auto rfit : ref_faces) {
718  Range ref_faces_nodes;
719  CHKERR moab.get_connectivity(&rfit, 1, ref_faces_nodes, true);
720  // Check if ref face is in coarse face
721  int ff = 0;
722  for (; ff < 4; ff++) {
723  // Check if refined triangle is contained by face of tetrahedral
724  if (intersect(faces_nodes[ff], ref_faces_nodes).size() == 3) {
725  EntityHandle face = tit_faces[ff];
726  CHKERR set_parent(rfit, face, refined_ents_ptr, cOre);
727  int side = 0;
728  // Set face side if it is on interface
729  CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
730  CHKERR moab.tag_set_data(th_interface_side, &rfit, 1, &side);
731  break;
732  }
733  }
734  if (ff < 4)
735  continue; // this face is contained by one of tetrahedrons
736  // check if ref face is in coarse tetrahedral
737  // this is ref face which is contained by tetrahedral volume
738  if (intersect(tet_nodes, ref_faces_nodes).size() == 3) {
739  CHKERR set_parent(rfit, tit, refined_ents_ptr, cOre);
740  continue;
741  }
742  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
743  "impossible refined face");
744  }
745  }
746  }
747 
748  if (debug)
749  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
750  CHKERR set_parent(refined_ents_ptr);
751  if (debug)
752  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
753  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(ents_to_set_bit,
754  bit, true, verb);
755  if (debug)
756  CHKERR check(ref_edges_ptr, refined_ents_ptr, cOre);
757 
759 }
Deprecated interface functions.
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
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:483
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:514
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
std::bitset< BITREFEDGES_SIZE > BitRefEdges
Definition: Types.hpp:44
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
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
#define CHKERR
Inline error check.
Definition: definitions.h:602
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
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
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: