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

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

#include <src/interfaces/PrismInterface.hpp>

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

Public Member Functions

MoFEMErrorCode query_interface (const MOFEMuuid &uuid, UnknownInterface **iface) const
 
 PrismInterface (const MoFEM::Core &core)
 
 ~PrismInterface ()=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 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
continuity_check_on_contact_prism_side_ele.cpp, 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)) {
38 
39  if (!LogManager::checkIfChannelExist("PRISM_INTERFACE")) {
40  auto core_log = logging::core::get();
41  core_log->add_sink(
42  LogManager::createSink(LogManager::getStrmWorld(), "PRISM_INTERFACE"));
43  LogManager::setLog("PRISM_INTERFACE");
44  MOFEM_LOG_TAG("PRISM_INTERFACE", "PrismInterface");
45  }
46 
47  MOFEM_LOG("PRISM_INTERFACE", Sev::noisy) << "Prism interface created";
48 }

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

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

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

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

◆ getSides() [2/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
continuity_check_on_contact_prism_side_ele.cpp, 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 50 of file PrismInterface.cpp.

53  {
54 
55  Interface &m_field = cOre;
56  MeshsetsManager *meshsets_manager_ptr;
58  CHKERR m_field.getInterface(meshsets_manager_ptr);
59  CubitMeshSet_multiIndex::index<
60  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator miit =
61  meshsets_manager_ptr->getMeshsetsMultindex()
62  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
63  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
64  if (miit != meshsets_manager_ptr->getMeshsetsMultindex()
65  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
66  .end()) {
67  CHKERR getSides(miit->meshset, mesh_bit_level, recursive, verb);
68  } else {
69  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "msId is not there");
70  }
72 }

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

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

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

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

551  {
552 
554  CHKERR splitSides(meshset, bit, BitRefLevel(), BitRefLevel(), sideset,
555  add_interface_entities, recursive, verb);
557 }

◆ 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, split_sideset.cpp, and test_jacobian_of_simple_contact_element.cpp.

Definition at line 520 of file PrismInterface.cpp.

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

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:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
EntityHandle
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
NOISY
@ NOISY
Definition: definitions.h:280
n
static Index< 'n', 3 > n
Definition: BasicFeTools.hpp:78
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:50
FTensor::d
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
MoFEM::reconstructMultiIndex
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:1116
MoFEM::PrismInterface::cOre
MoFEM::Core & cOre
Definition: PrismInterface.hpp:41
VERY_NOISY
@ VERY_NOISY
Definition: definitions.h:281
dim
const int dim
Definition: elec_phys_2D.cpp:12
MoFEM::IDD_MOFEMPrismInterface
static const MOFEMuuid IDD_MOFEMPrismInterface
Definition: PrismInterface.hpp:28
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
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:762
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:315
m
static Index< 'm', 3 > m
Definition: BasicFeTools.hpp:77
MoFEM::CoreTmp< 0 >::get_th_RefParentHandle
Tag get_th_RefParentHandle() const
Definition: Core.hpp:192
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
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:1080
HenkyOps::f
auto f
Definition: HenkyOps.hpp:19
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:779
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:283
MoFEM::LogManager::checkIfChannelExist
static bool checkIfChannelExist(const std::string channel)
Check if channel exist.
Definition: LogManager.cpp:378
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:126
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:343
THROW_MESSAGE
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:125
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:77
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:312
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:80
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:327
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:123
MoFEM::PrismInterface::PrismInterface
PrismInterface(const MoFEM::Core &core)
Definition: PrismInterface.cpp:36
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
QUIET
@ QUIET
Definition: definitions.h:277
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:368
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1129
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:520
VERY_VERBOSE
@ VERY_VERBOSE
Definition: definitions.h:279
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Definition: UnknownInterface.hpp:130
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:33
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