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

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

#include <src/interfaces/PrismInterface.hpp>

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

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 
 PrismInterface (const MoFEM::Core &core)
 
 ~PrismInterface ()=default
 destructor More...
 
MoFEMErrorCode getSides (const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
 Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshset. More...
 
MoFEMErrorCode getSides (const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, Range seed_side, const bool recursive, int verb=QUIET)
 Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshset. More...
 
MoFEMErrorCode getSides (const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
 Store tetrahedra from each side of the interface in two child meshsets of the parent meshset. More...
 
MoFEMErrorCode getSides (const EntityHandle sideset, const BitRefLevel mesh_bit_level, Range seed_side, const bool recursive, int verb=QUIET)
 Store tetrahedra from each side of the interface in two child meshsets of the parent meshset. More...
 
MoFEMErrorCode findFacesWithThreeNodesOnInternalSurfaceSkin (const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, Range &faces_with_three_nodes_on_front, int verb=QUIET)
 Find triangles which have three nodes on internal surface skin. More...
 
MoFEMErrorCode splitSides (const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
 Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in-between. More...
 
MoFEMErrorCode splitSides (const EntityHandle meshset, const BitRefLevel &bit, const EntityHandle sideset, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
 Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in-between. More...
 
MoFEMErrorCode splitSides (const EntityHandle meshset, const BitRefLevel &bit, const BitRefLevel &inhered_from_bit_level, const BitRefLevel &inhered_from_bit_level_mask, const EntityHandle sideset, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
 Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in-between. More...
 
- Public Member Functions inherited from MoFEM::UnknownInterface
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

Public Attributes

MoFEM::CorecOre
 

Additional Inherited Members

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

Detailed Description

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

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

Definition at line 23 of file PrismInterface.hpp.

Constructor & Destructor Documentation

◆ PrismInterface()

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

Definition at line 17 of file PrismInterface.cpp.

18 : cOre(const_cast<Core &>(core)) {
19
20 if (!LogManager::checkIfChannelExist("PRISM_INTERFACE")) {
21 auto core_log = logging::core::get();
22 core_log->add_sink(
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}
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
CoreTmp< 0 > Core
Definition: Core.hpp:1086
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

◆ ~PrismInterface()

MoFEM::PrismInterface::~PrismInterface ( )
default

destructor

Member Function Documentation

◆ findFacesWithThreeNodesOnInternalSurfaceSkin()

MoFEMErrorCode MoFEM::PrismInterface::findFacesWithThreeNodesOnInternalSurfaceSkin ( const EntityHandle  sideset,
const BitRefLevel  mesh_bit_level,
const bool  recursive,
Range faces_with_three_nodes_on_front,
int  verb = QUIET 
)

Find triangles which have three nodes on internal surface skin.

Internal surface skin is a set of all edges on the boundary of a given surface inside the body. This set of edges is also called the surface front. If a triangle has three nodes on the surface front, none of these nodes can be split. Therefore, such a triangle cannot be split and should be removed from the surface.

Parameters
sidesetmeshset with surface
mesh_bit_levelbit ref level of the volume mesh
recursiveif true search in sub-meshsets
faces_with_three_nodes_on_frontreturned faces
verbverbosity level
Returns
error code

Definition at line 393 of file PrismInterface.cpp.

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

◆ getSides() [1/4]

MoFEMErrorCode MoFEM::PrismInterface::getSides ( const EntityHandle  sideset,
const BitRefLevel  mesh_bit_level,
const bool  recursive,
int  verb = QUIET 
)

Store tetrahedra from each side of the interface in two child meshsets of the parent meshset.

Additional third child meshset contains nodes which can be split and skin edges

Parameters
sidesetparent meshset with the surface
mesh_bit_levelinterface is added on this bit level
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
if bit_level == BitRefLevel.set() then interface will be added on all bit levels
  1. Get tets adjacent to nodes of the interface meshset.
  2. Take skin faces from these tets and get edges from that skin.
  3. Take skin from triangles of the interface.
  4. Subtract edges of skin faces from skin of triangles in order to get edges in the volume of the body, and not on the interface boundary.
  5. Iterate between all triangles of the interface and find adjacent tets on each side of the interface

Definition at line 387 of file PrismInterface.cpp.

389 {
390 return getSides(sideset, mesh_bit_level, Range(), recursive, verb);
391}
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...

◆ getSides() [2/4]

MoFEMErrorCode MoFEM::PrismInterface::getSides ( const EntityHandle  sideset,
const BitRefLevel  mesh_bit_level,
Range  seed_side,
const bool  recursive,
int  verb = QUIET 
)

Store tetrahedra from each side of the interface in two child meshsets of the parent meshset.

Additional third child meshset contains nodes which can be split and skin edges

Parameters
sidesetparent meshset with the surface
mesh_bit_levelinterface is added on this bit level
seed_sideuse seed to decide which side to choose first
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
if bit_level == BitRefLevel.set() then interface will be added on all bit levels
  1. Get tets adjacent to nodes of the interface meshset.
  2. Take skin faces from these tets and get edges from that skin.
  3. Take skin from triangles of the interface.
  4. Subtract edges of skin faces from skin of triangles in order to get edges in the volume of the body, and not on the interface boundary.
  5. Iterate between all triangles of the interface and find adjacent tets on each side of the interface

Definition at line 64 of file PrismInterface.cpp.

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

◆ getSides() [3/4]

MoFEMErrorCode MoFEM::PrismInterface::getSides ( const int  msId,
const CubitBCType  cubit_bc_type,
const BitRefLevel  mesh_bit_level,
const bool  recursive,
int  verb = QUIET 
)

Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshset.

Additional third child meshset contains nodes which can be split and skin edges

Parameters
msIdId of meshset
cubit_bc_typetype of meshset (NODESET, SIDESET or BLOCKSET and more)
mesh_bit_levelinterface is added on this bit level
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
if bit_level == BitRefLevel.set() then interface will be added on all bit levels
Examples
continuity_check_on_contact_prism_side_ele.cpp, forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, simple_contact.cpp, and test_jacobian_of_simple_contact_element.cpp.

Definition at line 56 of file PrismInterface.cpp.

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

◆ getSides() [4/4]

MoFEMErrorCode MoFEM::PrismInterface::getSides ( const int  msId,
const CubitBCType  cubit_bc_type,
const BitRefLevel  mesh_bit_level,
Range  seed_side,
const bool  recursive,
int  verb = QUIET 
)

Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshset.

Additional third child meshset contains nodes which can be split and skin edges

Parameters
msIdId of meshset
cubit_bc_typetype of meshset (NODESET, SIDESET or BLOCKSET and more)
mesh_bit_levelinterface is added on this bit level
seed_sideuse seed to decide which side to choose first
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
if bit_level == BitRefLevel.set() then interface will be added on all bit levels

Definition at line 31 of file PrismInterface.cpp.

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

◆ query_interface()

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

Implements MoFEM::UnknownInterface.

Definition at line 11 of file PrismInterface.cpp.

12 {
13 *iface = const_cast<PrismInterface *>(this);
14 return 0;
15}
PrismInterface(const MoFEM::Core &core)

◆ splitSides() [1/3]

MoFEMErrorCode MoFEM::PrismInterface::splitSides ( const EntityHandle  meshset,
const BitRefLevel bit,
const BitRefLevel inhered_from_bit_level,
const BitRefLevel inhered_from_bit_level_mask,
const EntityHandle  sideset,
const bool  add_interface_entities,
const bool  recursive = false,
int  verb = QUIET 
)

Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in-between.

Parameters
meshsetvolume meshset containing 3D entities around the interface
bitbit ref level on which new entities will be stored
inhered_from_bit_levelinherit nodes and other entities form this bit level
inhered_from_bit_level_maskcorresponding mask
sidesetmeshset with surface
add_interface_entitiesif true add prism elements at interface
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
Parent meshset must have three child meshsets: two with tetrahedra from each side of the interface, third containing nodes which can be split and skin edges
inhered_from_bit_level needs to be specified for some meshsets with interfaces. Some nodes on some refinement levels are dividing edges but not splitting faces. Inheriting those nodes will not split faces.

Definition at line 558 of file PrismInterface.cpp.

562 {
563 Interface &m_field = cOre;
564 moab::Interface &moab = m_field.get_moab();
565 auto refined_ents_ptr = m_field.get_ref_ents();
567
568 std::vector<EntityHandle> children;
569 // get children meshsets
570 CHKERR moab.get_child_meshsets(sideset, children);
571 if (children.size() != 3)
572 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
573 "should be 3 child meshsets, each of them contains tets on two "
574 "sides of interface");
575
576 // 3d ents on "father" side
577 Range side_ents3d;
578 CHKERR moab.get_entities_by_handle(children[0], side_ents3d, false);
579 // 3d ents on "mother" side
580 Range other_ents3d;
581 CHKERR moab.get_entities_by_handle(children[1], other_ents3d, false);
582 // faces of interface
583 Range triangles;
584 CHKERR moab.get_entities_by_type(sideset, MBTRI, triangles, recursive);
585 Range side_ents3d_tris;
586 CHKERR moab.get_adjacencies(side_ents3d, 2, true, side_ents3d_tris,
587 moab::Interface::UNION);
588 triangles = intersect(triangles, side_ents3d_tris);
589 // nodes on interface but not on crack front (those should not be splitted)
590 Range nodes;
591 CHKERR moab.get_entities_by_type(children[2], MBVERTEX, nodes, false);
592 Range meshset_3d_ents, meshset_2d_ents;
593 CHKERR moab.get_entities_by_dimension(meshset, 3, meshset_3d_ents, true);
594 Range meshset_tets = meshset_3d_ents.subset_by_type(MBTET);
595 CHKERR moab.get_adjacencies(meshset_tets, 2, false, meshset_2d_ents,
596 moab::Interface::UNION);
597 side_ents3d = intersect(meshset_3d_ents, side_ents3d);
598 other_ents3d = intersect(meshset_3d_ents, other_ents3d);
599 triangles = intersect(meshset_2d_ents, triangles);
600
601 MOFEM_LOG_C("PRISM_INTERFACE", Sev::verbose, "Split sides triangles %u",
602 triangles.size());
603 MOFEM_LOG_C("PRISM_INTERFACE", Sev::verbose, "Split sides 3d entities %u",
604 side_ents3d.size());
605 MOFEM_LOG_C("PRISM_INTERFACE", Sev::verbose, "split sides nodes %u",
606 nodes.size());
607
608 struct PartentAndChild {
609 EntityHandle parent;
610 EntityHandle child;
611 };
612
613 typedef multi_index_container<
614 PartentAndChild,
615 indexed_by<
616
617 hashed_unique<
618 member<PartentAndChild, EntityHandle, &PartentAndChild::parent>>,
619
620 hashed_unique<
621 member<PartentAndChild, EntityHandle, &PartentAndChild::child>>
622
623 >>
624 ParentChildMI;
625
626 ParentChildMI parent_child;
627
628 typedef std::map<EntityHandle, /*node on "mother" side*/
629 EntityHandle /*node on "father" side*/
630 >
631 MapNodes;
632 MapNodes map_nodes, reverse_map_nodes;
633
634 // Map nodes on sides, set parent node and set bit ref level
635 {
636 struct CreateSideNodes {
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}
static Index< 'p', 3 > p
@ QUIET
Definition: definitions.h:208
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
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.
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
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1473
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:1432
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:1396
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ splitSides() [2/3]

MoFEMErrorCode MoFEM::PrismInterface::splitSides ( const EntityHandle  meshset,
const BitRefLevel bit,
const EntityHandle  sideset,
const bool  add_interface_entities,
const bool  recursive = false,
int  verb = QUIET 
)

Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in-between.

Parameters
meshsetvolume meshset containing 3D entities around the interface
bitbit ref level on which new entities will be stored
sidesetmeshset with surface
add_interface_entitiesif true add prism elements at interface
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
Parent meshset must have three child meshsets: two with tetrahedra from each side of the interface, third containing nodes which can be split and skin edges

Definition at line 546 of file PrismInterface.cpp.

550 {
551
553 CHKERR splitSides(meshset, bit, BitRefLevel(), BitRefLevel(), sideset,
554 add_interface_entities, recursive, verb);
556}
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...

◆ splitSides() [3/3]

MoFEMErrorCode MoFEM::PrismInterface::splitSides ( const EntityHandle  meshset,
const BitRefLevel bit,
const int  msId,
const CubitBCType  cubit_bc_type,
const bool  add_interface_entities,
const bool  recursive = false,
int  verb = QUIET 
)

Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in-between.

Parameters
meshsetvolume meshset containing 3D entities around the interface
bitbit ref level on which new entities will be stored
msIdmeshset ID of the surface
cubit_bc_typetype of meshset (NODESET, SIDESET or BLOCKSET and more)
add_interface_entitiesif true add prism elements at interface
recursiveif true parent meshset is searched recursively
verbverbosity level
Note
Parent meshset must have three child meshsets: two with tetrahedra from each side of the interface, third containing nodes which can be split and skin edges
Examples
continuity_check_on_contact_prism_side_ele.cpp, forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, simple_contact.cpp, and test_jacobian_of_simple_contact_element.cpp.

Definition at line 519 of file PrismInterface.cpp.

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

Member Data Documentation

◆ cOre

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

Definition at line 28 of file PrismInterface.hpp.


The documentation for this struct was generated from the following files: