v0.6.9
Public Member Functions | Public Attributes | List of all members
MoFEM::PrismInterface Struct Reference

Make interface on given faces and create flat prism in that space. More...

#include <src/interfaces/PrismInterface.hpp>

Inheritance diagram for MoFEM::PrismInterface:
[legend]
Collaboration diagram for MoFEM::PrismInterface:
[legend]

Public Member Functions

MoFEMErrorCode query_interface (const MOFEMuuid &uuid, UnknownInterface **iface) const
 
 PrismInterface (const MoFEM::Core &core)
 
 ~PrismInterface ()
 destructor More...
 
MoFEMErrorCode getSides (const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
 create two children meshsets in the meshset containing tetrahedral on two sides of faces More...
 
MoFEMErrorCode getSides (const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
 create two children meshsets in the meshset containing tetrahedral on two sides of faces More...
 
MoFEMErrorCode findIfTringleHasThreeNodesOnInternalSurfaceSkin (const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, Range &faces_with_three_nodes_on_front, int verb=QUIET)
 Find if triangle has three nodes on internal surface skin. More...
 
MoFEMErrorCode splitSides (const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
 split nodes and other entities of tetrahedral in children sets and add prism elements More...
 
MoFEMErrorCode splitSides (const EntityHandle meshset, const BitRefLevel &bit, const EntityHandle, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
 split nodes and other entities of tetrahedral in children sets and add prism elements More...
 
MoFEMErrorCode splitSides (const EntityHandle meshset, const BitRefLevel &bit, const BitRefLevel &inhered_from_bit_level, const BitRefLevel &inhered_from_bit_level_mask, const EntityHandle sideset, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
 split nodes and other entities of tetrahedrons in children sets and add prism elements 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>
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 >
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type >
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

Make interface on given faces and create flat prism in that space.

Todo:
FIXME Names of methods do not follow naming convention and are difficult to work with.
Examples:
split_sideset.cpp.

Definition at line 37 of file PrismInterface.hpp.

Constructor & Destructor Documentation

◆ PrismInterface()

MoFEM::PrismInterface::PrismInterface ( const MoFEM::Core core)

Definition at line 36 of file PrismInterface.cpp.

36  :
37 cOre(const_cast<Core&>(core)) {
38 }

◆ ~PrismInterface()

MoFEM::PrismInterface::~PrismInterface ( )

destructor

Definition at line 46 of file PrismInterface.hpp.

46 {}

Member Function Documentation

◆ findIfTringleHasThreeNodesOnInternalSurfaceSkin()

MoFEMErrorCode MoFEM::PrismInterface::findIfTringleHasThreeNodesOnInternalSurfaceSkin ( const EntityHandle  sideset,
const BitRefLevel  mesh_bit_level,
const bool  recursive,
Range &  faces_with_three_nodes_on_front,
int  verb = QUIET 
)

Find if triangle has three nodes on internal surface skin.

Internal surface skin is a set of edges in interface of the body on boundary of surface. This set of edges is called surface front. If surface face has three nodes on surface front, non of the face nodes is split and should be removed from surface if it is going to be split.

Parameters
sidesetmeshset with surface
mesh_bit_levelbit ref level of the volume mesh
recursivesearch in sub-meshsets
faces_with_three_nodes_on_frontreturned faces
verbverbosity level
Returns
error code

Definition at line 447 of file PrismInterface.cpp.

449  {
450  Interface &m_field = cOre;
451  moab::Interface &moab = m_field.get_moab();
453 
454  Range mesh_level_ents3d;
455  Range mesh_level_edges,mesh_level_tris;
456  if(mesh_bit_level.any()) {
457  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
458  mesh_bit_level, BitRefLevel().set(), MBTET, mesh_level_ents3d);
459  CHKERRG(ierr);
460  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
461  mesh_bit_level, BitRefLevel().set(), MBTRI, mesh_level_tris);
462  CHKERRG(ierr);
463  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
464  mesh_bit_level, BitRefLevel().set(), MBEDGE, mesh_level_edges);
465  CHKERRG(ierr);
466  }
467 
468  Skinner skin(&moab);
469  //get interface triangles from side set
470  Range triangles;
471  rval = moab.get_entities_by_type(sideset, MBTRI, triangles, recursive);
473  if(mesh_bit_level.any()) {
474  triangles = intersect(triangles,mesh_level_tris);
475  }
476  if(verb>=VERBOSE) {
477  PetscPrintf(m_field.get_comm(),"Nb. of triangles in set %u\n",triangles.size());
478  }
479  //get nodes, edges and 3d ents (i.e. tets and prisms)
480  Range nodes; // nodes from triangles
481  rval = moab.get_connectivity(triangles, nodes, true);
483  Range ents3d; // 3d ents form nodes
484  rval = moab.get_adjacencies(nodes, 3, false, ents3d, moab::Interface::UNION);
486  if(mesh_bit_level.any()) {
487  ents3d = intersect(ents3d,mesh_level_ents3d);
488  }
489  //take skin faces
490  Range skin_faces; // skin faces from 3d ents
491  rval = skin.find_skin(0, ents3d.subset_by_type(MBTET), false, skin_faces);
493  //take skin edges (boundary of surface if there is any)
494  Range skin_edges_boundary; //skin edges from triangles
495  rval = skin.find_skin(0, triangles, false, skin_edges_boundary);
497  if(verb>=VERY_VERBOSE) {
498  PetscPrintf(m_field.get_comm(), "skin_edges_boundary %u\n",
499  skin_edges_boundary.size());
500  }
501  //take all edges on skin faces (i.e. skin surface)
502  Range skin_faces_edges; //edges from skin faces of 3d ents
503  rval = moab.get_adjacencies(skin_faces, 1, false, skin_faces_edges,
504  moab::Interface::UNION);
506  if(mesh_bit_level.any()) {
507  skin_faces_edges = intersect(skin_faces_edges,mesh_level_edges);
508  }
509  if(verb>=VERY_VERBOSE) {
510  PetscPrintf(m_field.get_comm(),"skin_faces_edges %u\n",skin_faces_edges.size());
511  }
512  //note: that skin faces edges do not contain internal boundary
513  //note: that prisms are not included in ents3d, so if ents3d have border with other inteface is like external boundary
514  //skin edges boundary are internal edge <- skin_faces_edges contains edges which are on the body boundary <- that is the trick
515  skin_edges_boundary = subtract(skin_edges_boundary,skin_faces_edges); // from skin edges subtract edges from skin faces of 3d ents (only internal edges)
516 
517  if(verb>=VERY_VERBOSE) {
518  EntityHandle out_meshset;
519  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
520  rval = moab.add_entities(out_meshset,triangles); CHKERRQ_MOAB(rval);
521  rval = moab.write_file("triangles.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
522  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
523  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
524  rval = moab.add_entities(out_meshset,ents3d); CHKERRQ_MOAB(rval);
525  rval = moab.write_file("ents3d.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
526  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
527  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
528  rval = moab.add_entities(out_meshset,skin_edges_boundary); CHKERRQ_MOAB(rval);
529  rval = moab.write_file("skin_edges_boundary.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
530  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
531  }
532  if(verb>=VERY_VERBOSE) {
533  PetscPrintf(
534  m_field.get_comm(),"subtract skin_edges_boundary %u\n",skin_edges_boundary.size()
535  );
536  }
537 
538  //Get nodes on boundary edge
539  Range skin_nodes_boundary;
540  rval = moab.get_connectivity(skin_edges_boundary, skin_nodes_boundary, true);
542  //Remove node which are boundary with other existing interface
543  Range prisms_nodes;
544  rval =
545  moab.get_connectivity(ents3d.subset_by_type(MBPRISM), prisms_nodes, true);
547  skin_nodes_boundary = subtract(skin_nodes_boundary,prisms_nodes);
548  if(verb>=VERY_VERBOSE) {
549  PetscPrintf(m_field.get_comm(), "subtract skin_nodes_boundary %u\n",
550  skin_nodes_boundary.size());
551  }
552  //use nodes on body boundary and interface (without internal boundary) to find adjacent tets
553  Range nodes_without_front = subtract(nodes,skin_nodes_boundary); // nodes_without_front adjacent to all split face edges except those on internal edge
554  if(verb>=VERY_VERBOSE) {
555  PetscPrintf(m_field.get_comm(),
556  "adj. node if ents3d but not on the internal edge %u\n",
557  nodes_without_front.size());
558  }
559 
560  Range skin_nodes_boundary_tris;
561  rval = moab.get_adjacencies(
562  skin_nodes_boundary,2,false,skin_nodes_boundary_tris,moab::Interface::UNION
563  ); CHKERRG(ierr);
564  Range skin_nodes_boundary_tris_nodes;
565  rval = moab.get_connectivity(skin_nodes_boundary_tris,
566  skin_nodes_boundary_tris_nodes, true);
568  skin_nodes_boundary_tris_nodes =
569  subtract(skin_nodes_boundary_tris_nodes, skin_nodes_boundary);
570  Range skin_nodes_boundary_tris_nodes_tris;
571  rval = moab.get_adjacencies(skin_nodes_boundary_tris_nodes, 2, false,
572  skin_nodes_boundary_tris_nodes_tris,
573  moab::Interface::UNION);
574  CHKERRG(ierr);
575  // Triangle which has tree nodes on front boundary
576  skin_nodes_boundary_tris = intersect(
577  triangles,
578  subtract(skin_nodes_boundary_tris,skin_nodes_boundary_tris_nodes_tris)
579  );
580  faces_with_three_nodes_on_front.swap(skin_nodes_boundary_tris);
581 
583 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:482
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78

◆ getSides() [1/2]

MoFEMErrorCode MoFEM::PrismInterface::getSides ( const int  msId,
const CubitBCType  cubit_bc_type,
const BitRefLevel  mesh_bit_level,
const bool  recursive,
int  verb = QUIET 
)

create two children meshsets in the meshset containing tetrahedral on two sides of faces

Parameters
msIdId of meshset
cubit_bc_typetype of meshset (NODESET, SIDESET or BLOCKSET and more)
mesh_bit_leveladd interface on bit level is bit_level = BitRefLevel.set() then add interface on all bit levels
recursiveif true parent meshset is searched recursively
Examples:
split_sideset.cpp.

Definition at line 40 of file PrismInterface.cpp.

43  {
44 
45  Interface &m_field = cOre;
46  MeshsetsManager *meshsets_manager_ptr;
48  ierr = m_field.getInterface(meshsets_manager_ptr);
49  CHKERRG(ierr);
50  CubitMeshSet_multiIndex::index<
51  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator miit =
52  meshsets_manager_ptr->getMeshsetsMultindex()
53  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
54  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
55  if (miit !=
56  meshsets_manager_ptr->getMeshsetsMultindex()
57  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
58  .end()) {
59  ierr = getSides(miit->meshset, mesh_bit_level, recursive, verb);
60  CHKERRG(ierr);
61  } else {
62  SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_FOUND,"msId is not there");
63  }
65 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
create two children meshsets in the meshset containing tetrahedral on two sides of faces ...

◆ getSides() [2/2]

MoFEMErrorCode MoFEM::PrismInterface::getSides ( const EntityHandle  sideset,
const BitRefLevel  mesh_bit_level,
const bool  recursive,
int  verb = QUIET 
)

create two children meshsets in the meshset containing tetrahedral on two sides of faces

Get tets adj to faces. Take skin form tets and get edges from that skin. Take skin form triangles (the face). Subtrac skin faces edges form skin edges in order to get edges on the boundary of the face which is in the volume of the body, but is not on the boundary. Each child set has a child containing nodes which can be split and skin edges. After that simply iterate under all tets on one side which are adjacent to the face are found. Side tets are stored in to children meshsets of the SIDESET meshset.

Definition at line 67 of file PrismInterface.cpp.

69  {
70  Interface &m_field = cOre;
71  moab::Interface &moab = m_field.get_moab();
73  Range mesh_level_ents3d,mesh_level_ents3d_tris;
74  Range mesh_level_tris;
75  Range mesh_level_edges;
76  Range mesh_level_nodes;
77  if (mesh_bit_level.any()) {
78  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
79  mesh_bit_level, BitRefLevel().set(), MBTET, mesh_level_ents3d);
80  CHKERRG(ierr);
81  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
82  mesh_bit_level, BitRefLevel().set(), MBTRI, mesh_level_tris);
83  CHKERRG(ierr);
84  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
85  mesh_bit_level, BitRefLevel().set(), MBEDGE, mesh_level_edges);
86  CHKERRG(ierr);
87  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
88  mesh_bit_level, BitRefLevel().set(), MBVERTEX, mesh_level_nodes);
89  CHKERRG(ierr);
90  rval = moab.get_adjacencies(mesh_level_ents3d, 2, false,
91  mesh_level_ents3d_tris, moab::Interface::UNION);
93  }
94  Range mesh_level_prisms;
95  if(mesh_bit_level.any()) {
96  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
97  mesh_bit_level, BitRefLevel().set(), MBPRISM, mesh_level_prisms);
98  CHKERRG(ierr);
99  mesh_level_ents3d.merge(mesh_level_prisms);
100  }
101  Skinner skin(&moab);
102  //get interface triangles from side set
103  Range triangles;
104  rval = moab.get_entities_by_type(sideset, MBTRI, triangles, recursive);
106  if(mesh_bit_level.any()) {
107  triangles = intersect(triangles,mesh_level_ents3d_tris);
108  }
109  if(verb>=VERBOSE) {
110  PetscPrintf(m_field.get_comm(), "Nb. of triangles in set %u\n",
111  triangles.size());
112  }
113  //get nodes, edges and 3d ents (i.e. tets and prisms)
114  Range nodes; // nodes from triangles
115  rval = moab.get_connectivity(triangles,nodes,true); CHKERRQ_MOAB(rval);
116  Range ents3d,ents3d_with_prisms; // 3d ents form nodes
117  rval = moab.get_adjacencies(nodes, 3, false, ents3d_with_prisms,
118  moab::Interface::UNION);
120  if(mesh_bit_level.any()) {
121  ents3d_with_prisms = intersect(ents3d_with_prisms,mesh_level_ents3d);
122  }
123  ents3d = ents3d_with_prisms.subset_by_type(MBTET); // take only tets, add prism later
124  //take skin faces
125  Range skin_faces; // skin faces from 3d ents
126  rval = skin.find_skin(0,ents3d,false,skin_faces); CHKERRQ_MOAB(rval);
127  //take skin edges (boundary of surface if there is any)
128  Range skin_edges_boundary; //skin edges from triangles
129  rval = skin.find_skin(0, triangles, false, skin_edges_boundary);
131  if (verb >= VERY_VERBOSE)
132  PetscPrintf(m_field.get_comm(), "skin_edges_boundary %u\n",
133  skin_edges_boundary.size());
134  //take all edges on skin faces (i.e. skin surface)
135  Range skin_faces_edges; //edges from skin faces of 3d ents
136  rval = moab.get_adjacencies(skin_faces, 1, false, skin_faces_edges,
137  moab::Interface::UNION);
139  if(mesh_bit_level.any()) {
140  skin_faces_edges = intersect(skin_faces_edges,mesh_level_edges);
141  }
142  if (verb >= VERY_VERBOSE)
143  PetscPrintf(m_field.get_comm(), "skin_faces_edges %u\n",
144  skin_faces_edges.size());
145  //note: that skin faces edges do not contain internal boundary
146  //note: that prisms are not included in ents3d, so if ents3d have border with other inteface is like external boundary
147  //skin edges boundary are internal edge <- skin_faces_edges contains edges which are on the body boundary <- that is the trick
148  skin_edges_boundary = subtract(skin_edges_boundary,skin_faces_edges); // from skin edges subtract edges from skin faces of 3d ents (only internal edges)
149  if (verb >= NOISY) {
150  EntityHandle out_meshset;
151  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, out_meshset);
153  rval = moab.add_entities(out_meshset, triangles);
155  rval = moab.write_file("triangles.vtk", "VTK", "", &out_meshset, 1);
157  rval = moab.delete_entities(&out_meshset, 1);
159  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, out_meshset);
161  rval = moab.add_entities(out_meshset, ents3d);
163  rval = moab.write_file("ents3d.vtk", "VTK", "", &out_meshset, 1);
165  rval = moab.delete_entities(&out_meshset, 1);
167  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, out_meshset);
169  rval = moab.add_entities(out_meshset, skin_edges_boundary);
171  rval =
172  moab.write_file("skin_edges_boundary.vtk", "VTK", "", &out_meshset, 1);
174  rval = moab.delete_entities(&out_meshset, 1);
176  }
177  if (verb >= VERY_VERBOSE)
178  PetscPrintf(m_field.get_comm(), "subtract skin_edges_boundary %u\n",
179  skin_edges_boundary.size());
180  //Get nodes on boundary edge
181  Range skin_nodes_boundary;
182  rval = moab.get_connectivity(skin_edges_boundary, skin_nodes_boundary, true);
184  //Remove node which are boundary with other existing interface
185  Range prisms_nodes;
186  rval = moab.get_connectivity(ents3d_with_prisms.subset_by_type(MBPRISM),
187  prisms_nodes, true);
189  skin_nodes_boundary = subtract(skin_nodes_boundary,prisms_nodes);
190  if (verb >= VERY_VERBOSE)
191  PetscPrintf(m_field.get_comm(), "subtract skin_nodes_boundary %u\n",
192  skin_nodes_boundary.size());
193  //use nodes on body boundary and interface (without internal boundary) to find adjacent tets
194  Range nodes_without_front = subtract(nodes,skin_nodes_boundary); // nodes_without_front adjacent to all split face edges except those on internal edge
195  if (verb >= VERY_VERBOSE)
196  PetscPrintf(m_field.get_comm(),
197  "adj. node if ents3d but not on the internal edge %u\n",
198  nodes_without_front.size());
199  //ents3 that are adjacent to front nodes on split faces but not those which are on the front nodes on internal edge
200  ents3d.clear();
201  ents3d_with_prisms.clear();
202  rval = moab.get_adjacencies(
203  nodes_without_front,3,false,ents3d_with_prisms,moab::Interface::UNION
204  ); CHKERRQ_MOAB(rval);
205 
206  // Add tets adjacent to front and triangle which has all nodes on crack front
207  struct FindTrianglesOnFrontAndAdjacentTet {
208  static MoFEMErrorCode fUN(moab::Interface &moab, const Range &triangles,
209  const Range &skin_nodes_boundary, ///< nodes on front
210  const Range &skin_edges_boundary, ///< edges on front
211  Range &ents3d_with_prisms) {
213  // get all triangles adjacent to front
214  Range skin_nodes_boundary_tris;
215  rval = moab.get_adjacencies(skin_nodes_boundary, 2, false,
216  skin_nodes_boundary_tris,
217  moab::Interface::UNION);
218  CHKERRG(ierr);
219  // get nodes of triangles adjacent to front nodes
220  Range skin_nodes_boundary_tris_nodes;
221  rval = moab.get_connectivity(skin_nodes_boundary_tris,
222  skin_nodes_boundary_tris_nodes, true);
224  // get hanging nodes, i.e. nodes which are not on the front but adjacent
225  // to triangles adjacent to crack front
226  skin_nodes_boundary_tris_nodes =
227  subtract(skin_nodes_boundary_tris_nodes, skin_nodes_boundary);
228  // get triangles adjacent to hanging nodes
229  Range skin_nodes_boundary_tris_nodes_tris;
230  rval = moab.get_adjacencies(skin_nodes_boundary_tris_nodes, 2, false,
231  skin_nodes_boundary_tris_nodes_tris,
232  moab::Interface::UNION);
233  CHKERRG(ierr);
234  // triangles which have tree nodes on front boundary
235  skin_nodes_boundary_tris =
236  intersect(triangles, subtract(skin_nodes_boundary_tris,
237  skin_nodes_boundary_tris_nodes_tris));
238  if (!skin_nodes_boundary_tris.empty()) {
239  // Get internal edges of triangle which has three nodes on boundary
240  Range skin_nodes_boundary_tris_edges;
241  rval = moab.get_adjacencies(skin_nodes_boundary_tris, 1, false,
242  skin_nodes_boundary_tris_edges,
243  moab::Interface::UNION);
244  CHKERRG(rval);
245  skin_nodes_boundary_tris_edges =
246  subtract(skin_nodes_boundary_tris_edges, skin_edges_boundary);
247  // Get 3d elements adjacent to internal edge which has two nodes on
248  // boundary
249  rval = moab.get_adjacencies(skin_nodes_boundary_tris_edges, 3, false,
250  ents3d_with_prisms, moab::Interface::UNION);
251  CHKERRG(rval);
252  }
254  }
255  };
256  ierr = FindTrianglesOnFrontAndAdjacentTet::fUN(
257  moab,triangles,skin_nodes_boundary,skin_edges_boundary,ents3d_with_prisms
258  ); CHKERRG(ierr);
259 
260  // prism and tets on both side of interface
261  if(mesh_bit_level.any()) {
262  ents3d_with_prisms = intersect(ents3d_with_prisms,mesh_level_ents3d);
263  }
264  ents3d = ents3d_with_prisms.subset_by_type(MBTET);
265  if (verb >= VERY_VERBOSE)
266  PetscPrintf(m_field.get_comm(), "adj. ents3d to front nodes %u\n",
267  ents3d.size());
268 
269  Range side_ents3d;
270  unsigned int nb_side_ents3d = side_ents3d.size();
271  side_ents3d.insert(*ents3d.begin());
272  Range side_ents3d_tris_on_surface;
273 
274  // get all tets adjacent to crack surface, but only on one side of it
275  do {
276 
277  do {
278  Range adj_tris,adj_ents3d;
279  nb_side_ents3d = side_ents3d.size();
280  if (verb >= VERBOSE)
281  PetscPrintf(m_field.get_comm(), "nb_side_ents3d %u\n", nb_side_ents3d);
282  //get faces
283  rval = moab.get_adjacencies(side_ents3d.subset_by_type(MBTET), 2, false,
284  adj_tris, moab::Interface::UNION);
285  CHKERRQ_MOAB(rval);
286  if(mesh_bit_level.any()) {
287  adj_tris = intersect(adj_tris,mesh_level_tris);
288  }
289  //subtrace from faces interface
290  adj_tris = subtract(adj_tris,triangles);
291  if (verb >= VERBOSE)
292  PetscPrintf(m_field.get_comm(), "adj_tris %u\n", adj_tris.size());
293  //get tets adjacent to faces
294  rval = moab.get_adjacencies(
295  adj_tris,3,true,adj_ents3d,moab::Interface::UNION
296  ); CHKERRQ_MOAB(rval);
297  //intersect tets with tets adjacent to inetface
298  adj_ents3d = intersect(adj_ents3d,ents3d_with_prisms);
299  if (verb >= VERBOSE)
300  PetscPrintf(m_field.get_comm(), "adj_ents3d %u\n", adj_ents3d.size());
301  //add tets to side
302  side_ents3d.insert(adj_ents3d.begin(),adj_ents3d.end());
303  if (verb >= VERY_NOISY) {
304  EntityHandle out_meshset;
305  rval =
306  moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, out_meshset);
307  CHKERRQ_MOAB(rval);
308  rval = moab.add_entities(out_meshset,side_ents3d);
309  CHKERRQ_MOAB(rval);
310  std::string file = "side_ents3d_" +
311  boost::lexical_cast<std::string>(nb_side_ents3d) +
312  ".vtk";
313  rval = moab.write_file(file.c_str(), "VTK", "", &out_meshset, 1);
314  CHKERRQ_MOAB(rval);
315  rval = moab.delete_entities(&out_meshset, 1);
316  CHKERRQ_MOAB(rval);
317  }
318  } while (nb_side_ents3d != side_ents3d.size());
319 
320  Range side_ents3d_tris;
321  rval = moab.get_adjacencies(
322  side_ents3d,2,false,side_ents3d_tris,moab::Interface::UNION
323  ); CHKERRQ_MOAB(rval);
324  side_ents3d_tris_on_surface = intersect(side_ents3d_tris,triangles);
325 
326  if (verb >= VERY_NOISY) {
327  Range left_triangles = subtract(triangles, side_ents3d_tris_on_surface);
328  if(!left_triangles.empty()) {
329  EntityHandle out_meshset;
330  rval =
331  moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, out_meshset);
332  CHKERRQ_MOAB(rval);
333  rval = moab.add_entities(out_meshset, left_triangles);
334  CHKERRQ_MOAB(rval);
335  std::string file = "left_triangles_" +
336  boost::lexical_cast<std::string>(nb_side_ents3d) +
337  ".vtk";
338  rval = moab.write_file(file.c_str(), "VTK", "", &out_meshset, 1);
339  CHKERRQ_MOAB(rval);
340  rval = moab.delete_entities(&out_meshset, 1);
341  CHKERRQ_MOAB(rval);
342  }
343  }
344 
345  // This is a case when separate sub-domains are split, so wee need
346  // additional
347  // tetrahedron for seed process
348  if (side_ents3d_tris_on_surface.size() != triangles.size()) {
349  Range left_triangles = subtract(triangles, side_ents3d_tris_on_surface);
350  Range tets;
351  rval = moab.get_adjacencies(&*left_triangles.begin(), 1, 3, false, tets);
352  CHKERRQ_MOAB(rval);
353  tets = intersect(tets, ents3d_with_prisms);
354  if (tets.empty()) {
355  Range left_triangles_nodes;
356  rval = moab.get_connectivity(&*left_triangles.begin(), 1,
357  left_triangles_nodes, true);
358  CHKERRG(ierr);
359  EntityHandle meshset;
360  rval = moab.create_meshset(MESHSET_SET, meshset);
361  CHKERRQ_MOAB(rval);
362  rval = moab.add_entities(meshset, left_triangles);
363  CHKERRQ_MOAB(rval);
364  rval = moab.write_file("error.vtk", "VTK", "", &meshset, 1);
365  CHKERRQ_MOAB(rval);
366  rval = moab.delete_entities(&meshset, 1);
367  CHKERRQ_MOAB(rval);
368  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
369  "Not all faces on surface going to be split, see error.vtk for "
370  "problematic triangle. "
371  "It could be a case where triangle on front (part boundary of "
372  "surface in interior) "
373  "has three nodes front.");
374  }
375  side_ents3d.insert(*tets.begin());
376  }
377 
378  } while (side_ents3d_tris_on_surface.size()!=triangles.size());
379 
380  if(ents3d_with_prisms.size() == side_ents3d.size()) {
381  SETERRQ(
382  m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,
383  "All tets on one side, no-interface"
384  );
385  }
386  //other side ents
387  Range other_side = subtract(ents3d_with_prisms,side_ents3d);
388  //side nodes
389  Range side_nodes;
390  rval = moab.get_connectivity(side_ents3d.subset_by_type(MBTET), side_nodes,
391  true);
392  CHKERRQ_MOAB(rval);
393  //nodes on crack surface without front
394  nodes_without_front = intersect(nodes_without_front,side_nodes);
395  Range side_edges;
396  rval = moab.get_adjacencies(side_ents3d.subset_by_type(MBTET), 1, false,
397  side_edges, moab::Interface::UNION);
398  CHKERRQ_MOAB(rval);
399  skin_edges_boundary = intersect(skin_edges_boundary,side_edges);
400  //make child meshsets
401  std::vector<EntityHandle> children;
402  rval = moab.get_child_meshsets(sideset,children); CHKERRQ_MOAB(rval);
403  if(children.empty()) {
404  children.resize(3);
405  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, children[0]);
406  CHKERRQ_MOAB(rval);
407  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, children[1]);
408  CHKERRQ_MOAB(rval);
409  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, children[2]);
410  CHKERRQ_MOAB(rval);
411  rval = moab.add_child_meshset(sideset,children[0]); CHKERRQ_MOAB(rval);
412  rval = moab.add_child_meshset(sideset,children[1]); CHKERRQ_MOAB(rval);
413  rval = moab.add_child_meshset(sideset,children[2]); CHKERRQ_MOAB(rval);
414  } else {
415  if(children.size()!=3) {
416  SETERRQ(PETSC_COMM_SELF,1,"this meshset should have 3 children meshsets");
417  }
418  children.resize(3);
419  ierr = moab.clear_meshset(&children[0],3); CHKERRG(ierr);
420  }
421  EntityHandle &child_side = children[0];
422  EntityHandle &child_other_side = children[1];
423  EntityHandle &child_nodes_and_skin_edges = children[2];
424  rval = moab.add_entities(child_side,side_ents3d); CHKERRQ_MOAB(rval);
425  rval = moab.add_entities(child_other_side,other_side); CHKERRQ_MOAB(rval);
426  rval = moab.add_entities(child_nodes_and_skin_edges, nodes_without_front);
427  CHKERRQ_MOAB(rval);
428  rval = moab.add_entities(child_nodes_and_skin_edges, skin_edges_boundary);
429  CHKERRQ_MOAB(rval);
430  if(verb>=VERBOSE) {
431  PetscPrintf(m_field.get_comm(), "Nb. of side ents3d in set %u\n",
432  side_ents3d.size());
433  PetscPrintf(m_field.get_comm(), "Nb. of other side ents3d in set %u\n",
434  other_side.size());
435  PetscPrintf(m_field.get_comm(), "Nb. of boundary edges %u\n",
436  skin_edges_boundary.size());
437  }
438  if(verb>=NOISY) {
439  ierr = moab.write_file("side.vtk", "VTK", "", &children[0], 1);
440  CHKERRG(ierr);
441  ierr = moab.write_file("other_side.vtk", "VTK", "", &children[1], 1);
442  CHKERRG(ierr);
443  }
445 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:482
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Common.hpp:60
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78

◆ query_interface()

MoFEMErrorCode MoFEM::PrismInterface::query_interface ( const MOFEMuuid uuid,
UnknownInterface **  iface 
) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 24 of file PrismInterface.cpp.

25  {
27  *iface = NULL;
28  if(uuid == IDD_MOFEMPrismInterface) {
29  *iface = const_cast<PrismInterface*>(this);
31  }
32  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"unknown interface");
34 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
PrismInterface(const MoFEM::Core &core)
static const MOFEMuuid IDD_MOFEMPrismInterface

◆ splitSides() [1/3]

MoFEMErrorCode MoFEM::PrismInterface::splitSides ( const EntityHandle  meshset,
const BitRefLevel bit,
const int  msId,
const CubitBCType  cubit_bc_type,
const bool  add_interface_entities,
const bool  recursive = false,
int  verb = QUIET 
)

split nodes and other entities of tetrahedral in children sets and add prism elements

The all new entities (prisms, tets) are added to refinement level given by bit

Parameters
meshsetmeshset to get entities from
BitRefLevelnew level where refinement would be stored
msIdmeshset ID imported from cubit
cubit_bc_typetype of meshset (NODESET, SIDESET or BLOCKSET and more)
add_interface_entitiesmeshset which contain the interface
recursiveif true parent meshset is searched recursively

Each inteface face has two tags, const int def_side[] = {0}; rval = moab.tag_get_handle("INTERFACE_SIDE",1,MB_TYPE_INTEGER, th_interface_side,MB_TAG_CREAT|MB_TAG_SPARSE,def_side); CHKERRQ_MOAB(rval);

const EntityHandle def_node[] = {0}; rval = moab.tag_get_handle("SIDE_INTFACE_ELEMENT",1,MB_TYPE_HANDLE, th_side_elem,MB_TAG_CREAT|MB_TAG_SPARSE,def_node); CHKERRQ_MOAB(rval);

First tag inform about inteface side, second tag inform about side adjacent inteface element.

Examples:
split_sideset.cpp.

Definition at line 585 of file PrismInterface.cpp.

590  {
591 
592  Interface &m_field = cOre;
593  MeshsetsManager *meshsets_manager_ptr;
595  ierr = m_field.getInterface(meshsets_manager_ptr); CHKERRG(ierr);
596  CubitMeshSet_multiIndex::index<
597  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator miit =
598  meshsets_manager_ptr->getMeshsetsMultindex()
599  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
600  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
601  if (miit !=
602  meshsets_manager_ptr->getMeshsetsMultindex()
603  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
604  .end()) {
605  ierr = splitSides(meshset, bit, miit->meshset, add_interface_entities,
606  recursive, verb);
607  CHKERRG(ierr);
608  } else {
609  SETERRQ(PETSC_COMM_SELF,1,"msId is not there");
610  }
612 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
split nodes and other entities of tetrahedral in children sets and add prism elements ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80

◆ splitSides() [2/3]

MoFEMErrorCode MoFEM::PrismInterface::splitSides ( const EntityHandle  meshset,
const BitRefLevel bit,
const EntityHandle  sideset,
const bool  add_interface_entities,
const bool  recursive = false,
int  verb = QUIET 
)

split nodes and other entities of tetrahedral in children sets and add prism elements

The all new entities (prisms, tets) are added to refinement level given by bit

Definition at line 613 of file PrismInterface.cpp.

617  {
618 
620  ierr = splitSides(meshset, bit, BitRefLevel(), BitRefLevel(), sideset,
621  add_interface_entities, recursive, verb);
622  CHKERRG(ierr);
624 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
split nodes and other entities of tetrahedral in children sets and add prism elements ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80

◆ splitSides() [3/3]

MoFEMErrorCode MoFEM::PrismInterface::splitSides ( const EntityHandle  meshset,
const BitRefLevel bit,
const BitRefLevel inhered_from_bit_level,
const BitRefLevel inhered_from_bit_level_mask,
const EntityHandle  sideset,
const bool  add_interface_entities,
const bool  recursive = false,
int  verb = QUIET 
)

split nodes and other entities of tetrahedrons in children sets and add prism elements

The all new entities (prisms, tets) are added to refinement level given by bit

Parameters
meshset
Refinementbit level of new mesh
inhered_from_bit_levelinhered nodes and other entities form this bit level.
add_interface_entitiesadd prism elements at interface
recursivedo meshesets in the meshset

note inhered_from_bit_level is need to be specified to some meshsets with interfaces. Some nodes on some refinement levels dividing edges but not splitting faces. Inheriting those nodes will not split faces.

Definition at line 625 of file PrismInterface.cpp.

629  {
630 
631  Interface &m_field = cOre;
632  moab::Interface &moab = m_field.get_moab();
634  std::vector<EntityHandle> children;
635  //get children meshsets
636  rval = moab.get_child_meshsets(sideset,children); CHKERRQ_MOAB(rval);
637  if(children.size()!=3) {
638  SETERRQ(
639  m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,
640  "should be 3 child meshsets, each of them contains tets on two sides of interface"
641  );
642  }
643  //3d ents on "father" side
644  Range side_ents3d;
645  rval = moab.get_entities_by_handle(children[0], side_ents3d, false);
647  //3d ents on "mather" side
648  Range other_ents3d;
649  rval = moab.get_entities_by_handle(children[1], other_ents3d, false);
651  //faces of interface
652  Range triangles;
653  rval = moab.get_entities_by_type(sideset, MBTRI, triangles, recursive);
655  Range side_ents3d_tris;
656  rval = moab.get_adjacencies(side_ents3d, 2, true, side_ents3d_tris,
657  moab::Interface::UNION);
659  triangles = intersect(triangles,side_ents3d_tris);
660  //nodes on interface but not on crack front (those should not be splitted)
661  Range nodes;
662  rval = moab.get_entities_by_type(children[2], MBVERTEX, nodes, false);
664  Range meshset_3d_ents,meshset_2d_ents;
665  rval = moab.get_entities_by_dimension(meshset, 3, meshset_3d_ents, true);
667  Range meshset_tets = meshset_3d_ents.subset_by_type(MBTET);
668  rval = moab.get_adjacencies(meshset_tets, 2, false, meshset_2d_ents,
669  moab::Interface::UNION);
671  side_ents3d = intersect(meshset_3d_ents,side_ents3d);
672  other_ents3d = intersect(meshset_3d_ents,other_ents3d);
673  triangles = intersect(meshset_2d_ents,triangles);
674  if(verb>3) {
675  PetscPrintf(m_field.get_comm(),"triangles %u\n",triangles.size());
676  PetscPrintf(m_field.get_comm(),"side_ents3d %u\n",side_ents3d.size());
677  PetscPrintf(m_field.get_comm(),"nodes %u\n",nodes.size());
678  }
679 
680  const RefEntity_multiIndex *refined_ents_ptr;
681  ierr = m_field.get_ref_ents(&refined_ents_ptr); CHKERRG(ierr);
682  typedef RefEntity_multiIndex::index<Ent_mi_tag>::type RefEntsByEntType;
683  const RefEntsByEntType &ref_ents_by_type = refined_ents_ptr->get<Ent_mi_tag>();
684  RefEntity_multiIndex_view_by_parent_entity ref_parent_ents_view;
685  //create view index by parent entity
686  {
687  typedef RefEntity_multiIndex::index<
688  Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
689  const RefEntsByComposite &ref_ents =
690  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
691 
692  RefEntsByComposite::iterator miit;
693  RefEntsByComposite::iterator hi_miit;
694  //view by parent type (VERTEX)
695  miit = ref_ents.lower_bound(boost::make_tuple(MBVERTEX,MBVERTEX));
696  hi_miit = ref_ents.upper_bound(boost::make_tuple(MBVERTEX,MBVERTEX));
697  for(;miit!=hi_miit;miit++) {
698  if(((*miit)->getBitRefLevel()&inhered_from_bit_level_mask) == (*miit)->getBitRefLevel()) {
699  if(((*miit)->getBitRefLevel()&inhered_from_bit_level).any()) {
700  std::pair<RefEntity_multiIndex_view_by_parent_entity::iterator,bool> p_ref_ent_view;
701  p_ref_ent_view = ref_parent_ents_view.insert(*miit);
702  if(!p_ref_ent_view.second) {
703  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"non unique insertion");
704  }
705  }
706  }
707  }
708  }
709  //maps nodes on "father" and "mather" side
710  std::map<
711  EntityHandle, /*node on "mather" side*/
712  EntityHandle /*node on "father" side*/
713  > map_nodes;
714  //add new nodes on interface and create map
715  Range::iterator nit = nodes.begin();
716  double coord[3];
717  for(;nit!=nodes.end();nit++) {
718  //find ref enet
719  RefEntsByEntType::iterator miit_ref_ent = ref_ents_by_type.find(*nit);
720  if(miit_ref_ent == ref_ents_by_type.end()) {
721  SETERRQ(PETSC_COMM_SELF,1,"can not find node in MoFEM database");
722  }
723  EntityHandle child_entity = 0;
724  RefEntity_multiIndex::iterator child_it;
725  RefEntity_multiIndex_view_by_parent_entity::iterator child_iit;
726  child_iit = ref_parent_ents_view.find(*nit);
727  if(child_iit != ref_parent_ents_view.end()) {
728  child_it = refined_ents_ptr->find((*child_iit)->getRefEnt());
729  BitRefLevel bit_child = (*child_it)->getBitRefLevel();
730  if( (inhered_from_bit_level&bit_child).none() ) {
731  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"data inconsistency");
732  }
733  child_entity = (*child_it)->getRefEnt();
734  }
735  //
736  bool success;
737  if(child_entity == 0) {
738  rval = moab.get_coords(&*nit,1,coord); CHKERRQ_MOAB(rval);
739  EntityHandle new_node;
740  rval = moab.create_vertex(coord,new_node); CHKERRQ_MOAB(rval);
741  map_nodes[*nit] = new_node;
742  //create new node on "father" side
743  //parent is node on "mather" side
744  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(),&new_node,1,&*nit); CHKERRQ_MOAB(rval);
745  std::pair<RefEntity_multiIndex::iterator,bool> p_ref_ent =
746  const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->insert(
747  boost::shared_ptr<RefEntity>
748  (new RefEntity(m_field.get_basic_entity_data_ptr(),new_node))
749  );
750  //set ref bit level to node on "father" side
751  success = const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->
752  modify(p_ref_ent.first,RefEntity_change_add_bit(bit));
753  if (!success)
754  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
755  "modification unsuccessful");
756  } else {
757  map_nodes[*nit] = child_entity;
758  //set ref bit level to node on "father" side
759  success = const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->
760  modify(child_it,RefEntity_change_add_bit(bit));
761  if (!success)
762  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
763  "modification unsuccessful");
764  }
765  //set ref bit level to node on "mather" side
766  success = const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->
767  modify(miit_ref_ent,RefEntity_change_add_bit(bit));
768  if (!success)
769  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
770  "modification unsuccessful");
771  }
772  //crete meshset for new mesh bit level
773  EntityHandle meshset_for_bit_level;
774  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
775  meshset_for_bit_level);
777  //subtract those elements which will be refined, i.e. disconnected form other side elements, and connected to new prisms, if they area created
778  meshset_3d_ents = subtract(meshset_3d_ents,side_ents3d);
779  rval = moab.add_entities(meshset_for_bit_level, meshset_3d_ents);
781  for(int dd = 0;dd<3;dd++) {
782  Range ents_dd;
783  rval = moab.get_adjacencies(meshset_3d_ents, dd, false, ents_dd,
784  moab::Interface::UNION);
786  rval = moab.add_entities(meshset_for_bit_level, ents_dd);
788  }
789  //
790  //typedef RefEntity_multiIndex::index<Composite_ParentEnt_And_EntType_mi_tag>::type ref_ent_by_composite;
791  //ref_ent_by_composite &by_composite = refinedEntities.get<Composite_ParentEnt_And_EntType_mi_tag>();
792  //create new 3d ents on "father" side
793  Range new_3d_ents;
794  for (Range::iterator eit3d = side_ents3d.begin(); eit3d != side_ents3d.end();
795  eit3d++) {
796  RefEntsByEntType::iterator miit_ref_ent = ref_ents_by_type.find(*eit3d);
797  if(miit_ref_ent==ref_ents_by_type.end()) {
798  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
799  "tet not in database");
800  }
801  int num_nodes;
802  const EntityHandle* conn;
803  rval = moab.get_connectivity(*eit3d, conn, num_nodes, true);
805  EntityHandle new_conn[num_nodes];
806  int nb_new_conn = 0;
807  for (int ii = 0; ii < num_nodes; ii++) {
808  std::map<EntityHandle, EntityHandle>::iterator mit =
809  map_nodes.find(conn[ii]);
810  if (mit != map_nodes.end()) {
811  new_conn[ii] = mit->second;
812  nb_new_conn++;
813  if (verb >= VERY_NOISY) {
814  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
815  new_conn[ii]);
816  }
817  } else {
818  new_conn[ii] = conn[ii];
819  }
820  }
821  if (nb_new_conn == 0) {
822  // Add this tet to bit ref level
823  rval = moab.add_entities(meshset_for_bit_level, &*eit3d, 1);
824  rval = moab.add_entities(meshset_for_bit_level, conn, num_nodes);
825  continue;
826  }
827 
828  const RefElement_multiIndex *refined_finite_elements_ptr;
829  ierr = m_field.get_ref_finite_elements(&refined_finite_elements_ptr);
830  CHKERRG(ierr);
831 
832  //here is created new or prism is on interface
833  EntityHandle existing_ent = 0;
834  /* check if tet element with new connectivity is in database*/
835  RefElement_multiIndex::index<Ent_Ent_mi_tag>::type::iterator child_iit,hi_child_iit;
836  child_iit =
837  refined_finite_elements_ptr->get<Ent_Ent_mi_tag>().lower_bound(*eit3d);
838  hi_child_iit =
839  refined_finite_elements_ptr->get<Ent_Ent_mi_tag>().upper_bound(*eit3d);
840  for(;child_iit!=hi_child_iit;child_iit++) {
841  const EntityHandle* conn_ref_tet;
842  rval = moab.get_connectivity(child_iit->getRefEnt(), conn_ref_tet,
843  num_nodes, true);
845  int nn = 0;
846  for(;nn<num_nodes;nn++) {
847  if(conn_ref_tet[nn]!=new_conn[nn]) {
848  break;
849  }
850  }
851  if(nn == num_nodes) {
852  if(existing_ent != 0) {
853  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
854  "database inconsistency");
855  }
856  existing_ent = child_iit->getRefEnt();
857  }
858  }
859  switch (moab.type_from_handle(*eit3d)) {
860  case MBTET: {
861  RefEntsByEntType::iterator child_it;
862  EntityHandle tet;
863  if(existing_ent == 0) {
864  Range new_conn_tet;
865  rval = moab.get_adjacencies(new_conn, 4, 3, false, new_conn_tet);
867  if(new_conn_tet.empty()) {
868  rval = moab.create_element(MBTET, new_conn, 4, tet);
870  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(), &tet, 1,
871  &*eit3d);
873  } else {
874  RefElement_multiIndex::index<Ent_mi_tag>::type::iterator rit,
875  new_rit;
876  rit = refined_finite_elements_ptr->get<Ent_mi_tag>().find(*eit3d);
877  if(rit==refined_finite_elements_ptr->get<Ent_mi_tag>().end()) {
878  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
879  "can't find this in database");
880  }
881  new_rit = refined_finite_elements_ptr->get<Ent_mi_tag>().find(
882  *new_conn_tet.begin());
883  if(new_rit==refined_finite_elements_ptr->get<Ent_mi_tag>().end()) {
884  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
885  "can't find this in database");
886  }
887  tet = *new_conn_tet.begin();
888  /*std::ostringstream ss;
889  ss << "nb new conns: " << nb_new_conn << std::endl;
890  ss << "new_conn_tets.size() " << new_conn_tet.size() << std::endl;
891  ss << "data inconsistency\n";
892  ss << "this ent:\n";
893  ss << *rit->ref_ptr << std::endl;
894  ss << "found this ent:\n";
895  ss << *new_rit->ref_ptr << std::endl;
896  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());*/
897  }
898  } else {
899  tet = existing_ent;
900  }
901  rval = moab.add_entities(meshset_for_bit_level, &tet, 1);
903  rval = moab.add_entities(meshset_for_bit_level, new_conn, 4);
905  new_3d_ents.insert(tet);
906  } break;
907  case MBPRISM: {
908  EntityHandle prism;
909  if(verb>3) {
910  PetscPrintf(m_field.get_comm(), "prims nb_new_nodes %d\n",
911  nb_new_conn);
912  }
913  if(existing_ent == 0) {
914  Range new_conn_prism;
915  rval = moab.get_adjacencies(new_conn, 6, 3, false, new_conn_prism);
917  if(new_conn_prism.empty()) {
918  rval = moab.create_element(MBPRISM, new_conn, 6, prism);
920  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(), &prism, 1,
921  &*eit3d);
923  } else {
924  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
925  "data inconsistency");
926  }
927  } else {
928  prism = existing_ent;
929  }
930  rval = moab.add_entities(meshset_for_bit_level, &prism, 1);
932  rval = moab.add_entities(meshset_for_bit_level, new_conn, 4);
934  new_3d_ents.insert(prism);
935  } break;
936  default:
937  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"not implemented");
938  }
939  }
940  Range new_ents;
941  //create new entities by adjacencies form new tets
942  rval = moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET), 2, true,
943  new_ents, moab::Interface::UNION);
945  rval = moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET), 1, true,
946  new_ents, moab::Interface::UNION);
948  //Tags for setting side
949  Tag th_interface_side;
950  const int def_side[] = {0};
951  rval = moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
952  th_interface_side, MB_TAG_CREAT | MB_TAG_SPARSE,
953  def_side);
955  //add new edges and triangles to mofem database
956  Range ents;
957  rval =
958  moab.get_adjacencies(triangles, 1, false, ents, moab::Interface::UNION);
960  ents.insert(triangles.begin(),triangles.end());
961  Range new_ents_in_database; //this range contains all new entities
962  for(Range::iterator eit = ents.begin();eit!=ents.end();eit++) {
963  int num_nodes;
964  const EntityHandle* conn;
965  rval = moab.get_connectivity(*eit, conn, num_nodes, true);
967  EntityHandle new_conn[num_nodes];
968  int nb_new_conn = 0;
969  int ii = 0;
970  for(;ii<num_nodes; ii++) {
971  std::map<EntityHandle, EntityHandle>::iterator mit =
972  map_nodes.find(conn[ii]);
973  if(mit != map_nodes.end()) {
974  new_conn[ii] = mit->second;
975  nb_new_conn++;
976  if (verb >= VERY_NOISY) {
977  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
978  new_conn[ii]);
979  }
980  } else {
981  new_conn[ii] = conn[ii];
982  }
983  }
984  if(nb_new_conn==0) continue;
985  RefEntsByEntType::iterator miit_ref_ent = ref_ents_by_type.find(*eit);
986  if(miit_ref_ent == ref_ents_by_type.end()) {
987  SETERRQ(PETSC_COMM_SELF, 1,
988  "this entity (edge or tri) should be already in database");
989  }
990  Range new_ent; //contains all entities (edges or triangles) added to mofem database
991  switch (moab.type_from_handle(*eit)) {
992  case MBTRI: {
993  //get entity based on its connectivity
994  rval = moab.get_adjacencies(new_conn, 3, 2, false, new_ent);
996  if (new_ent.size() != 1)
997  SETERRQ(PETSC_COMM_SELF, 1, "this tri should be in moab database");
998  int new_side = 1;
999  rval = moab.tag_set_data(th_interface_side, &*new_ent.begin(), 1,
1000  &new_side);
1001  CHKERRQ_MOAB(rval);
1002  if (verb >= VERY_VERBOSE)
1003  PetscPrintf(m_field.get_comm(), "new_ent %u\n", new_ent.size());
1004  //add prism element
1005  if(add_interface_entities) {
1006  if(inhered_from_bit_level.any()) {
1007  SETERRQ(PETSC_COMM_SELF, 1,
1008  "not implemented for inhered_from_bit_level");
1009  }
1010  //set prism connectivity
1011  EntityHandle prism_conn[6] = {
1012  conn[0],conn[1],conn[2],
1013  new_conn[0],new_conn[1],new_conn[2]
1014  };
1015  EntityHandle prism;
1016  rval = moab.create_element(MBPRISM, prism_conn, 6, prism);
1017  CHKERRQ_MOAB(rval);
1018  ierr = cOre.addPrismToDatabase(prism, verb);
1019  CHKERRG(ierr);
1020  rval = moab.add_entities(meshset_for_bit_level, &prism, 1);
1021  CHKERRQ_MOAB(rval);
1022  }
1023  } break;
1024  case MBEDGE: {
1025  rval = moab.get_adjacencies(new_conn, 2, 1, false, new_ent);
1026  CHKERRQ_MOAB(rval);
1027  if (new_ent.size() != 1) {
1028  if (m_field.get_comm_rank() == 0) {
1029  EntityHandle out_meshset;
1030  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
1031  out_meshset);
1032  CHKERRQ_MOAB(rval);
1033  rval = moab.add_entities(out_meshset, &*eit, 1);
1034  CHKERRQ_MOAB(rval);
1035  rval = moab.write_file("debug_splitSides.vtk", "VTK", "",
1036  &out_meshset, 1);
1037  CHKERRQ_MOAB(rval);
1038  rval = moab.delete_entities(&out_meshset, 1);
1039  CHKERRQ_MOAB(rval);
1040  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
1041  out_meshset);
1042  CHKERRQ_MOAB(rval);
1043  rval = moab.add_entities(out_meshset, side_ents3d);
1044  CHKERRQ_MOAB(rval);
1045  rval = moab.write_file("debug_splitSides_side_ents3d.vtk", "VTK",
1046  "", &out_meshset, 1);
1047  CHKERRQ_MOAB(rval);
1048  rval = moab.delete_entities(&out_meshset, 1);
1049  CHKERRQ_MOAB(rval);
1050  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
1051  out_meshset);
1052  CHKERRQ_MOAB(rval);
1053  rval = moab.add_entities(out_meshset, other_ents3d);
1054  CHKERRQ_MOAB(rval);
1055  rval = moab.write_file("debug_splitSides_other_ents3d.vtk", "VTK",
1056  "", &out_meshset, 1);
1057  CHKERRQ_MOAB(rval);
1058  rval = moab.delete_entities(&out_meshset, 1);
1059  CHKERRQ_MOAB(rval);
1060  rval = moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
1061  out_meshset);
1062  CHKERRQ_MOAB(rval);
1063  rval = moab.add_entities(out_meshset, triangles);
1064  CHKERRQ_MOAB(rval);
1065  rval = moab.write_file("debug_get_msId_3dENTS_split_triangles.vtk",
1066  "VTK", "", &out_meshset, 1);
1067  CHKERRQ_MOAB(rval);
1068  rval = moab.delete_entities(&out_meshset, 1);
1069  CHKERRQ_MOAB(rval);
1070  }
1071  SETERRQ2(PETSC_COMM_SELF, 1, "this edge should be in moab database "
1072  "new_ent.size() = %u nb_new_conn = %d",
1073  new_ent.size(), nb_new_conn);
1074  }
1075  } break;
1076  default:
1077  SETERRQ(PETSC_COMM_SELF, 1, "Houston we have a problem !!!");
1078  }
1079  if (new_ent.size() != 1) {
1080  SETERRQ1(PETSC_COMM_SELF, 1,
1081  "new_ent.size() = %u, size always should be 1",
1082  new_ent.size());
1083  }
1084  //set parent
1085  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(), &*new_ent.begin(),
1086  1, &*eit);
1087  CHKERRQ_MOAB(rval);
1088  // add to database
1089  std::pair<RefEntity_multiIndex::iterator, bool> p_ref_ent =
1090  const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
1091  ->insert(boost::shared_ptr<RefEntity>(new RefEntity(
1092  m_field.get_basic_entity_data_ptr(), new_ent[0])));
1093  const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
1094  ->modify(p_ref_ent.first, RefEntity_change_add_bit(bit));
1095  new_ents_in_database.insert(new_ent.begin(), new_ent.end());
1096  }
1097  //all other entities, some ents like triangles and faces on the side of tets
1098  Range side_adj_faces_and_edges;
1099  rval = moab.get_adjacencies(side_ents3d.subset_by_type(MBTET), 1, true,
1100  side_adj_faces_and_edges, moab::Interface::UNION);
1101  CHKERRQ_MOAB(rval);
1102  rval = moab.get_adjacencies(side_ents3d.subset_by_type(MBTET), 2, true,
1103  side_adj_faces_and_edges, moab::Interface::UNION);
1104  CHKERRQ_MOAB(rval);
1105  //subtract entities already added to mofem database
1106  side_adj_faces_and_edges =
1107  subtract(side_adj_faces_and_edges, new_ents_in_database);
1108  for (Range::iterator eit = side_adj_faces_and_edges.begin();
1109  eit != side_adj_faces_and_edges.end(); eit++) {
1110  int num_nodes;
1111  const EntityHandle* conn;
1112  rval = moab.get_connectivity(*eit, conn, num_nodes, true);
1113  CHKERRQ_MOAB(rval);
1114  EntityHandle new_conn[num_nodes];
1115  int nb_new_conn = 0;
1116  int ii = 0;
1117  for (; ii < num_nodes; ii++) {
1118  std::map<EntityHandle, EntityHandle>::iterator mit =
1119  map_nodes.find(conn[ii]);
1120  if (mit != map_nodes.end()) {
1121  new_conn[ii] = mit->second;
1122  nb_new_conn++;
1123  if (verb >= VERY_NOISY) {
1124  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
1125  new_conn[ii]);
1126  }
1127  } else {
1128  new_conn[ii] = conn[ii];
1129  }
1130  }
1131  if (nb_new_conn == 0)
1132  continue;
1133  RefEntsByEntType::iterator miit_ref_ent = ref_ents_by_type.find(*eit);
1134  if (miit_ref_ent == ref_ents_by_type.end()) {
1135  SETERRQ1(PETSC_COMM_SELF, 1,
1136  "entity should be in MoFem database, num_nodes = %d", num_nodes);
1137  }
1138  Range new_ent;
1139  switch (moab.type_from_handle(*eit)) {
1140  case MBTRI: {
1141  rval = moab.get_adjacencies(new_conn, 3, 2, false, new_ent);
1142  CHKERRQ_MOAB(rval);
1143  } break;
1144  case MBEDGE: {
1145  rval = moab.get_adjacencies(new_conn, 2, 1, false, new_ent);
1146  CHKERRQ_MOAB(rval);
1147  } break;
1148  default:
1149  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1150  "Houston we have a problem");
1151  }
1152  if (new_ent.size() != 1) {
1153  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1154  "database inconsistency, new_ent.size() = %u", new_ent.size());
1155  }
1156  // add entity to mofem database
1157  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(), &*new_ent.begin(),
1158  1, &*eit);
1159  CHKERRQ_MOAB(rval);
1160  std::pair<RefEntity_multiIndex::iterator, bool> p_ref_ent =
1161  const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
1162  ->insert(boost::shared_ptr<RefEntity>(new RefEntity(
1163  m_field.get_basic_entity_data_ptr(), new_ent[0])));
1164  const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
1165  ->modify(p_ref_ent.first, RefEntity_change_add_bit(bit));
1166  if (verb > VERY_VERBOSE)
1167  PetscPrintf(m_field.get_comm(), "new_ent %u\n", new_ent.size());
1168  new_ents_in_database.insert(new_ent.begin(), new_ent.end());
1169  }
1170  // add new prisms which parents are part of other intefaces
1171  Range new_3d_prims = new_3d_ents.subset_by_type(MBPRISM);
1172  for (Range::iterator pit = new_3d_prims.begin(); pit != new_3d_prims.end();
1173  pit++) {
1174  ierr = cOre.addPrismToDatabase(*pit, verb);
1175  CHKERRG(ierr);
1176  // get parent entity
1177  EntityHandle parent_prism;
1178  rval = moab.tag_get_data(cOre.get_th_RefParentHandle(), &*pit, 1,
1179  &parent_prism);
1180  CHKERRQ_MOAB(rval);
1181  const EntityHandle root_meshset = moab.get_root_set();
1182  if (parent_prism == root_meshset) {
1183  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1184  "this prism should have parent");
1185  }
1186  if (moab.type_from_handle(parent_prism) != MBPRISM) {
1187  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1188  "this prism should have parent which is prism as well");
1189  }
1190  int num_nodes;
1191  // parent prism
1192  const EntityHandle *conn_parent;
1193  rval = moab.get_connectivity(parent_prism, conn_parent, num_nodes, true);
1194  MOAB_THROW(rval);
1195  Range face_side3_parent, face_side4_parent;
1196  rval = moab.get_adjacencies(conn_parent, 3, 2, false, face_side3_parent);
1197  CHKERRQ_MOAB(rval);
1198  rval =
1199  moab.get_adjacencies(&conn_parent[3], 3, 2, false, face_side4_parent);
1200  CHKERRQ_MOAB(rval);
1201  if (face_side3_parent.size() != 1) {
1202  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1203  "parent face3.size() = %u", face_side3_parent.size());
1204  }
1205  if (face_side4_parent.size() != 1) {
1206  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1207  "parent face4.size() = %u", face_side4_parent.size());
1208  }
1209  // new prism
1210  const EntityHandle *conn;
1211  rval = moab.get_connectivity(*pit, conn, num_nodes, true);
1212  MOAB_THROW(rval);
1213  Range face_side3, face_side4;
1214  rval = moab.get_adjacencies(conn, 3, 2, false, face_side3);
1215  CHKERRQ_MOAB(rval);
1216  rval = moab.get_adjacencies(&conn[3], 3, 2, false, face_side4);
1217  CHKERRQ_MOAB(rval);
1218  if (face_side3.size() != 1) {
1219  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "face3 is missing");
1220  }
1221  if (face_side4.size() != 1) {
1222  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "face4 is missing");
1223  }
1224  //
1225  std::vector<EntityHandle> face(2), parent_face(2);
1226  face[0] = *face_side3.begin();
1227  face[1] = *face_side4.begin();
1228  parent_face[0] = *face_side3_parent.begin();
1229  parent_face[1] = *face_side4_parent.begin();
1230  for (int ff = 0; ff < 2; ff++) {
1231  if (parent_face[ff] == face[ff])
1232  continue;
1233  int interface_side;
1234  rval = moab.tag_get_data(th_interface_side, &parent_face[ff], 1,
1235  &interface_side);
1236  CHKERRQ_MOAB(rval);
1237  rval =
1238  moab.tag_set_data(th_interface_side, &face[ff], 1, &interface_side);
1239  CHKERRQ_MOAB(rval);
1240  EntityHandle parent_tri;
1241  rval = moab.tag_get_data(cOre.get_th_RefParentHandle(), &face[ff], 1,
1242  &parent_tri);
1243  CHKERRQ_MOAB(rval);
1244  if (parent_tri != parent_face[ff]) {
1245  SETERRQ1(PETSC_COMM_SELF, 1, "wrong parent %lu", parent_tri);
1246  }
1247  if (new_ents_in_database.find(face[ff]) == new_ents_in_database.end()) {
1248  RefEntity_multiIndex::index<Ent_mi_tag>::type::iterator miit_ref_ent =
1249  refined_ents_ptr->get<Ent_mi_tag>().find(face[ff]);
1250  if (miit_ref_ent == refined_ents_ptr->get<Ent_mi_tag>().end()) {
1251  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1252  "this is not in database, but should not be");
1253  }
1254  }
1255  }
1256  }
1257  // finalise by adding new tets and prism ti bit level
1258  ierr = m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
1259  meshset_for_bit_level, 3, bit);
1260  CHKERRG(ierr);
1261  rval = moab.delete_entities(&meshset_for_bit_level, 1);
1262  CHKERRQ_MOAB(rval);
1263  ierr = moab.clear_meshset(&children[0], 3);
1264  CHKERRG(ierr);
1266 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:482
MoFEMErrorCode addPrismToDatabase(const EntityHandle prism, int verb=-1)
add prim element
Definition: Core.cpp:233
Tag get_th_RefParentHandle() const
Definition: Core.hpp:138
multi_index_container< ptrWrapperRefElement, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityHandle, &ptrWrapperRefElement::getRefEnt > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityHandle, &ptrWrapperRefElement::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityType, &ptrWrapperRefElement::getEntType > >, ordered_non_unique< tag< Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag >, composite_key< ptrWrapperRefElement, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityHandle, &ptrWrapperRefElement::getParentEnt >, const_mem_fun< ptrWrapperRefElement::interface_type_RefElement, int, &ptrWrapperRefElement::getBitRefEdgesUlong > > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< ptrWrapperRefElement, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityHandle, &ptrWrapperRefElement::getRefEnt >, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityHandle, &ptrWrapperRefElement::getParentEnt > > > >> RefElement_multiIndex
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:468
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< hashed_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_parent_entity
multi-index view of RefEntity by parent entity
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:511
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:474
#define MOAB_THROW(a)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:569
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
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:28
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, EntityHandle,&RefEntity::getParentEnt >, const_mem_fun< RefEntity::BasicEntity, EntityType,&RefEntity::getEntType > > > >> RefEntity_multiIndex

Member Data Documentation

◆ cOre

MoFEM::Core& MoFEM::PrismInterface::cOre

Definition at line 42 of file PrismInterface.hpp.


The documentation for this struct was generated from the following files: