v0.8.20
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_verices_in_the_middel_of_edges (const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
 make vertices in the middle of edges in meshset and add them to refinement levels defined by bit More...
 
virtual MoFEMErrorCode add_verices_in_the_middel_of_edges (const Range &edges, const BitRefLevel &bit, int verb=QUIET, EntityHandle start_v=0)
 make vertices in the middle of edges in meshset and add them to Refinement levels defined by bit More...
 
virtual MoFEMErrorCode refine_TET (const EntityHandle meshset, const BitRefLevel &bit, const bool respect_interface=false, int verb=QUIET, Range *ref_edges=NULL)
 refine TET in the meshset More...
 
virtual MoFEMErrorCode refine_TET (const Range &tets, const BitRefLevel &bit, const bool respect_interface=false, int verb=QUIET, Range *ref_edges=NULL)
 refine TET in the meshset More...
 
virtual MoFEMErrorCode refine_PRISM (const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET)
 refine PRISM in the meshset More...
 
virtual MoFEMErrorCode refine_MESHSET (const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET)
 refine meshset, i.e. add child of refined entities to meshset More...
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (const MOFEMuuid &uuid, bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE , bool VERIFY = false>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, IFACE *&iface) const
 Get interface by uuid and return reference to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()
 
virtual MoFEMErrorCode getLibVersion (Version &version) const
 Get library version. More...
 
virtual const MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version) const
 Get database major version. More...
 
virtual MoFEMErrorCode getInterfaceVersion (Version &version) const
 Get database major version. More...
 
template<>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, UnknownInterface *&iface) const
 

Public Attributes

MoFEM::CorecOre
 

Additional Inherited Members

- Protected Member Functions inherited from MoFEM::UnknownInterface
boost::typeindex::type_index getClassIdx (const MOFEMuuid &uid) const
 Get type name for interface Id. More...
 
MOFEMuuid getUId (const boost::typeindex::type_index &class_idx) const
 Get interface Id for class name. More...
 

Detailed Description

Mesh refinement interface.

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

Bug:

Not working on partitioned meshes

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

If outsourced, class member functions should follow name convention

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

Definition at line 41 of file MeshRefinement.hpp.

Constructor & Destructor Documentation

◆ MeshRefinement()

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

Definition at line 34 of file MeshRefinement.cpp.

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

◆ ~MeshRefinement()

MoFEM::MeshRefinement::~MeshRefinement ( )

destructor

Definition at line 50 of file MeshRefinement.hpp.

50 {}

Member Function Documentation

◆ add_verices_in_the_middel_of_edges() [1/2]

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 
)
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_verices_in_the_middel_of_edges(edges, bit, verb, start_v);
90 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
virtual 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)
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]

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

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

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

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

Definition at line 91 of file MeshRefinement.cpp.

93  {
94  Interface &m_field = cOre;
95  moab::Interface &moab = m_field.get_moab();
96  const RefEntity_multiIndex *refined_ents_ptr;
98  CHKERR m_field.get_ref_ents(&refined_ents_ptr);
99  Range edges = _edges.subset_by_type(MBEDGE);
100  typedef const RefEntity_multiIndex::index<
101  Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
102  RefEntsByComposite &ref_ents =
103  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
104  RefEntsByComposite::iterator miit =
105  ref_ents.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE));
106  RefEntsByComposite::iterator hi_miit =
107  ref_ents.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE));
109  ref_parent_ents_view.insert(miit,hi_miit);
110  if (verb >= VERBOSE) {
111  std::ostringstream ss;
112  ss << "ref level " << bit << " nb. edges to refine " << edges.size()
113  << std::endl;
114  PetscPrintf(m_field.get_comm(), ss.str().c_str());
115  }
116  std::vector<double> vert_coords[3];
117  for (int dd = 0; dd != 3; dd++) {
118  vert_coords[dd].reserve(edges.size());
119  }
120  std::vector<EntityHandle> parent_edge;
121  parent_edge.reserve(edges.size());
122  Range add_bit;
123  for (Range::iterator eit = edges.begin(); eit != edges.end(); ++eit) {
124  RefEntity_multiIndex_view_by_hashed_parent_entity::iterator miit_view =
125  ref_parent_ents_view.find(*eit);
126  if (miit_view == ref_parent_ents_view.end()) {
127  const EntityHandle *conn;
128  int num_nodes;
129  CHKERR moab.get_connectivity(*eit, conn, num_nodes, true);
130  if (num_nodes != 2) {
131  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
132  "edge should have 2 edges");
133  }
134  double coords[6];
135  CHKERR moab.get_coords(conn, num_nodes, coords);
136  cblas_daxpy(3, 1., &coords[3], 1, coords, 1);
137  cblas_dscal(3, 0.5, coords, 1);
138  for(int dd = 0;dd!=3;++dd)
139  vert_coords[dd].push_back(coords[dd]);
140  parent_edge.push_back(*eit);
141  } else {
142  if ((*miit_view)->getEntType() != MBVERTEX) {
143  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
144  "child of edge should be vertex");
145  }
146  const EntityHandle node = (*miit_view)->getRefEnt();
147  RefEntity_multiIndex::iterator nit = refined_ents_ptr->find(node);
148  if(nit->get()->getParentEnt() != *eit) {
149  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
150  "parent edge does not much");
151  }
152  add_bit.insert(node);
153  }
154  }
155  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit, bit);
156  if (!vert_coords[0].empty()) {
157  int num_nodes = vert_coords[0].size();
158  vector<double *> arrays_coord;
159  EntityHandle startv;
160  {
161  ReadUtilIface *iface;
162  CHKERR moab.query_interface(iface);
163  CHKERR iface->get_node_coords(3, num_nodes, start_v, startv,
164  arrays_coord);
165  }
166  Range verts(startv, startv + num_nodes - 1);
167  for (int dd = 0; dd != 3; ++dd) {
168  std::copy(vert_coords[dd].begin(), vert_coords[dd].end(),
169  arrays_coord[dd]);
170  }
171  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), verts,
172  &*parent_edge.begin());
173  CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
174  verts, bit, verb);
175  }
177 }
Tag get_th_RefParentHandle() const
Definition: Core.hpp:150
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
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
void cblas_dscal(const int N, const double alpha, double *X, const int incX)
Definition: cblas_dscal.c:11
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< hashed_non_unique< const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< boost::shared_ptr< RefEntity >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getRefEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_hashed_parent_entity
multi-index view of RefEntity by parent entity
#define CHKERR
Inline error check.
Definition: definitions.h:594
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ 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 }
MeshRefinement(const MoFEM::Core &core)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506
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 903 of file MeshRefinement.cpp.

905  {
906  Interface &m_field = cOre;
907  const RefEntity_multiIndex *refined_ents_ptr;
909  CHKERR m_field.get_ref_ents(&refined_ents_ptr);
910  typedef const RefEntity_multiIndex::index<Ent_mi_tag>::type RefEntsByEnt;
911  RefEntsByEnt::iterator miit = refined_ents_ptr->find(meshset);
912  if (miit == refined_ents_ptr->end()) {
913  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
914  "this meshset is not in ref database");
915  }
916  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
917  meshset, bit, meshset, MBEDGE, recursive, verb);
918  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
919  meshset, bit, meshset, MBTRI, recursive, verb);
920  CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
921  meshset, bit, meshset, MBTET, recursive, verb);
922  *(const_cast<RefEntity *>(miit->get())->getBitRefLevelPtr()) |= bit;
924 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

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

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

◆ refine_TET() [1/2]

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

refine TET in the meshset

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

Definition at line 179 of file MeshRefinement.cpp.

182  {
183  Interface &m_field = cOre;
184  moab::Interface &moab = m_field.get_moab();
186  Range tets;
187  CHKERR moab.get_entities_by_type(meshset, MBTET, tets, false);
188  CHKERR refine_TET(tets, bit, respect_interface, verb, ref_edges_ptr);
190 }
virtual MoFEMErrorCode refine_TET(const EntityHandle meshset, const BitRefLevel &bit, const bool respect_interface=false, int verb=QUIET, Range *ref_edges=NULL)
refine TET in the meshset
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ refine_TET() [2/2]

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

refine TET in the meshset

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

Definition at line 192 of file MeshRefinement.cpp.

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