v0.14.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 (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 
 PrismInterface (const MoFEM::Core &core)
 
 ~PrismInterface ()=default
 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 int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, Range seed_side, 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 getSides (const EntityHandle sideset, const BitRefLevel mesh_bit_level, Range seed_side, 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 (bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface reference 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
 

Public Attributes

MoFEM::CorecOre
 

Additional Inherited Members

- Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode getLibVersion (Version &version)
 Get library version. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

Detailed Description

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

Examples
continuity_check_on_contact_prism_side_ele.cpp, forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, simple_contact.cpp, simple_contact_thermal.cpp, split_sideset.cpp, and test_jacobian_of_simple_contact_element.cpp.

Definition at line 23 of file PrismInterface.hpp.

Constructor & Destructor Documentation

◆ PrismInterface()

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

Definition at line 17 of file PrismInterface.cpp.

18  : cOre(const_cast<Core &>(core)) {
19 
20  if (!LogManager::checkIfChannelExist("PRISM_INTERFACE")) {
21  auto core_log = logging::core::get();
22  core_log->add_sink(
23  LogManager::createSink(LogManager::getStrmWorld(), "PRISM_INTERFACE"));
24  LogManager::setLog("PRISM_INTERFACE");
25  MOFEM_LOG_TAG("PRISM_INTERFACE", "PrismInterface");
26  }
27 
28  MOFEM_LOG("PRISM_INTERFACE", Sev::noisy) << "Prism interface created";
29 }

◆ ~PrismInterface()

MoFEM::PrismInterface::~PrismInterface ( )
default

destructor

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

395  {
396  Interface &m_field = cOre;
397  moab::Interface &moab = m_field.get_moab();
399 
400  Range mesh_level_ents3d;
401  Range mesh_level_edges, mesh_level_tris;
402  if (mesh_bit_level.any()) {
403  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
404  mesh_bit_level, BitRefLevel().set(), MBTET, mesh_level_ents3d);
405 
406  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
407  mesh_bit_level, BitRefLevel().set(), MBTRI, mesh_level_tris);
408 
409  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
410  mesh_bit_level, BitRefLevel().set(), MBEDGE, mesh_level_edges);
411  }
412 
413  Skinner skin(&moab);
414  // get interface triangles from side set
415  Range triangles;
416  CHKERR moab.get_entities_by_type(sideset, MBTRI, triangles, recursive);
417 
418  if (mesh_bit_level.any()) {
419  triangles = intersect(triangles, mesh_level_tris);
420  }
421 
422  MOFEM_LOG_C("PRISM_INTERFACE", Sev::verbose, "Nb. of triangles in set %u",
423  triangles.size());
424 
425  // get nodes, edges and 3d ents (i.e. tets and prisms)
426  Range nodes; // nodes from triangles
427  CHKERR moab.get_connectivity(triangles, nodes, true);
428 
429  Range ents3d; // 3d ents from nodes
430  CHKERR moab.get_adjacencies(nodes, 3, false, ents3d, moab::Interface::UNION);
431 
432  if (mesh_bit_level.any()) {
433  ents3d = intersect(ents3d, mesh_level_ents3d);
434  }
435  // take skin faces
436  Range skin_faces; // skin faces from 3d ents
437  CHKERR skin.find_skin(0, ents3d.subset_by_type(MBTET), false, skin_faces);
438 
439  // take skin edges (boundary of surface if there is any)
440  Range skin_edges_boundary; // skin edges from triangles
441  CHKERR skin.find_skin(0, triangles, false, skin_edges_boundary);
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 
452  // note: that skin faces edges do not contain internal boundary
453  // note: that prisms are not included in ents3d, so if ents3d have border with
454  // other inteface is like external boundary skin edges boundary are internal
455  // edge <- skin_faces_edges contains edges which are on the body boundary <-
456  // that is the trick
457  skin_edges_boundary =
458  subtract(skin_edges_boundary,
459  skin_faces_edges); // from skin edges subtract edges from skin
460  // faces of 3d ents (only internal edges)
461 
462  if (verb >= VERY_VERBOSE) {
463  EntityHandle out_meshset;
464  CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
465  CHKERR moab.add_entities(out_meshset, triangles);
466  CHKERR moab.write_file("triangles.vtk", "VTK", "", &out_meshset, 1);
467  CHKERR moab.delete_entities(&out_meshset, 1);
468  CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
469  CHKERR moab.add_entities(out_meshset, ents3d);
470  CHKERR moab.write_file("ents3d.vtk", "VTK", "", &out_meshset, 1);
471  CHKERR moab.delete_entities(&out_meshset, 1);
472  CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
473  CHKERR moab.add_entities(out_meshset, skin_edges_boundary);
474  CHKERR moab.write_file("skin_edges_boundary.vtk", "VTK", "", &out_meshset,
475  1);
476  CHKERR moab.delete_entities(&out_meshset, 1);
477  }
478 
479  // Get nodes on boundary edge
480  Range skin_nodes_boundary;
481  CHKERR moab.get_connectivity(skin_edges_boundary, skin_nodes_boundary, true);
482 
483  // Remove node which are boundary with other existing interface
484  Range prisms_nodes;
485  CHKERR
486  moab.get_connectivity(ents3d.subset_by_type(MBPRISM), prisms_nodes, true);
487 
488  skin_nodes_boundary = subtract(skin_nodes_boundary, prisms_nodes);
489 
490  // use nodes on body boundary and interface (without internal boundary) to
491  // find adjacent tets
492  Range nodes_without_front = subtract(
493  nodes, skin_nodes_boundary); // nodes_without_front adjacent to all split
494  // face edges except those on internal edge
495 
496  Range skin_nodes_boundary_tris;
497  CHKERR moab.get_adjacencies(skin_nodes_boundary, 2, false,
498  skin_nodes_boundary_tris, moab::Interface::UNION);
499  Range skin_nodes_boundary_tris_nodes;
500  CHKERR moab.get_connectivity(skin_nodes_boundary_tris,
501  skin_nodes_boundary_tris_nodes, true);
502 
503  skin_nodes_boundary_tris_nodes =
504  subtract(skin_nodes_boundary_tris_nodes, skin_nodes_boundary);
505  Range skin_nodes_boundary_tris_nodes_tris;
506  CHKERR moab.get_adjacencies(skin_nodes_boundary_tris_nodes, 2, false,
507  skin_nodes_boundary_tris_nodes_tris,
508  moab::Interface::UNION);
509 
510  // Triangle which has tree nodes on front boundary
511  skin_nodes_boundary_tris =
512  intersect(triangles, subtract(skin_nodes_boundary_tris,
513  skin_nodes_boundary_tris_nodes_tris));
514  faces_with_three_nodes_on_front.swap(skin_nodes_boundary_tris);
515 
517 } // namespace MoFEM

◆ getSides() [1/4]

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

389  {
390  return getSides(sideset, mesh_bit_level, Range(), recursive, verb);
391 }

◆ getSides() [2/4]

MoFEMErrorCode MoFEM::PrismInterface::getSides ( const EntityHandle  sideset,
const BitRefLevel  mesh_bit_level,
Range  seed_side,
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
seed_sideuse seed to decide which side to choose first
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 64 of file PrismInterface.cpp.

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

◆ getSides() [3/4]

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
continuity_check_on_contact_prism_side_ele.cpp, forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, simple_contact.cpp, and test_jacobian_of_simple_contact_element.cpp.

Definition at line 56 of file PrismInterface.cpp.

59  {
60  return getSides(msId, cubit_bc_type, mesh_bit_level, Range(), recursive,
61  verb);
62 }

◆ getSides() [4/4]

MoFEMErrorCode MoFEM::PrismInterface::getSides ( const int  msId,
const CubitBCType  cubit_bc_type,
const BitRefLevel  mesh_bit_level,
Range  seed_side,
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
seed_sideuse seed to decide which side to choose first
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
if bit_level == BitRefLevel.set() then interface will be added on all bit levels

Definition at line 31 of file PrismInterface.cpp.

35  {
36 
37  Interface &m_field = cOre;
38  MeshsetsManager *meshsets_manager_ptr;
40  CHKERR m_field.getInterface(meshsets_manager_ptr);
41  CubitMeshSet_multiIndex::index<
42  Composite_Cubit_msId_And_MeshsetType_mi_tag>::type::iterator miit =
43  meshsets_manager_ptr->getMeshsetsMultindex()
44  .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
45  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
46  if (miit != meshsets_manager_ptr->getMeshsetsMultindex()
47  .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
48  .end()) {
49  CHKERR getSides(miit->meshset, mesh_bit_level, recursive, verb);
50  } else {
51  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "msId is not there");
52  }
54 }

◆ query_interface()

MoFEMErrorCode MoFEM::PrismInterface::query_interface ( boost::typeindex::type_index  type_index,
UnknownInterface **  iface 
) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 11 of file PrismInterface.cpp.

12  {
13  *iface = const_cast<PrismInterface *>(this);
14  return 0;
15 }

◆ splitSides() [1/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 558 of file PrismInterface.cpp.

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

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

550  {
551 
553  CHKERR splitSides(meshset, bit, BitRefLevel(), BitRefLevel(), sideset,
554  add_interface_entities, recursive, verb);
556 }

◆ splitSides() [3/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
continuity_check_on_contact_prism_side_ele.cpp, forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, simple_contact.cpp, and test_jacobian_of_simple_contact_element.cpp.

Definition at line 519 of file PrismInterface.cpp.

524  {
525 
526  Interface &m_field = cOre;
527  MeshsetsManager *meshsets_manager_ptr;
529  CHKERR m_field.getInterface(meshsets_manager_ptr);
530  CubitMeshSet_multiIndex::index<
531  Composite_Cubit_msId_And_MeshsetType_mi_tag>::type::iterator miit =
532  meshsets_manager_ptr->getMeshsetsMultindex()
533  .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
534  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
535  if (miit != meshsets_manager_ptr->getMeshsetsMultindex()
536  .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
537  .end()) {
538  CHKERR splitSides(meshset, bit, miit->meshset, add_interface_entities,
539  recursive, verb);
540  } else {
541  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "msId is not there");
542  }
544 }

Member Data Documentation

◆ cOre

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

Definition at line 28 of file PrismInterface.hpp.


The documentation for this struct was generated from the following files:
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::LogManager::checkIfChannelExist
static bool checkIfChannelExist(const std::string channel)
Check if channel exist.
Definition: LogManager.cpp:404
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
EntityHandle
NOISY
@ NOISY
Definition: definitions.h:211
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VERY_NOISY
@ VERY_NOISY
Definition: definitions.h:212
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
MoFEM::RefEntity_multiIndex
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity, 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, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
Definition: RefEntsMultiIndices.hpp:760
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
sdf.r
int r
Definition: sdf.py:8
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::reconstructMultiIndex
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:1831
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::RefEntity_multiIndex_view_by_hashed_parent_entity
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::getEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >> > > RefEntity_multiIndex_view_by_hashed_parent_entity
multi-index view of RefEntity by parent entity
Definition: RefEntsMultiIndices.hpp:777
MoFEM::insertOrdered
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:1795
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::CoreTmp< 0 >::get_th_RefParentHandle
Tag get_th_RefParentHandle() const
Definition: Core.hpp:197
MoFEM::PrismInterface::splitSides
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...
Definition: PrismInterface.cpp:519
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1876
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
convert.n
n
Definition: convert.py:82
Range
MoFEM::PrismInterface::cOre
MoFEM::Core & cOre
Definition: PrismInterface.hpp:28
FTensor::dd
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
save_range
auto save_range
Definition: thermo_elastic.cpp:144
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
MoFEM::PrismInterface::getSides
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...
Definition: PrismInterface.cpp:56
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::CoreInterface::get_ref_ents
virtual const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
QUIET
@ QUIET
Definition: definitions.h:208
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEM::PrismInterface::PrismInterface
PrismInterface(const MoFEM::Core &core)
Definition: PrismInterface.cpp:17
VERY_VERBOSE
@ VERY_VERBOSE
Definition: definitions.h:210
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MoFEM::BcManagerImplTools::get_adj_ents
auto get_adj_ents(moab::Interface &moab, const Range &ents)
Definition: BcManager.cpp:17