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

Create interface from given surface and insert flat prisms in-between. 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)
 Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshset. More...
 
MoFEMErrorCode getSides (const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
 Store tetrahedra from each side of the interface in two child meshsets of the parent meshset. More...
 
MoFEMErrorCode findFacesWithThreeNodesOnInternalSurfaceSkin (const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, Range &faces_with_three_nodes_on_front, int verb=QUIET)
 Find triangles which have 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 tetrahedra on both sides of the interface and insert flat prisms in-between. More...
 
MoFEMErrorCode 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 tetrahedra on both sides of the interface and insert flat prisms in-between. 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 tetrahedra on both sides of the interface and insert flat prisms in-between. More...
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (const MOFEMuuid &uuid, bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE , bool VERIFY = false>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, IFACE *&iface) const
 Get interface by uuid and return reference to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 
virtual MoFEMErrorCode getLibVersion (Version &version) const
 Get library version. More...
 
virtual const MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version) const
 Get database major version. More...
 
virtual MoFEMErrorCode getInterfaceVersion (Version &version) const
 Get database major version. More...
 
template<>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, UnknownInterface *&iface) const
 

Public Attributes

MoFEM::CorecOre
 

Additional Inherited Members

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

Detailed Description

Create interface from given surface and insert flat prisms in-between.

Examples
forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, and split_sideset.cpp.

Definition at line 36 of file PrismInterface.hpp.

Constructor & Destructor Documentation

◆ PrismInterface()

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

Definition at line 36 of file PrismInterface.cpp.

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

◆ ~PrismInterface()

MoFEM::PrismInterface::~PrismInterface ( )

destructor

Definition at line 45 of file PrismInterface.hpp.

45 {}

Member Function Documentation

◆ findFacesWithThreeNodesOnInternalSurfaceSkin()

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

Find triangles which have three nodes on internal surface skin.

Internal surface skin is a set of all edges on the boundary of a given surface inside the body. This set of edges is also called the surface front. If a triangle has three nodes on the surface front, none of these nodes can be split. Therefore, such a triangle cannot be split and should be removed from the surface.

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

Definition at line 397 of file PrismInterface.cpp.

399  {
400  Interface &m_field = cOre;
401  moab::Interface &moab = m_field.get_moab();
403 
404  Range mesh_level_ents3d;
405  Range mesh_level_edges, mesh_level_tris;
406  if (mesh_bit_level.any()) {
407  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
408  mesh_bit_level, BitRefLevel().set(), MBTET, mesh_level_ents3d);
409 
410  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
411  mesh_bit_level, BitRefLevel().set(), MBTRI, mesh_level_tris);
412 
413  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
414  mesh_bit_level, BitRefLevel().set(), MBEDGE, mesh_level_edges);
415  }
416 
417  Skinner skin(&moab);
418  // get interface triangles from side set
419  Range triangles;
420  CHKERR moab.get_entities_by_type(sideset, MBTRI, triangles, recursive);
421 
422  if (mesh_bit_level.any()) {
423  triangles = intersect(triangles, mesh_level_tris);
424  }
425  if (verb >= VERBOSE) {
426  PetscPrintf(m_field.get_comm(), "Nb. of triangles in set %u\n",
427  triangles.size());
428  }
429  // get nodes, edges and 3d ents (i.e. tets and prisms)
430  Range nodes; // nodes from triangles
431  CHKERR moab.get_connectivity(triangles, nodes, true);
432 
433  Range ents3d; // 3d ents from nodes
434  CHKERR moab.get_adjacencies(nodes, 3, false, ents3d, moab::Interface::UNION);
435 
436  if (mesh_bit_level.any()) {
437  ents3d = intersect(ents3d, mesh_level_ents3d);
438  }
439  // take skin faces
440  Range skin_faces; // skin faces from 3d ents
441  CHKERR skin.find_skin(0, ents3d.subset_by_type(MBTET), false, skin_faces);
442 
443  // take skin edges (boundary of surface if there is any)
444  Range skin_edges_boundary; // skin edges from triangles
445  CHKERR skin.find_skin(0, triangles, false, skin_edges_boundary);
446 
447  if (verb >= VERY_VERBOSE) {
448  PetscPrintf(m_field.get_comm(), "skin_edges_boundary %u\n",
449  skin_edges_boundary.size());
450  }
451  // take all edges on skin faces (i.e. skin surface)
452  Range skin_faces_edges; // edges from skin faces of 3d ents
453  CHKERR moab.get_adjacencies(skin_faces, 1, false, skin_faces_edges,
454  moab::Interface::UNION);
455 
456  if (mesh_bit_level.any()) {
457  skin_faces_edges = intersect(skin_faces_edges, mesh_level_edges);
458  }
459  if (verb >= VERY_VERBOSE) {
460  PetscPrintf(m_field.get_comm(), "skin_faces_edges %u\n",
461  skin_faces_edges.size());
462  }
463  // note: that skin faces edges do not contain internal boundary
464  // note: that prisms are not included in ents3d, so if ents3d have border with
465  // other inteface is like external boundary skin edges boundary are internal
466  // edge <- skin_faces_edges contains edges which are on the body boundary <-
467  // that is the trick
468  skin_edges_boundary =
469  subtract(skin_edges_boundary,
470  skin_faces_edges); // from skin edges subtract edges from skin
471  // faces of 3d ents (only internal edges)
472 
473  if (verb >= VERY_VERBOSE) {
474  EntityHandle out_meshset;
475  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, out_meshset);
476  CHKERR moab.add_entities(out_meshset, triangles);
477  CHKERR moab.write_file("triangles.vtk", "VTK", "", &out_meshset, 1);
478  CHKERR moab.delete_entities(&out_meshset, 1);
479  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, out_meshset);
480  CHKERR moab.add_entities(out_meshset, ents3d);
481  CHKERR moab.write_file("ents3d.vtk", "VTK", "", &out_meshset, 1);
482  CHKERR moab.delete_entities(&out_meshset, 1);
483  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, out_meshset);
484  CHKERR moab.add_entities(out_meshset, skin_edges_boundary);
485  CHKERR moab.write_file("skin_edges_boundary.vtk", "VTK", "", &out_meshset,
486  1);
487  CHKERR moab.delete_entities(&out_meshset, 1);
488  }
489  if (verb >= VERY_VERBOSE) {
490  PetscPrintf(m_field.get_comm(), "subtract skin_edges_boundary %u\n",
491  skin_edges_boundary.size());
492  }
493 
494  // Get nodes on boundary edge
495  Range skin_nodes_boundary;
496  CHKERR moab.get_connectivity(skin_edges_boundary, skin_nodes_boundary, true);
497 
498  // Remove node which are boundary with other existing interface
499  Range prisms_nodes;
500  CHKERR
501  moab.get_connectivity(ents3d.subset_by_type(MBPRISM), prisms_nodes, true);
502 
503  skin_nodes_boundary = subtract(skin_nodes_boundary, prisms_nodes);
504  if (verb >= VERY_VERBOSE) {
505  PetscPrintf(m_field.get_comm(), "subtract skin_nodes_boundary %u\n",
506  skin_nodes_boundary.size());
507  }
508  // use nodes on body boundary and interface (without internal boundary) to
509  // find adjacent tets
510  Range nodes_without_front = subtract(
511  nodes, skin_nodes_boundary); // nodes_without_front adjacent to all split
512  // face edges except those on internal edge
513  if (verb >= VERY_VERBOSE) {
514  PetscPrintf(m_field.get_comm(),
515  "adj. node if ents3d but not on the internal edge %u\n",
516  nodes_without_front.size());
517  }
518 
519  Range skin_nodes_boundary_tris;
520  CHKERR moab.get_adjacencies(skin_nodes_boundary, 2, false,
521  skin_nodes_boundary_tris, moab::Interface::UNION);
522  Range skin_nodes_boundary_tris_nodes;
523  CHKERR moab.get_connectivity(skin_nodes_boundary_tris,
524  skin_nodes_boundary_tris_nodes, true);
525 
526  skin_nodes_boundary_tris_nodes =
527  subtract(skin_nodes_boundary_tris_nodes, skin_nodes_boundary);
528  Range skin_nodes_boundary_tris_nodes_tris;
529  CHKERR moab.get_adjacencies(skin_nodes_boundary_tris_nodes, 2, false,
530  skin_nodes_boundary_tris_nodes_tris,
531  moab::Interface::UNION);
532 
533  // Triangle which has tree nodes on front boundary
534  skin_nodes_boundary_tris =
535  intersect(triangles, subtract(skin_nodes_boundary_tris,
536  skin_nodes_boundary_tris_nodes_tris));
537  faces_with_three_nodes_on_front.swap(skin_nodes_boundary_tris);
538 
540 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ 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 
)

Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshset.

Additional third child meshset contains nodes which can be split and skin edges

Parameters
msIdId of meshset
cubit_bc_typetype of meshset (NODESET, SIDESET or BLOCKSET and more)
mesh_bit_levelinterface is added on this bit level
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
if bit_level == BitRefLevel.set() then interface will be added on all bit levels
Examples
forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, and split_sideset.cpp.

Definition at line 39 of file PrismInterface.cpp.

42  {
43 
44  Interface &m_field = cOre;
45  MeshsetsManager *meshsets_manager_ptr;
47  CHKERR m_field.getInterface(meshsets_manager_ptr);
48  CubitMeshSet_multiIndex::index<
49  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator miit =
50  meshsets_manager_ptr->getMeshsetsMultindex()
51  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
52  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
53  if (miit != meshsets_manager_ptr->getMeshsetsMultindex()
54  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
55  .end()) {
56  CHKERR getSides(miit->meshset, mesh_bit_level, recursive, verb);
57  } else {
58  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "msId is not there");
59  }
61 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ getSides() [2/2]

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

Store tetrahedra from each side of the interface in two child meshsets of the parent meshset.

Additional third child meshset contains nodes which can be split and skin edges

Parameters
sidesetparent meshset with the surface
mesh_bit_levelinterface is added on this bit level
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
if bit_level == BitRefLevel.set() then interface will be added on all bit levels
  1. Get tets adjacent to nodes of the interface meshset.
  2. Take skin faces from these tets and get edges from that skin.
  3. Take skin from triangles of the interface.
  4. Subtract edges of skin faces from skin of triangles in order to get edges in the volume of the body, and not on the interface boundary.
  5. Iterate between all triangles of the interface and find adjacent tets on each side of the interface

Definition at line 63 of file PrismInterface.cpp.

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

◆ 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:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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 tetrahedra on both sides of the interface and insert flat prisms in-between.

Parameters
meshsetvolume meshset containing 3D entities around the interface
bitbit ref level on which new entities will be stored
msIdmeshset ID of the surface
cubit_bc_typetype of meshset (NODESET, SIDESET or BLOCKSET and more)
add_interface_entitiesif true add prism elements at interface
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
Parent meshset must have three child meshsets: two with tetrahedra from each side of the interface, third containing nodes which can be split and skin edges
Examples
forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, and split_sideset.cpp.

Definition at line 542 of file PrismInterface.cpp.

547  {
548 
549  Interface &m_field = cOre;
550  MeshsetsManager *meshsets_manager_ptr;
552  CHKERR m_field.getInterface(meshsets_manager_ptr);
553  CubitMeshSet_multiIndex::index<
554  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator miit =
555  meshsets_manager_ptr->getMeshsetsMultindex()
556  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
557  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
558  if (miit != meshsets_manager_ptr->getMeshsetsMultindex()
559  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
560  .end()) {
561  CHKERR splitSides(meshset, bit, miit->meshset, add_interface_entities,
562  recursive, verb);
563  } else {
564  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "msId is not there");
565  }
567 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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 tetrahedra on both sides of the interface and insert flat prisms in...
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ 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 tetrahedra on both sides of the interface and insert flat prisms in-between.

Parameters
meshsetvolume meshset containing 3D entities around the interface
bitbit ref level on which new entities will be stored
sidesetmeshset with surface
add_interface_entitiesif true add prism elements at interface
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
Parent meshset must have three child meshsets: two with tetrahedra from each side of the interface, third containing nodes which can be split and skin edges

Definition at line 569 of file PrismInterface.cpp.

573  {
574 
576  CHKERR splitSides(meshset, bit, BitRefLevel(), BitRefLevel(), sideset,
577  add_interface_entities, recursive, verb);
579 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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 tetrahedra on both sides of the interface and insert flat prisms in...
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ 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 tetrahedra on both sides of the interface and insert flat prisms in-between.

Parameters
meshsetvolume meshset containing 3D entities around the interface
bitbit ref level on which new entities will be stored
inhered_from_bit_levelinherit nodes and other entities form this bit level
inhered_from_bit_level_maskcorresponding mask
sidesetmeshset with surface
add_interface_entitiesif true add prism elements at interface
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
Parent meshset must have three child meshsets: two with tetrahedra from each side of the interface, third containing nodes which can be split and skin edges
inhered_from_bit_level needs to be specified for some meshsets with interfaces. Some nodes on some refinement levels are dividing edges but not splitting faces. Inheriting those nodes will not split faces.

Definition at line 581 of file PrismInterface.cpp.

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

Member Data Documentation

◆ cOre

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

Definition at line 41 of file PrismInterface.hpp.


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