v0.8.13
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=-1, 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=0, 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=0, 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=0, Range *ref_edges=NULL)
 refine TET in the meshset More...
 
virtual MoFEMErrorCode refine_PRISM (const EntityHandle meshset, const BitRefLevel &bit, int verb=0)
 refine PRISM in the meshset More...
 
virtual MoFEMErrorCode refine_MESHSET (const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=0)
 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 = -1,
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:495
#define CHKERR
Inline error check.
Definition: definitions.h:614
virtual MoFEMErrorCode add_verices_in_the_middel_of_edges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=-1, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit ...
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

◆ 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 = 0,
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:495
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:614
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:439

◆ 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:519
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:526
static const MOFEMuuid IDD_MOFEMMeshRefine

◆ refine_MESHSET()

MoFEMErrorCode MoFEM::MeshRefinement::refine_MESHSET ( const EntityHandle  meshset,
const BitRefLevel bit,
const bool  recursive = false,
int  verb = 0 
)
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 902 of file MeshRefinement.cpp.

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

◆ refine_PRISM()

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

refine PRISM in the meshset

Parameters
EntityHandlemeshset
BitRefLevelbitLevel

Definition at line 694 of file MeshRefinement.cpp.

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

refine TET in the meshset

Parameters
EntityHandlemeshset
BitRefLevelbitLevel
IfTRUE, interface elements would be refined too

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 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
virtual MoFEMErrorCode refine_TET(const EntityHandle meshset, const BitRefLevel &bit, const bool respect_interface=false, int verb=0, Range *ref_edges=NULL)
refine TET in the meshset
#define CHKERR
Inline error check.
Definition: definitions.h:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

◆ refine_TET() [2/2]

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

refine TET in the meshset

Parameters
Rangeof tets to refine
BitRefLevelbitLevel
IfTRUE, interface elements would be refined too

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