v0.8.14
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 = 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

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 392 of file PrismInterface.cpp.

394  {
395  Interface &m_field = cOre;
396  moab::Interface &moab = m_field.get_moab();
398 
399  Range mesh_level_ents3d;
400  Range mesh_level_edges,mesh_level_tris;
401  if(mesh_bit_level.any()) {
402  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
403  mesh_bit_level, BitRefLevel().set(), MBTET, mesh_level_ents3d);
404  CHKERRG(ierr);
405  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
406  mesh_bit_level, BitRefLevel().set(), MBTRI, mesh_level_tris);
407  CHKERRG(ierr);
408  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
409  mesh_bit_level, BitRefLevel().set(), MBEDGE, mesh_level_edges);
410  CHKERRG(ierr);
411  }
412 
413  Skinner skin(&moab);
414  //get interface triangles from side set
415  Range triangles;
416  rval = moab.get_entities_by_type(sideset, MBTRI, triangles, recursive);
418  if(mesh_bit_level.any()) {
419  triangles = intersect(triangles,mesh_level_tris);
420  }
421  if(verb>=VERBOSE) {
422  PetscPrintf(m_field.get_comm(),"Nb. of triangles in set %u\n",triangles.size());
423  }
424  //get nodes, edges and 3d ents (i.e. tets and prisms)
425  Range nodes; // nodes from triangles
426  rval = moab.get_connectivity(triangles, nodes, true);
428  Range ents3d; // 3d ents form nodes
429  rval = moab.get_adjacencies(nodes, 3, false, ents3d, moab::Interface::UNION);
431  if(mesh_bit_level.any()) {
432  ents3d = intersect(ents3d,mesh_level_ents3d);
433  }
434  //take skin faces
435  Range skin_faces; // skin faces from 3d ents
436  rval = skin.find_skin(0, ents3d.subset_by_type(MBTET), false, skin_faces);
438  //take skin edges (boundary of surface if there is any)
439  Range skin_edges_boundary; //skin edges from triangles
440  rval = skin.find_skin(0, triangles, false, skin_edges_boundary);
442  if(verb>=VERY_VERBOSE) {
443  PetscPrintf(m_field.get_comm(), "skin_edges_boundary %u\n",
444  skin_edges_boundary.size());
445  }
446  //take all edges on skin faces (i.e. skin surface)
447  Range skin_faces_edges; //edges from skin faces of 3d ents
448  rval = moab.get_adjacencies(skin_faces, 1, false, skin_faces_edges,
449  moab::Interface::UNION);
451  if(mesh_bit_level.any()) {
452  skin_faces_edges = intersect(skin_faces_edges,mesh_level_edges);
453  }
454  if(verb>=VERY_VERBOSE) {
455  PetscPrintf(m_field.get_comm(),"skin_faces_edges %u\n",skin_faces_edges.size());
456  }
457  //note: that skin faces edges do not contain internal boundary
458  //note: that prisms are not included in ents3d, so if ents3d have border with other inteface is like external boundary
459  //skin edges boundary are internal edge <- skin_faces_edges contains edges which are on the body boundary <- that is the trick
460  skin_edges_boundary = subtract(skin_edges_boundary,skin_faces_edges); // from skin edges subtract edges from skin faces of 3d ents (only internal edges)
461 
462  if(verb>=VERY_VERBOSE) {
463  EntityHandle out_meshset;
464  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
465  rval = moab.add_entities(out_meshset,triangles); CHKERRQ_MOAB(rval);
466  rval = moab.write_file("triangles.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
467  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
468  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
469  rval = moab.add_entities(out_meshset,ents3d); CHKERRQ_MOAB(rval);
470  rval = moab.write_file("ents3d.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
471  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
472  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
473  rval = moab.add_entities(out_meshset,skin_edges_boundary); CHKERRQ_MOAB(rval);
474  rval = moab.write_file("skin_edges_boundary.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
475  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
476  }
477  if(verb>=VERY_VERBOSE) {
478  PetscPrintf(
479  m_field.get_comm(),"subtract skin_edges_boundary %u\n",skin_edges_boundary.size()
480  );
481  }
482 
483  //Get nodes on boundary edge
484  Range skin_nodes_boundary;
485  rval = moab.get_connectivity(skin_edges_boundary, skin_nodes_boundary, true);
487  //Remove node which are boundary with other existing interface
488  Range prisms_nodes;
489  rval =
490  moab.get_connectivity(ents3d.subset_by_type(MBPRISM), prisms_nodes, true);
492  skin_nodes_boundary = subtract(skin_nodes_boundary,prisms_nodes);
493  if(verb>=VERY_VERBOSE) {
494  PetscPrintf(m_field.get_comm(), "subtract skin_nodes_boundary %u\n",
495  skin_nodes_boundary.size());
496  }
497  //use nodes on body boundary and interface (without internal boundary) to find adjacent tets
498  Range nodes_without_front = subtract(nodes,skin_nodes_boundary); // nodes_without_front adjacent to all split face edges except those on internal edge
499  if(verb>=VERY_VERBOSE) {
500  PetscPrintf(m_field.get_comm(),
501  "adj. node if ents3d but not on the internal edge %u\n",
502  nodes_without_front.size());
503  }
504 
505  Range skin_nodes_boundary_tris;
506  rval = moab.get_adjacencies(
507  skin_nodes_boundary,2,false,skin_nodes_boundary_tris,moab::Interface::UNION
508  ); CHKERRG(ierr);
509  Range skin_nodes_boundary_tris_nodes;
510  rval = moab.get_connectivity(skin_nodes_boundary_tris,
511  skin_nodes_boundary_tris_nodes, true);
513  skin_nodes_boundary_tris_nodes =
514  subtract(skin_nodes_boundary_tris_nodes, skin_nodes_boundary);
515  Range skin_nodes_boundary_tris_nodes_tris;
516  rval = moab.get_adjacencies(skin_nodes_boundary_tris_nodes, 2, false,
517  skin_nodes_boundary_tris_nodes_tris,
518  moab::Interface::UNION);
519  CHKERRG(ierr);
520  // Triangle which has tree nodes on front boundary
521  skin_nodes_boundary_tris = intersect(
522  triangles,
523  subtract(skin_nodes_boundary_tris,skin_nodes_boundary_tris_nodes_tris)
524  );
525  faces_with_three_nodes_on_front.swap(skin_nodes_boundary_tris);
526 
528 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:497
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:483
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:526
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:490
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  CHKERR m_field.getInterface(meshsets_manager_ptr);
49  CubitMeshSet_multiIndex::index<
50  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator miit =
51  meshsets_manager_ptr->getMeshsetsMultindex()
52  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
53  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
54  if (miit !=
55  meshsets_manager_ptr->getMeshsetsMultindex()
56  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
57  .end()) {
58  CHKERR getSides(miit->meshset, mesh_bit_level, recursive, verb);
59  } else {
60  SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_FOUND,"msId is not there");
61  }
63 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
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 ...
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403

◆ 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 65 of file PrismInterface.cpp.

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

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

535  {
536 
537  Interface &m_field = cOre;
538  MeshsetsManager *meshsets_manager_ptr;
540  CHKERR m_field.getInterface(meshsets_manager_ptr); CHKERRG(ierr);
541  CubitMeshSet_multiIndex::index<
542  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator miit =
543  meshsets_manager_ptr->getMeshsetsMultindex()
544  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
545  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
546  if (miit !=
547  meshsets_manager_ptr->getMeshsetsMultindex()
548  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
549  .end()) {
550  CHKERR splitSides(meshset, bit, miit->meshset, add_interface_entities,
551  recursive, verb);
552  } else {
553  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"msId is not there");
554  }
556 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
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:526
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403

◆ 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 558 of file PrismInterface.cpp.

562  {
563 
565  CHKERR splitSides(meshset, bit, BitRefLevel(), BitRefLevel(), sideset,
566  add_interface_entities, recursive, verb);
568 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
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 ...
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
#define CHKERR
Inline error check.
Definition: definitions.h:578
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:403

◆ 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 569 of file PrismInterface.cpp.

573  {
574 
575  Interface &m_field = cOre;
576  moab::Interface &moab = m_field.get_moab();
577  const RefEntity_multiIndex *refined_ents_ptr;
579 
580  CHKERR m_field.get_ref_ents(&refined_ents_ptr);
581 
582  std::vector<EntityHandle> children;
583  //get children meshsets
584  CHKERR moab.get_child_meshsets(sideset,children);
585  if(children.size()!=3) {
586  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
587  "should be 3 child meshsets, each of them contains tets on two "
588  "sides of interface");
589  }
590  //3d ents on "father" side
591  Range side_ents3d;
592  CHKERR moab.get_entities_by_handle(children[0], side_ents3d, false);
593  //3d ents on "mather" side
594  Range other_ents3d;
595  CHKERR moab.get_entities_by_handle(children[1], other_ents3d, false);
596  //faces of interface
597  Range triangles;
598  CHKERR moab.get_entities_by_type(sideset, MBTRI, triangles, recursive);
599  Range side_ents3d_tris;
600  CHKERR moab.get_adjacencies(side_ents3d, 2, true, side_ents3d_tris,
601  moab::Interface::UNION);
602  triangles = intersect(triangles,side_ents3d_tris);
603  //nodes on interface but not on crack front (those should not be splitted)
604  Range nodes;
605  CHKERR moab.get_entities_by_type(children[2], MBVERTEX, nodes, false);
606  Range meshset_3d_ents,meshset_2d_ents;
607  CHKERR moab.get_entities_by_dimension(meshset, 3, meshset_3d_ents, true);
608  Range meshset_tets = meshset_3d_ents.subset_by_type(MBTET);
609  CHKERR moab.get_adjacencies(meshset_tets, 2, false, meshset_2d_ents,
610  moab::Interface::UNION);
611  side_ents3d = intersect(meshset_3d_ents,side_ents3d);
612  other_ents3d = intersect(meshset_3d_ents,other_ents3d);
613  triangles = intersect(meshset_2d_ents,triangles);
614  if (verb >= VERY_VERBOSE) {
615  PetscPrintf(m_field.get_comm(), "split sides triangles %u\n",
616  triangles.size());
617  PetscPrintf(m_field.get_comm(), "split sides side_ents3d %u\n",
618  side_ents3d.size());
619  PetscPrintf(m_field.get_comm(), "split sides nodes %u\n", nodes.size());
620  }
621 
622  typedef std::map<EntityHandle, /*node on "mather" side*/
623  EntityHandle /*node on "father" side*/
624  >
625  MapNodes;
626  MapNodes map_nodes;
627 
628  // Map nodes on sides, set parent node and set bit ref level
629  {
630  struct CreateSideNodes {
631  MoFEM::Core &cOre;
632  MoFEM::Interface &m_field;
633  std::vector<EntityHandle> splitNodes;
634  std::vector<double> splitCoords[3];
635 
637 
638  CreateSideNodes(MoFEM::Core &core, int split_size = 0)
639  : cOre(core), m_field(core) {
640  splitNodes.reserve(split_size);
641  for (int dd = 0; dd != 3; ++dd) {
642  splitCoords[dd].reserve(split_size);
643  }
644  }
645  MoFEMErrorCode operator()(const double coords[], const EntityHandle n) {
647  splitNodes.push_back(n);
648  for (int dd = 0; dd != 3; ++dd) {
649  splitCoords[dd].push_back(coords[dd]);
650  }
652  }
653 
654  MoFEMErrorCode operator()(const BitRefLevel &bit, MapNodes &map_nodes) {
656  int num_nodes = splitNodes.size();
657  vector<double *> arrays_coord;
658  EntityHandle startv;
659  {
660  ReadUtilIface *iface;
661  CHKERR m_field.get_moab().query_interface(iface);
662  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays_coord);
663  }
664  Range verts(startv, startv + num_nodes - 1);
665  for (int dd = 0; dd != 3; ++dd) {
666  std::copy(splitCoords[dd].begin(), splitCoords[dd].end(),
667  arrays_coord[dd]);
668  }
669  for (int nn = 0; nn != num_nodes; ++nn) {
670  map_nodes[splitNodes[nn]] = verts[nn];
671  }
672  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
673  verts, &*splitNodes.begin());
674  CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
675  verts, bit, QUIET);
676 
678  }
679  };
680  CreateSideNodes create_side_nodes(cOre, nodes.size());
681 
683  struct CreateParentEntView {
684  typedef RefEntity_multiIndex::index<
685  Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
686  MoFEMErrorCode operator()(const BitRefLevel &bit, const BitRefLevel &mask,
687  const RefEntity_multiIndex *refined_ents_ptr,
689  &ref_parent_ents_view) const {
691  const RefEntsByComposite &ref_ents =
692  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
693  RefEntsByComposite::iterator miit;
694  RefEntsByComposite::iterator hi_miit;
695  // view by parent type (VERTEX)
696  miit = ref_ents.lower_bound(boost::make_tuple(MBVERTEX, MBVERTEX));
697  hi_miit = ref_ents.upper_bound(boost::make_tuple(MBVERTEX, MBVERTEX));
698  for (; miit != hi_miit; miit++) {
699  if (((*miit)->getBitRefLevel() & mask) == (*miit)->getBitRefLevel()) {
700  if (((*miit)->getBitRefLevel() & bit).any()) {
701  std::pair<RefEntity_multiIndex_view_by_hashed_parent_entity::iterator,
702  bool>
703  p_ref_ent_view;
704  p_ref_ent_view = ref_parent_ents_view.insert(*miit);
705  if (!p_ref_ent_view.second) {
706  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
707  "non unique insertion");
708  }
709  }
710  }
711  }
713  }
714  };
715  CHKERR CreateParentEntView()(inhered_from_bit_level,
716  inhered_from_bit_level_mask, refined_ents_ptr,
717  ref_parent_ents_view);
718 
719  // add new nodes on interface and create map
720  Range add_bit_nodes;
721  for (Range::const_iterator nit = nodes.begin(); nit != nodes.end(); ++nit) {
722 
723  // find ref enet by parent node
724  RefEntity_multiIndex::iterator miit_ref_ent =
725  refined_ents_ptr->find(*nit);
726  if (miit_ref_ent == refined_ents_ptr->end()) {
727  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
728  "can not find node in MoFEM database");
729  }
730  EntityHandle child_entity = 0;
731  RefEntity_multiIndex_view_by_hashed_parent_entity::iterator child_it =
732  ref_parent_ents_view.find(*nit);
733  if (child_it != ref_parent_ents_view.end()) {
734  child_entity = (*child_it)->getRefEnt();
735  }
736 
737  if (child_entity == 0) {
738  double coords[3];
739  CHKERR moab.get_coords(&*nit, 1, coords);
740  CHKERR create_side_nodes(coords, *nit);
741  } else {
742  map_nodes[*nit] = child_entity;
743  add_bit_nodes.insert(child_entity);
744  }
745  }
746  add_bit_nodes.merge(nodes);
747  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit_nodes,
748  bit);
749  CHKERR create_side_nodes(bit, map_nodes);
750  }
751 
752  // crete meshset for new mesh bit level
753  EntityHandle meshset_for_bit_level;
754  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
755  meshset_for_bit_level);
756  // subtract those elements which will be refined, i.e. disconnected form other
757  // side elements, and connected to new prisms, if they area created
758  meshset_3d_ents = subtract(meshset_3d_ents,side_ents3d);
759  CHKERR moab.add_entities(meshset_for_bit_level, meshset_3d_ents);
760 
761  // create new 3d ents on "father" side
762  Range new_3d_ents;
763  for (Range::iterator eit3d = side_ents3d.begin(); eit3d != side_ents3d.end();
764  eit3d++) {
765  RefEntity_multiIndex::iterator miit_ref_ent =
766  refined_ents_ptr->find(*eit3d);
767  if (miit_ref_ent == refined_ents_ptr->end()) {
768  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
769  "tet not in database");
770  }
771  int num_nodes;
772  const EntityHandle* conn;
773  CHKERR moab.get_connectivity(*eit3d, conn, num_nodes, true);
774  EntityHandle new_conn[num_nodes];
775  int nb_new_conn = 0;
776  for (int ii = 0; ii < num_nodes; ii++) {
777  std::map<EntityHandle, EntityHandle>::iterator mit =
778  map_nodes.find(conn[ii]);
779  if (mit != map_nodes.end()) {
780  new_conn[ii] = mit->second;
781  nb_new_conn++;
782  if (verb >= VERY_NOISY) {
783  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
784  new_conn[ii]);
785  }
786  } else {
787  new_conn[ii] = conn[ii];
788  }
789  }
790  if (nb_new_conn == 0) {
791  // Add this tet to bit ref level
792  CHKERR moab.add_entities(meshset_for_bit_level, &*eit3d, 1);
793  continue;
794  }
795 
796  // here is created new or prism is on interface
797  EntityHandle existing_ent = 0;
798  /* check if tet element with new connectivity is in database*/
799  RefEntity_multiIndex::index<Ent_Ent_mi_tag>::type::iterator child_iit,
800  hi_child_iit;
801  child_iit =
802  refined_ents_ptr->get<Ent_Ent_mi_tag>().lower_bound(*eit3d);
803  hi_child_iit =
804  refined_ents_ptr->get<Ent_Ent_mi_tag>().upper_bound(*eit3d);
805 
806  // Check if child entity has the same connectivity
807  for(;child_iit!=hi_child_iit;child_iit++) {
808  const EntityHandle* conn_ref_tet;
809  CHKERR moab.get_connectivity(child_iit->get()->getRefEnt(), conn_ref_tet,
810  num_nodes, true);
811  int nn = 0;
812  for(;nn<num_nodes;nn++) {
813  if(conn_ref_tet[nn]!=new_conn[nn]) {
814  break;
815  }
816  }
817  if(nn == num_nodes) {
818  if(existing_ent != 0) {
819  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
820  "database inconsistency");
821  }
822  existing_ent = child_iit->get()->getRefEnt();
823  }
824  }
825 
826  switch (moab.type_from_handle(*eit3d)) {
827  case MBTET: {
828  RefEntity_multiIndex::iterator child_it;
829  EntityHandle tet;
830  if(existing_ent == 0) {
831  Range new_conn_tet;
832  CHKERR moab.get_adjacencies(new_conn, 4, 3, false, new_conn_tet);
833  if(new_conn_tet.empty()) {
834  CHKERR moab.create_element(MBTET, new_conn, 4, tet);
835  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &tet, 1,
836  &*eit3d);
837  } else {
838  // RefEntity_multiIndex::index<Ent_mi_tag>::type::iterator rit;
839  // rit = refined_ents_ptr->get<Ent_mi_tag>().find(*eit3d);
840  // if(rit==refined_ents_ptr->get<Ent_mi_tag>().end()) {
841  // SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
842  // "can't find this in database");
843  // }
844  // FIXME: That takes firs element form the list. Should throw error
845  // if is more than one or handle it properly.
846  RefEntity_multiIndex::index<Ent_mi_tag>::type::iterator new_rit =
847  refined_ents_ptr->get<Ent_mi_tag>().find(*new_conn_tet.begin());
848  if(new_rit==refined_ents_ptr->get<Ent_mi_tag>().end()) {
849  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
850  "can't find this in database");
851  }
852  tet = *new_conn_tet.begin();
853  }
854  } else {
855  tet = existing_ent;
856  }
857  CHKERR moab.add_entities(meshset_for_bit_level, &tet, 1);
858  new_3d_ents.insert(tet);
859  } break;
860  case MBPRISM: {
861  EntityHandle prism;
862  if(existing_ent == 0) {
863  Range new_conn_prism;
864  CHKERR moab.get_adjacencies(new_conn, 6, 3, false, new_conn_prism);
865  if(new_conn_prism.empty()) {
866  CHKERR moab.create_element(MBPRISM, new_conn, 6, prism);
867  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &prism, 1,
868  &*eit3d);
869  } else {
870  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
871  "It is prism with such connectivity, that case has to be "
872  "handled but this is not implemented");
873  }
874  } else {
875  prism = existing_ent;
876  }
877  CHKERR moab.add_entities(meshset_for_bit_level, &prism, 1);
878  new_3d_ents.insert(prism);
879  } break;
880  default:
881  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
882  "not implemented element");
883  }
884  }
885 
886  Range new_ents;
887  //create new entities by adjacencies form new tets
888  CHKERR moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET), 2, true,
889  new_ents, moab::Interface::UNION);
890  CHKERR moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET), 1, true,
891  new_ents, moab::Interface::UNION);
892 
893  //Tags for setting side
894  Tag th_interface_side;
895  const int def_side[] = {0};
896  CHKERR moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
897  th_interface_side, MB_TAG_CREAT | MB_TAG_SPARSE,
898  def_side);
899 
900  struct SetParent {
901  map<EntityHandle, EntityHandle> parentsToChange;
902  MoFEMErrorCode operator()(const EntityHandle ent, const EntityHandle parent,
903  const RefEntity_multiIndex *ref_ents_ptr,
904  MoFEM::Core &cOre) {
905  MoFEM::Interface &m_field = cOre;
907  RefEntity_multiIndex::iterator it = ref_ents_ptr->find(ent);
908  if (it != ref_ents_ptr->end()) {
909  if (it->get()->getParentEnt() != parent && ent != parent) {
910  parentsToChange[ent] = parent;
911  }
912  } else {
913  if (ent != parent) {
914  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
915  &ent, 1, &parent);
916  }
917  }
919  }
920  MoFEMErrorCode operator()(const RefEntity_multiIndex *ref_ents_ptr) {
922  for (map<EntityHandle, EntityHandle>::iterator mit =
923  parentsToChange.begin();
924  mit != parentsToChange.end(); ++mit) {
925  RefEntity_multiIndex::iterator it = ref_ents_ptr->find(mit->first);
926  bool success = const_cast<RefEntity_multiIndex *>(ref_ents_ptr)
927  ->modify(it, RefEntity_change_parent(mit->second));
928  if (!success) {
929  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
930  "impossible to set parent");
931  }
932  }
934  }
935  };
936  SetParent set_parent;
937 
938  // add new edges and triangles to mofem database
939  Range ents;
940  CHKERR moab.get_adjacencies(triangles, 1, false, ents,
941  moab::Interface::UNION);
942  ents.insert(triangles.begin(),triangles.end());
943  for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
944  int num_nodes;
945  const EntityHandle *conn;
946  CHKERR moab.get_connectivity(*eit, conn, num_nodes, true);
947  EntityHandle new_conn[num_nodes];
948  int nb_new_conn = 0;
949  for (int ii = 0; ii != num_nodes; ++ii) {
950  std::map<EntityHandle, EntityHandle>::iterator mit =
951  map_nodes.find(conn[ii]);
952  if (mit != map_nodes.end()) {
953  new_conn[ii] = mit->second;
954  nb_new_conn++;
955  if (verb >= VERY_NOISY) {
956  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
957  new_conn[ii]);
958  }
959  } else {
960  new_conn[ii] = conn[ii];
961  }
962  }
963  if (nb_new_conn == 0)
964  continue;
965  RefEntity_multiIndex::iterator miit_ref_ent = refined_ents_ptr->find(*eit);
966  if (miit_ref_ent == refined_ents_ptr->end()) {
967  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
968  "this entity (edge or tri) should be already in database");
969  }
970  Range new_ent; // contains all entities (edges or triangles) added to mofem
971  // database
972  switch (moab.type_from_handle(*eit)) {
973  case MBTRI: {
974  // get entity based on its connectivity
975  CHKERR moab.get_adjacencies(new_conn, 3, 2, false, new_ent);
976  if (new_ent.size() != 1)
977  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
978  "this tri should be in moab database");
979  int new_side = 1;
980  CHKERR moab.tag_set_data(th_interface_side, &*new_ent.begin(), 1,
981  &new_side);
982  if (verb >= VERY_VERBOSE)
983  PetscPrintf(m_field.get_comm(), "new_ent %u\n", new_ent.size());
984  // add prism element
985  if (add_interface_entities) {
986  if (inhered_from_bit_level.any()) {
987  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
988  "not implemented for inhered_from_bit_level");
989  }
990  // set prism connectivity
991  EntityHandle prism_conn[6] = {conn[0], conn[1], conn[2],
992  new_conn[0], new_conn[1], new_conn[2]};
993  EntityHandle prism;
994  CHKERR moab.create_element(MBPRISM, prism_conn, 6, prism);
995  CHKERR moab.add_entities(meshset_for_bit_level, &prism, 1);
996  }
997  } break;
998  case MBEDGE: {
999  CHKERR moab.get_adjacencies(new_conn, 2, 1, false, new_ent);
1000  if (new_ent.size() != 1) {
1001  SETERRQ2(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1002  "this edge should be in moab database "
1003  "new_ent.size() = %u nb_new_conn = %d",
1004  new_ent.size(), nb_new_conn);
1005  }
1006  } break;
1007  default:
1008  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1009  "Houston we have a problem !!!");
1010  }
1011  if (new_ent.size() != 1) {
1012  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1013  "new_ent.size() = %u, size always should be 1", new_ent.size());
1014  }
1015  CHKERR set_parent(new_ent[0], *eit, refined_ents_ptr, cOre);
1016  }
1017 
1018  // all other entities, some ents like triangles and faces on the side of tets
1019  Range side_adj_faces_and_edges;
1020  CHKERR moab.get_adjacencies(side_ents3d.subset_by_type(MBTET), 1, true,
1021  side_adj_faces_and_edges, moab::Interface::UNION);
1022  CHKERR moab.get_adjacencies(side_ents3d.subset_by_type(MBTET), 2, true,
1023  side_adj_faces_and_edges, moab::Interface::UNION);
1024  for (Range::iterator eit = side_adj_faces_and_edges.begin();
1025  eit != side_adj_faces_and_edges.end(); ++eit) {
1026  int num_nodes;
1027  const EntityHandle* conn;
1028  CHKERR moab.get_connectivity(*eit, conn, num_nodes, true);
1029  EntityHandle new_conn[num_nodes];
1030  int nb_new_conn = 0;
1031  for (int ii = 0; ii < num_nodes; ii++) {
1032  std::map<EntityHandle, EntityHandle>::iterator mit =
1033  map_nodes.find(conn[ii]);
1034  if (mit != map_nodes.end()) {
1035  new_conn[ii] = mit->second;
1036  nb_new_conn++;
1037  if (verb >= VERY_NOISY) {
1038  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
1039  new_conn[ii]);
1040  }
1041  } else {
1042  new_conn[ii] = conn[ii];
1043  }
1044  }
1045  if (nb_new_conn == 0)
1046  continue;
1047  RefEntity_multiIndex::iterator miit_ref_ent = refined_ents_ptr->find(*eit);
1048  if (miit_ref_ent == refined_ents_ptr->end()) {
1049  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1050  "entity should be in MoFem database, num_nodes = %d", num_nodes);
1051  }
1052  Range new_ent;
1053  switch (moab.type_from_handle(*eit)) {
1054  case MBEDGE:
1055  CHKERR moab.get_adjacencies(new_conn, 2, 1, false, new_ent);
1056  break;
1057  case MBTRI:
1058  CHKERR moab.get_adjacencies(new_conn, 3, 2, false, new_ent);
1059  break;
1060  default:
1061  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1062  "Houston we have a problem");
1063  }
1064  if (new_ent.size() != 1) {
1065  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1066  "database inconsistency, new_ent.size() = %u", new_ent.size());
1067  }
1068  CHKERR set_parent(new_ent[0], *eit, refined_ents_ptr, cOre);
1069  }
1070 
1071  // add new prisms which parents are part of other intefaces
1072  Range new_3d_prims = new_3d_ents.subset_by_type(MBPRISM);
1073  for (Range::iterator pit = new_3d_prims.begin(); pit != new_3d_prims.end();
1074  ++pit) {
1075  // get parent entity
1076  EntityHandle parent_prism;
1077  CHKERR moab.tag_get_data(cOre.get_th_RefParentHandle(), &*pit, 1,
1078  &parent_prism);
1079  if (moab.type_from_handle(parent_prism) != MBPRISM) {
1080  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1081  "this prism should have parent which is prism as well");
1082  }
1083  int num_nodes;
1084  // parent prism
1085  const EntityHandle *conn_parent;
1086  CHKERR moab.get_connectivity(parent_prism, conn_parent, num_nodes, true);
1087  Range face_side3_parent, face_side4_parent;
1088  CHKERR moab.get_adjacencies(conn_parent, 3, 2, false, face_side3_parent);
1089  CHKERR moab.get_adjacencies(&conn_parent[3], 3, 2, false, face_side4_parent);
1090  if (face_side3_parent.size() != 1) {
1091  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1092  "parent face3.size() = %u", face_side3_parent.size());
1093  }
1094  if (face_side4_parent.size() != 1) {
1095  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1096  "parent face4.size() = %u", face_side4_parent.size());
1097  }
1098  // new prism
1099  const EntityHandle *conn;
1100  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
1101  Range face_side3, face_side4;
1102  CHKERR moab.get_adjacencies(conn, 3, 2, false, face_side3);
1103  CHKERR moab.get_adjacencies(&conn[3], 3, 2, false, face_side4);
1104  if (face_side3.size() != 1) {
1105  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "face3 is missing");
1106  }
1107  if (face_side4.size() != 1) {
1108  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "face4 is missing");
1109  }
1110  std::vector<EntityHandle> face(2), parent_face(2);
1111  face[0] = *face_side3.begin();
1112  face[1] = *face_side4.begin();
1113  parent_face[0] = *face_side3_parent.begin();
1114  parent_face[1] = *face_side4_parent.begin();
1115  for (int ff = 0; ff != 2; ++ff) {
1116  if (parent_face[ff] == face[ff])
1117  continue;
1118  int interface_side;
1119  CHKERR moab.tag_get_data(th_interface_side, &parent_face[ff], 1,
1120  &interface_side);
1121  CHKERR moab.tag_set_data(th_interface_side, &face[ff], 1,
1122  &interface_side);
1123  EntityHandle parent_tri;
1124  CHKERR moab.tag_get_data(cOre.get_th_RefParentHandle(), &face[ff], 1,
1125  &parent_tri);
1126  if (parent_tri != parent_face[ff]) {
1127  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1128  "wrong parent %lu", parent_tri);
1129  }
1130  }
1131  }
1132  // finalise by adding new tets and prism ti bit level
1133 
1134  CHKERR set_parent(refined_ents_ptr);
1135  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
1136  meshset_for_bit_level, 3, bit);
1137  CHKERR moab.delete_entities(&meshset_for_bit_level, 1);
1138  CHKERR moab.clear_meshset(&children[0], 3);
1140 }
Tag get_th_RefParentHandle() const
Definition: Core.hpp:150
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Common.hpp:60
virtual moab::Interface & get_moab()=0
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
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:459
Core (interface) class.
Definition: Core.hpp:50
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
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< 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:578
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:403
virtual MPI_Comm & get_comm() const =0

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: