v0.14.0
PrismInterface.cpp
Go to the documentation of this file.
1 /** \file PrismInterface.cpp
2  * \brief Insert prisms in the interface between two surfaces
3  * \todo FIXME this is not so good implementation
4  *
5  * \ingroup mofem_prism_interface
6  */
7 
8 namespace MoFEM {
9 
11 PrismInterface::query_interface(boost::typeindex::type_index type_index,
12  UnknownInterface **iface) const {
13  *iface = const_cast<PrismInterface *>(this);
14  return 0;
15 }
16 
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 }
30 
32  const CubitBCType cubit_bc_type,
33  const BitRefLevel mesh_bit_level,
34  Range seed_side, const bool recursive,
35  int verb) {
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()
45  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
46  if (miit != meshsets_manager_ptr->getMeshsetsMultindex()
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 }
55 
57  const CubitBCType cubit_bc_type,
58  const BitRefLevel mesh_bit_level,
59  const bool recursive, int verb) {
60  return getSides(msId, cubit_bc_type, mesh_bit_level, Range(), recursive,
61  verb);
62 }
63 
65  const BitRefLevel mesh_bit_level,
66  Range seed_side, const bool recursive,
67  int verb) {
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 }
386 
388  const BitRefLevel mesh_bit_level,
389  const bool recursive, int verb) {
390  return getSides(sideset, mesh_bit_level, Range(), recursive, verb);
391 }
392 
394  const EntityHandle sideset, const BitRefLevel mesh_bit_level,
395  const bool recursive, Range &faces_with_three_nodes_on_front, int verb) {
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
518 
520  const BitRefLevel &bit,
521  const int msId,
522  const CubitBCType cubit_bc_type,
523  const bool add_interface_entities,
524  const bool recursive, int verb) {
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()
534  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
535  if (miit != meshsets_manager_ptr->getMeshsetsMultindex()
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 }
545 
547  const BitRefLevel &bit,
548  const EntityHandle sideset,
549  const bool add_interface_entities,
550  const bool recursive, int verb) {
551 
553  CHKERR splitSides(meshset, bit, BitRefLevel(), BitRefLevel(), sideset,
554  add_interface_entities, recursive, verb);
556 }
557 
559  const EntityHandle meshset, const BitRefLevel &bit,
560  const BitRefLevel &inhered_from_bit_level,
561  const BitRefLevel &inhered_from_bit_level_mask, const EntityHandle sideset,
562  const bool add_interface_entities, const bool recursive, int verb) {
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 }
1220 } // namespace MoFEM
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::Ent_mi_tag
Definition: TagMultiIndices.hpp:21
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::PrismInterface
Create interface from given surface and insert flat prisms in-between.
Definition: PrismInterface.hpp:23
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::Ent_Ent_mi_tag
Definition: TagMultiIndices.hpp:55
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
MoFEM::Composite_EntType_and_ParentEntType_mi_tag
Definition: TagMultiIndices.hpp:77
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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::RefEntExtractor
Extract entity handle form multi-index container.
Definition: Templates.hpp:1769
MoFEM::Composite_Cubit_msId_And_MeshsetType_mi_tag
Definition: TagMultiIndices.hpp:15
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
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
Range
MoFEM::MeshsetsManager::getMeshsetsMultindex
CubitMeshSet_multiIndex & getMeshsetsMultindex()
Definition: MeshsetsManager.hpp:229
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::Types::CubitBCType
std::bitset< 32 > CubitBCType
Definition: Types.hpp:52
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
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::RefEntity_change_parent
change parent
Definition: RefEntsMultiIndices.hpp:812
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
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::PrismInterface::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: PrismInterface.cpp:11
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
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
MoFEM::PrismInterface::findFacesWithThreeNodesOnInternalSurfaceSkin
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.
Definition: PrismInterface.cpp:393
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