v0.9.1
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, simple_contact.cpp, split_sideset.cpp, and test_jacobian_of_simple_contact_element.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 389 of file PrismInterface.cpp.

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

◆ 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, simple_contact.cpp, split_sideset.cpp, and test_jacobian_of_simple_contact_element.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:483
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:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413

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

◆ 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:507
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
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, simple_contact.cpp, split_sideset.cpp, and test_jacobian_of_simple_contact_element.cpp.

Definition at line 534 of file PrismInterface.cpp.

539  {
540 
541  Interface &m_field = cOre;
542  MeshsetsManager *meshsets_manager_ptr;
544  CHKERR m_field.getInterface(meshsets_manager_ptr);
545  CubitMeshSet_multiIndex::index<
546  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator miit =
547  meshsets_manager_ptr->getMeshsetsMultindex()
548  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
549  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
550  if (miit != meshsets_manager_ptr->getMeshsetsMultindex()
551  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
552  .end()) {
553  CHKERR splitSides(meshset, bit, miit->meshset, add_interface_entities,
554  recursive, verb);
555  } else {
556  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "msId is not there");
557  }
559 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413

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

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

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

577  {
578  Interface &m_field = cOre;
579  moab::Interface &moab = m_field.get_moab();
580  auto refined_ents_ptr = m_field.get_ref_ents();
582 
583  std::vector<EntityHandle> children;
584  // get children meshsets
585  CHKERR moab.get_child_meshsets(sideset, children);
586  if (children.size() != 3)
587  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
588  "should be 3 child meshsets, each of them contains tets on two "
589  "sides of interface");
590 
591  // 3d ents on "father" side
592  Range side_ents3d;
593  CHKERR moab.get_entities_by_handle(children[0], side_ents3d, false);
594  // 3d ents on "mother" side
595  Range other_ents3d;
596  CHKERR moab.get_entities_by_handle(children[1], other_ents3d, false);
597  // faces of interface
598  Range triangles;
599  CHKERR moab.get_entities_by_type(sideset, MBTRI, triangles, recursive);
600  Range side_ents3d_tris;
601  CHKERR moab.get_adjacencies(side_ents3d, 2, true, side_ents3d_tris,
602  moab::Interface::UNION);
603  triangles = intersect(triangles, side_ents3d_tris);
604  // nodes on interface but not on crack front (those should not be splitted)
605  Range nodes;
606  CHKERR moab.get_entities_by_type(children[2], MBVERTEX, nodes, false);
607  Range meshset_3d_ents, meshset_2d_ents;
608  CHKERR moab.get_entities_by_dimension(meshset, 3, meshset_3d_ents, true);
609  Range meshset_tets = meshset_3d_ents.subset_by_type(MBTET);
610  CHKERR moab.get_adjacencies(meshset_tets, 2, false, meshset_2d_ents,
611  moab::Interface::UNION);
612  side_ents3d = intersect(meshset_3d_ents, side_ents3d);
613  other_ents3d = intersect(meshset_3d_ents, other_ents3d);
614  triangles = intersect(meshset_2d_ents, triangles);
615  if (verb >= VERY_VERBOSE) {
616  PetscPrintf(m_field.get_comm(), "split sides triangles %u\n",
617  triangles.size());
618  PetscPrintf(m_field.get_comm(), "split sides side_ents3d %u\n",
619  side_ents3d.size());
620  PetscPrintf(m_field.get_comm(), "split sides nodes %u\n", nodes.size());
621  }
622 
623 
624  struct PartentAndChild {
625  EntityHandle parent;
626  EntityHandle child;
627  };
628 
629  typedef multi_index_container<
630  PartentAndChild,
631  indexed_by<
632 
633  hashed_unique<
634  member<PartentAndChild, EntityHandle, &PartentAndChild::parent>>,
635 
636  hashed_unique<
637  member<PartentAndChild, EntityHandle, &PartentAndChild::child>>
638 
639  >>
640  ParentChildMI;
641 
642  ParentChildMI parent_child;
643 
644 
645  typedef std::map<EntityHandle, /*node on "mother" side*/
646  EntityHandle /*node on "father" side*/
647  >
648  MapNodes;
649  MapNodes map_nodes, reverse_map_nodes;
650 
651  // Map nodes on sides, set parent node and set bit ref level
652  {
653  struct CreateSideNodes {
654  MoFEM::Core &cOre;
655  MoFEM::Interface &m_field;
656  std::vector<EntityHandle> splitNodes;
657  std::vector<double> splitCoords[3];
658 
660 
661  CreateSideNodes(MoFEM::Core &core, int split_size = 0)
662  : cOre(core), m_field(core) {
663  splitNodes.reserve(split_size);
664  for (auto dd : {0, 1, 2})
665  splitCoords[dd].reserve(split_size);
666  }
667  MoFEMErrorCode operator()(const double coords[], const EntityHandle n) {
669  splitNodes.emplace_back(n);
670  for (auto dd : {0, 1, 2})
671  splitCoords[dd].emplace_back(coords[dd]);
673  }
674 
675  MoFEMErrorCode operator()(const BitRefLevel &bit, MapNodes &map_nodes,
676  MapNodes &reverse_map_nodes) {
677  ReadUtilIface *iface;
679  int num_nodes = splitNodes.size();
680  std::vector<double *> arrays_coord;
681  EntityHandle startv;
682  CHKERR m_field.get_moab().query_interface(iface);
683  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays_coord);
684  Range verts(startv, startv + num_nodes - 1);
685  for (int dd = 0; dd != 3; ++dd)
686  std::copy(splitCoords[dd].begin(), splitCoords[dd].end(),
687  arrays_coord[dd]);
688  for (int nn = 0; nn != num_nodes; ++nn) {
689  map_nodes[splitNodes[nn]] = verts[nn];
690  reverse_map_nodes[verts[nn]] = splitNodes[nn];
691  }
692  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
693  verts, &*splitNodes.begin());
694  CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
695  verts, bit, QUIET);
697  }
698  };
699  CreateSideNodes create_side_nodes(cOre, nodes.size());
700 
702  struct CreateParentEntView {
704  operator()(const BitRefLevel &bit, const BitRefLevel &mask,
705  const RefEntity_multiIndex *refined_ents_ptr,
707  &ref_parent_ents_view) const {
709  auto &ref_ents =
710  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
711  // view by parent type (VERTEX)
712  auto miit = ref_ents.lower_bound(boost::make_tuple(MBVERTEX, MBVERTEX));
713  auto hi_miit =
714  ref_ents.upper_bound(boost::make_tuple(MBVERTEX, MBVERTEX));
715  for (; miit != hi_miit; miit++) {
716  const auto &ent_bit = (*miit)->getBitRefLevel();
717  if ((ent_bit & bit).any() && (ent_bit & mask) == ent_bit) {
718  auto p_ref_ent_view = ref_parent_ents_view.insert(*miit);
719  if (!p_ref_ent_view.second)
720  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
721  "non unique insertion");
722  }
723  }
725  }
726  };
727  if (inhered_from_bit_level.any() && inhered_from_bit_level_mask.any())
728  CHKERR CreateParentEntView()(inhered_from_bit_level,
729  inhered_from_bit_level_mask,
730  refined_ents_ptr, ref_parent_ents_view);
731 
732  // add new nodes on interface and create map
733  Range add_bit_nodes;
734  for (auto pnit = nodes.pair_begin(); pnit != nodes.pair_end(); ++pnit) {
735  auto lo = refined_ents_ptr->lower_bound(pnit->first);
736  auto hi = refined_ents_ptr->upper_bound(pnit->second);
737  if (std::distance(lo, hi) != (pnit->second - pnit->first + 1))
738  SETERRQ(
739  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
740  "Can not find some nodes in database that are split on interface");
741  Range nodes_in_range;
742  insertOrdered(nodes_in_range, RefEntExtractor(), lo, hi);
743  std::vector<double> coords_range(nodes_in_range.size() * 3);
744  CHKERR moab.get_coords(nodes_in_range, &*coords_range.begin());
745  int pos = 0;
746  for (; lo != hi; ++lo, pos += 3) {
747  const EntityHandle node = (*lo)->getRefEnt();
748  EntityHandle child_entity = 0;
749  auto child_it = ref_parent_ents_view.find(node);
750  if (child_it != ref_parent_ents_view.end())
751  child_entity = (*child_it)->getRefEnt();
752  if (child_entity == 0) {
753  CHKERR create_side_nodes(&coords_range[pos], node);
754  } else {
755  map_nodes[node] = child_entity;
756  add_bit_nodes.insert(child_entity);
757  }
758  }
759  }
760  add_bit_nodes.merge(nodes);
761  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit_nodes,
762  bit);
763  CHKERR create_side_nodes(bit, map_nodes, reverse_map_nodes);
764  }
765 
766  // crete meshset for new mesh bit level
767  EntityHandle meshset_for_bit_level;
768  CHKERR moab.create_meshset(MESHSET_SET, meshset_for_bit_level);
769 
770  // subtract those elements which will be refined, i.e. disconnected from other
771  // side elements, and connected to new prisms, if they are created
772  meshset_3d_ents = subtract(meshset_3d_ents, side_ents3d);
773  CHKERR moab.add_entities(meshset_for_bit_level, meshset_3d_ents);
774 
775  // create new 3d ents on "father" side
776  Range new_3d_ents;
777  for (Range::iterator eit3d = side_ents3d.begin(); eit3d != side_ents3d.end();
778  eit3d++) {
779  auto miit_ref_ent = refined_ents_ptr->find(*eit3d);
780  if (miit_ref_ent == refined_ents_ptr->end())
781  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
782  "Tetrahedron not in database");
783 
784  int num_nodes;
785  const EntityHandle *conn;
786  CHKERR moab.get_connectivity(*eit3d, conn, num_nodes, true);
787  EntityHandle new_conn[num_nodes];
788  int nb_new_conn = 0;
789  for (int ii = 0; ii < num_nodes; ii++) {
790  std::map<EntityHandle, EntityHandle>::iterator mit =
791  map_nodes.find(conn[ii]);
792  if (mit != map_nodes.end()) {
793  new_conn[ii] = mit->second;
794  nb_new_conn++;
795  if (verb >= VERY_NOISY) {
796  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
797  new_conn[ii]);
798  }
799  } else {
800  new_conn[ii] = conn[ii];
801  }
802  }
803  if (nb_new_conn == 0) {
804  // Add this tet to bit ref level
805  CHKERR moab.add_entities(meshset_for_bit_level, &*eit3d, 1);
806  continue;
807  }
808 
809  // here is created new or prism is on interface
810  EntityHandle existing_ent = 0;
811  /* check if tet element with new connectivity is in database*/
812  auto child_iit =
813  refined_ents_ptr->get<Ent_Ent_mi_tag>().lower_bound(*eit3d);
814  auto hi_child_iit =
815  refined_ents_ptr->get<Ent_Ent_mi_tag>().upper_bound(*eit3d);
816 
817  // Check if child entity has the same connectivity
818  for (; child_iit != hi_child_iit; child_iit++) {
819  const EntityHandle *conn_ref_tet;
820  CHKERR moab.get_connectivity(child_iit->get()->getRefEnt(), conn_ref_tet,
821  num_nodes, true);
822  int nn = 0;
823  for (; nn < num_nodes; nn++) {
824  if (conn_ref_tet[nn] != new_conn[nn]) {
825  break;
826  }
827  }
828  if (nn == num_nodes) {
829  if (existing_ent != 0)
830  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
831  "Should be only one child entity with the same connectivity");
832  existing_ent = child_iit->get()->getRefEnt();
833  }
834  }
835 
836  switch (moab.type_from_handle(*eit3d)) {
837  case MBTET: {
838 
839  RefEntity_multiIndex::iterator child_it;
840  EntityHandle tet;
841  if (existing_ent == 0) {
842  Range new_conn_tet;
843  CHKERR moab.get_adjacencies(new_conn, 4, 3, false, new_conn_tet);
844  if (new_conn_tet.empty()) {
845  CHKERR moab.create_element(MBTET, new_conn, 4, tet);
846  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &tet, 1,
847  &*eit3d);
848 
849  } else {
850 
851  auto new_rit = refined_ents_ptr->get<Ent_mi_tag>().equal_range(
852  *new_conn_tet.begin());
853 
854  size_t nb_elems = std::distance(new_rit.first, new_rit.second);
855  if (nb_elems != 1)
856  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
857  "Can't find entity in database, size is %d", nb_elems);
858  tet = *new_conn_tet.begin();
859  }
860  } else {
861  tet = existing_ent;
862  }
863 
864  CHKERR moab.add_entities(meshset_for_bit_level, &tet, 1);
865  new_3d_ents.insert(tet);
866 
867  } break;
868  case MBPRISM: {
869  EntityHandle prism;
870  if (existing_ent == 0) {
871  Range new_conn_prism;
872  CHKERR moab.get_adjacencies(new_conn, 6, 3, false, new_conn_prism);
873  if (new_conn_prism.empty()) {
874  CHKERR moab.create_element(MBPRISM, new_conn, 6, prism);
875  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &prism, 1,
876  &*eit3d);
877  } else {
878  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
879  "It is prism with such connectivity, that case has to be "
880  "handled but this is not implemented");
881  }
882  } else {
883  prism = existing_ent;
884  }
885  CHKERR moab.add_entities(meshset_for_bit_level, &prism, 1);
886  new_3d_ents.insert(prism);
887  } break;
888  default:
889  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
890  "Not implemented element");
891  }
892  }
893 
894  struct SetParent {
895 
896  SetParent(MoFEM::Core &core) : cOre(core), mField(core) {}
897 
898  MoFEMErrorCode operator()(const EntityHandle ent, const EntityHandle parent,
899  const RefEntity_multiIndex *ref_ents_ptr) {
901  if (ent != parent) {
902  auto it = ref_ents_ptr->find(ent);
903  if (it != ref_ents_ptr->end()) {
904  parentsToChange[ent] = parent;
905  } else {
906  CHKERR mField.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
907  &ent, 1, &parent);
908  }
909  }
911  }
912 
913  MoFEMErrorCode override_parents(const RefEntity_multiIndex *ref_ents_ptr) {
915  for (auto &m : parentsToChange) {
916  auto it = ref_ents_ptr->find(m.first);
917  if (it != ref_ents_ptr->end()) {
918 
919  bool success = const_cast<RefEntity_multiIndex *>(ref_ents_ptr)
920  ->modify(it, RefEntity_change_parent(m.second));
921  if (!success)
922  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
923  "Impossible to set parent");
924  } else
925  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
926  "Entity not in database");
927  }
929  }
930 
931  private:
932  MoFEM::Core &cOre;
933  MoFEM::Interface &mField;
934  map<EntityHandle, EntityHandle> parentsToChange;
935  };
936 
937  SetParent set_parent(cOre);
938 
939  auto get_adj_ents = [&](const bool create) {
940  Range adj;
941  for (auto d : {1, 2}) {
942  // create new entities by adjacencies form new tets
943  CHKERR moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET), d, create,
944  adj, moab::Interface::UNION);
945  }
946  return adj;
947  };
948 
949  // Create entities
950  get_adj_ents(true);
951 
952  auto get_conn = [&](const auto e) {
953  int num_nodes;
954  const EntityHandle *conn;
955  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
956  return std::make_pair(conn, num_nodes);
957  };
958 
959  auto get_new_conn = [&](auto conn) {
960  std::array<EntityHandle, 8> new_conn;
961  int nb_new_conn = 0;
962  for (int ii = 0; ii != conn.second; ++ii) {
963  auto mit = map_nodes.find(conn.first[ii]);
964  if (mit != map_nodes.end()) {
965  new_conn[ii] = mit->second;
966  nb_new_conn++;
967  if (verb >= VERY_NOISY)
968  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn.first[ii],
969  new_conn[ii]);
970 
971  } else {
972  new_conn[ii] = conn.first[ii];
973  }
974  }
975  return std::make_pair(new_conn, nb_new_conn);
976  };
977 
978  auto get_reverse_conn = [&](auto conn) {
979  std::array<EntityHandle, 8> rev_conn;
980  int nb_new_conn = 0;
981  for (int ii = 0; ii != conn.second; ++ii) {
982  auto mit = reverse_map_nodes.find(conn.first[ii]);
983  if (mit != reverse_map_nodes.end()) {
984  rev_conn[ii] = mit->second;
985  nb_new_conn++;
986  if (verb >= VERY_NOISY)
987  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn.first[ii],
988  rev_conn[ii]);
989 
990  } else {
991  rev_conn[ii] = conn.first[ii];
992  }
993  }
994  return std::make_pair(rev_conn, nb_new_conn);
995  };
996 
997  auto get_new_ent = [&](auto new_conn, auto nb_nodes, int dim) {
998  Range new_ent;
999  CHKERR moab.get_adjacencies(&(new_conn.first[0]), nb_nodes, dim, false,
1000  new_ent);
1001  if (new_ent.size() != 1)
1002  THROW_MESSAGE("this entity should be in moab database");
1003  return new_ent.front();
1004  };
1005 
1006  auto create_prisms = [&]() {
1008 
1009  // Tags for setting side
1010  Tag th_interface_side;
1011  const int def_side[] = {0};
1012  CHKERR moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
1013  th_interface_side, MB_TAG_CREAT | MB_TAG_SPARSE,
1014  def_side);
1015 
1016 
1017  for (auto p = triangles.pair_begin(); p != triangles.pair_end(); ++p) {
1018  auto f = p->first;
1019  auto s = p->second;
1020 
1021  auto lo = refined_ents_ptr->lower_bound(f);
1022  auto hi = refined_ents_ptr->upper_bound(s);
1023  if (std::distance(lo, hi) != (s - f + 1))
1024  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1025  "Some triangles are not in database");
1026 
1027  for (; f <= s; ++f) {
1028 
1029  auto conn = get_conn(f);
1030  auto new_conn = get_new_conn(conn);
1031 
1032  if (new_conn.second) {
1033 
1034  auto set_side_tag = [&](auto new_triangle) {
1035  int new_side = 1;
1036  CHKERR moab.tag_set_data(th_interface_side, &new_triangle, 1,
1037  &new_side);
1038  };
1039 
1040  auto get_ent3d = [&](auto e) {
1041  Range ents_3d;
1042  CHKERR moab.get_adjacencies(&e, 1, 3, false, ents_3d);
1043  ents_3d = intersect(ents_3d, side_ents3d);
1044 
1045  switch (ents_3d.size()) {
1046  case 0:
1047  THROW_MESSAGE(
1048  "Did not find adjacent tets on one side of the interface, "
1049  "check its definition and try creating separate sidesets for "
1050  "each surface");
1051  case 2:
1052  THROW_MESSAGE(
1053  "Found both adjacent tets on one side of the interface, "
1054  "check "
1055  "its "
1056  "definition and try creating separate sidesets for each "
1057  "surface");
1058  default:
1059  break;
1060  }
1061 
1062  return ents_3d.front();
1063  };
1064 
1065  auto get_sense = [&](auto e, auto ent3d) {
1066  int sense, side, offset;
1067  CHKERR moab.side_number(ent3d, e, side, sense, offset);
1068  if (sense != 1 && sense != -1) {
1069  SETERRQ(
1071  "Undefined sense of a trinagle on the interface, check its "
1072  "definition and try creating separate sidesets for each "
1073  "surface");
1074  }
1075  return sense;
1076  };
1077 
1078  auto new_triangle = get_new_ent(new_conn, 3, 2);
1079  set_side_tag(new_triangle);
1080 
1081  if (add_interface_entities) {
1082 
1083  if (inhered_from_bit_level.any())
1084  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1085  "not implemented for inhered_from_bit_level");
1086 
1087  // set prism connectivity
1088  EntityHandle prism_conn[6] = {
1089  conn.first[0], conn.first[1], conn.first[2],
1090 
1091  new_conn.first[0], new_conn.first[1], new_conn.first[2]};
1092  if (get_sense(f, get_ent3d(f)) == -1) {
1093  // swap nodes in triangles for correct prism creation
1094  std::swap(prism_conn[1], prism_conn[2]);
1095  std::swap(prism_conn[4], prism_conn[5]);
1096  }
1097 
1098  EntityHandle prism;
1099  CHKERR moab.create_element(MBPRISM, prism_conn, 6, prism);
1100  CHKERR moab.add_entities(meshset_for_bit_level, &prism, 1);
1101  }
1102  }
1103  }
1104  }
1105 
1107  };
1108 
1109  auto set_parnets = [&](auto side_adj_faces_and_edges) {
1111 
1112  for (auto p = side_adj_faces_and_edges.pair_begin();
1113  p != side_adj_faces_and_edges.pair_end(); ++p) {
1114  auto f = p->first;
1115  auto s = p->second;
1116 
1117  for (; f <= s; ++f) {
1118  auto conn = get_conn(f);
1119  auto rev_conn = get_reverse_conn(conn);
1120  if (rev_conn.second) {
1121  auto rev_ent = get_new_ent(rev_conn, conn.second, conn.second - 1);
1122  CHKERR set_parent(f, rev_ent, refined_ents_ptr);
1123 
1124  }
1125  }
1126  };
1127 
1129  };
1130 
1131  auto all_new_adj_entities = [&](const bool create) {
1132  Range adj;
1133  for (auto d : {1, 2})
1134  CHKERR moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET), d, create,
1135  adj, moab::Interface::UNION);
1136  return adj;
1137  };
1138 
1139  auto add_new_prisms_which_parents_are_part_of_other_intefaces = [&]() {
1141 
1142  Tag th_interface_side;
1143  CHKERR moab.tag_get_handle("INTERFACE_SIDE", th_interface_side);
1144 
1145  Range new_3d_prims = new_3d_ents.subset_by_type(MBPRISM);
1146  for (Range::iterator pit = new_3d_prims.begin(); pit != new_3d_prims.end();
1147  ++pit) {
1148 
1149  // get parent entity
1150  EntityHandle parent_prism;
1151  CHKERR moab.tag_get_data(cOre.get_th_RefParentHandle(), &*pit, 1,
1152  &parent_prism);
1153  if (moab.type_from_handle(parent_prism) != MBPRISM)
1154  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1155  "this prism should have parent which is prism as well");
1156 
1157  int num_nodes;
1158  // parent prism
1159  const EntityHandle *conn_parent;
1160  CHKERR moab.get_connectivity(parent_prism, conn_parent, num_nodes, true);
1161  Range face_side3_parent, face_side4_parent;
1162  CHKERR moab.get_adjacencies(conn_parent, 3, 2, false, face_side3_parent);
1163  CHKERR moab.get_adjacencies(&conn_parent[3], 3, 2, false,
1164  face_side4_parent);
1165  if (face_side3_parent.size() != 1)
1166  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1167  "parent face3.size() = %u", face_side3_parent.size());
1168 
1169  if (face_side4_parent.size() != 1)
1170  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1171  "parent face4.size() = %u", face_side4_parent.size());
1172 
1173  // new prism
1174  const EntityHandle *conn;
1175  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
1176  Range face_side3, face_side4;
1177  CHKERR moab.get_adjacencies(conn, 3, 2, false, face_side3);
1178  CHKERR moab.get_adjacencies(&conn[3], 3, 2, false, face_side4);
1179  if (face_side3.size() != 1)
1180  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1181  "face3 is missing");
1182 
1183  if (face_side4.size() != 1)
1184  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1185  "face4 is missing");
1186 
1187  std::vector<EntityHandle> face(2), parent_face(2);
1188  face[0] = *face_side3.begin();
1189  face[1] = *face_side4.begin();
1190  parent_face[0] = *face_side3_parent.begin();
1191  parent_face[1] = *face_side4_parent.begin();
1192  for (int ff = 0; ff != 2; ++ff) {
1193  if (parent_face[ff] == face[ff])
1194  continue;
1195  int interface_side;
1196  CHKERR moab.tag_get_data(th_interface_side, &parent_face[ff], 1,
1197  &interface_side);
1198  CHKERR moab.tag_set_data(th_interface_side, &face[ff], 1,
1199  &interface_side);
1200  EntityHandle parent_tri;
1201  CHKERR moab.tag_get_data(cOre.get_th_RefParentHandle(), &face[ff], 1,
1202  &parent_tri);
1203  if (parent_tri != parent_face[ff]) {
1204  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1205  "Wrong parent %lu", parent_tri);
1206  }
1207  }
1208  }
1210  };
1211 
1212  CHKERR create_prisms();
1213 
1214  CHKERR set_parnets(all_new_adj_entities(true));
1215  CHKERR add_new_prisms_which_parents_are_part_of_other_intefaces();
1216 
1217  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
1218  meshset_for_bit_level, 3, bit);
1219  CHKERR moab.delete_entities(&meshset_for_bit_level, 1);
1220  CHKERR moab.clear_meshset(&children[0], 3);
1221 
1222  auto reconstruct_refined_ents = [&]() {
1226  };
1227 
1228  // Finalise by adding new tets and prism ti bit level
1229  CHKERR set_parent.override_parents(refined_ents_ptr);
1230 
1231  // Add function which reconstruct core multi-index. Node merging is messy
1232  // process and entity parent could be changed without notification to
1233  // multi-index. TODO Issue has to be tracked down better, however in principle
1234  // is better not to modify multi-index each time parent is changed, that makes
1235  // code slow. Is better to do it in the bulk as below.
1236  CHKERR reconstruct_refined_ents();
1237 
1239 }
Deprecated interface functions.
Tag get_th_RefParentHandle() const
Definition: Core.hpp:150
virtual moab::Interface & get_moab()=0
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
virtual const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:626
Core (interface) class.
Definition: Core.hpp:50
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
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:657
const int dim
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
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:693
#define CHKERR
Inline error check.
Definition: definitions.h:602
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
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
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
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: