v0.13.1
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
8namespace MoFEM {
9
11PrismInterface::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(
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);
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
64MoFEMErrorCode PrismInterface::getSides(const EntityHandle sideset,
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
387MoFEMErrorCode PrismInterface::getSides(const EntityHandle sideset,
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);
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 {
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:
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:
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:
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) {
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
static Index< 'p', 3 > p
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
constexpr double a
@ VERY_VERBOSE
Definition: definitions.h:210
@ QUIET
Definition: definitions.h:208
@ NOISY
Definition: definitions.h:211
@ VERY_NOISY
Definition: definitions.h:212
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
const int dim
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
auto save_range
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
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
virtual const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
auto bit
set bit
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
auto f
Definition: HenckyOps.hpp:5
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< 32 > CubitBCType
Definition: Types.hpp:52
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1479
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:1438
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:1402
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double r
rate factor
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:82
Tag get_th_RefParentHandle() const
Definition: Core.hpp:197
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
static bool checkIfChannelExist(const std::string channel)
Check if channel exist.
Definition: LogManager.cpp:374
Interface for managing meshsets containing materials and boundary conditions.
CubitMeshSet_multiIndex & getMeshsetsMultindex()
Create interface from given surface and insert flat prisms in-between.
PrismInterface(const MoFEM::Core &core)
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.
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
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...
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...
Extract entity handle form multi-index container.
Definition: Templates.hpp:1378
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.