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

Make interface on given faces and create flat prism in that space. More...

#include <src/interfaces/PrismInterface.hpp>

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

Public Member Functions

MoFEMErrorCode query_interface (const MOFEMuuid &uuid, UnknownInterface **iface) const
 
 PrismInterface (const MoFEM::Core &core)
 
 ~PrismInterface ()
 destructor More...
 
MoFEMErrorCode getSides (const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
 create two children meshsets in the meshset containing tetrahedral on two sides of faces More...
 
MoFEMErrorCode getSides (const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
 create two children meshsets in the meshset containing tetrahedral on two sides of faces More...
 
MoFEMErrorCode findIfTringleHasThreeNodesOnInternalSurfaceSkin (const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, Range &faces_with_three_nodes_on_front, int verb=QUIET)
 Find if triangle has 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 tetrahedral in children sets and add prism elements More...
 
MoFEMErrorCode splitSides (const EntityHandle meshset, const BitRefLevel &bit, const EntityHandle, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
 split nodes and other entities of tetrahedral in children sets and add prism elements 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 tetrahedrons in children sets and add prism elements More...
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (const MOFEMuuid &uuid, bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE , bool VERIFY = false>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, IFACE *&iface) const
 Get interface by uuid and return reference to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()
 
virtual MoFEMErrorCode getLibVersion (Version &version) const
 Get library version. More...
 
virtual const MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version) const
 Get database major version. More...
 
virtual MoFEMErrorCode getInterfaceVersion (Version &version) const
 Get database major version. More...
 
template<>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, UnknownInterface *&iface) const
 

Public Attributes

MoFEM::CorecOre
 

Additional Inherited Members

- Protected Member Functions inherited from MoFEM::UnknownInterface
boost::typeindex::type_index getClassIdx (const MOFEMuuid &uid) const
 Get type name for interface Id. More...
 
MOFEMuuid getUId (const boost::typeindex::type_index &class_idx) const
 Get interface Id for class name. More...
 

Detailed Description

Make interface on given faces and create flat prism in that space.

Todo:
FIXME Names of methods do not follow naming convention and are difficult to work with.
Examples
forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, and split_sideset.cpp.

Definition at line 37 of file PrismInterface.hpp.

Constructor & Destructor Documentation

◆ PrismInterface()

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

Definition at line 36 of file PrismInterface.cpp.

37  : cOre(const_cast<Core &>(core)) {}

◆ ~PrismInterface()

MoFEM::PrismInterface::~PrismInterface ( )

destructor

Definition at line 46 of file PrismInterface.hpp.

46 {}

Member Function Documentation

◆ findIfTringleHasThreeNodesOnInternalSurfaceSkin()

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

Find if triangle has three nodes on internal surface skin.

Internal surface skin is a set of edges in the body on the boundary of the surface. This set of edges is called a surface front. If the surface the face has three nodes on the surface front, none of the face nodes is split and should be removed from the surface. Otherwise, such a triangle cannot be split.

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

Definition at line 398 of file PrismInterface.cpp.

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

◆ getSides() [1/2]

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

create two children meshsets in the meshset containing tetrahedral on two sides of faces

Parameters
msIdId of meshset
cubit_bc_typetype of meshset (NODESET, SIDESET or BLOCKSET and more)
mesh_bit_leveladd interface on bit level is bit_level = BitRefLevel.set() then add interface on all bit levels
recursiveif true parent meshset is searched recursively
Examples
forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, and split_sideset.cpp.

Definition at line 39 of file PrismInterface.cpp.

42  {
43 
44  Interface &m_field = cOre;
45  MeshsetsManager *meshsets_manager_ptr;
47  CHKERR m_field.getInterface(meshsets_manager_ptr);
48  CubitMeshSet_multiIndex::index<
49  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator miit =
50  meshsets_manager_ptr->getMeshsetsMultindex()
51  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
52  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
53  if (miit != meshsets_manager_ptr->getMeshsetsMultindex()
54  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
55  .end()) {
56  CHKERR getSides(miit->meshset, mesh_bit_level, recursive, verb);
57  } else {
58  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "msId is not there");
59  }
61 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
create two children meshsets in the meshset containing tetrahedral on two sides of faces
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ getSides() [2/2]

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

create two children meshsets in the meshset containing tetrahedral on two sides of faces

Get tets adj to faces. Take skin form tets and get edges from that skin. Take skin form triangles (the face). Subtrac skin faces edges form skin edges in order to get edges on the boundary of the face which is in the volume of the body, but is not on the boundary. Each child set has a child containing nodes which can be split and skin edges. After that simply iterate under all tets on one side which are adjacent to the face are found. Side tets are stored in to children meshsets of the SIDESET meshset.

Definition at line 63 of file PrismInterface.cpp.

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

◆ query_interface()

MoFEMErrorCode MoFEM::PrismInterface::query_interface ( const MOFEMuuid uuid,
UnknownInterface **  iface 
) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 24 of file PrismInterface.cpp.

25  {
27  *iface = NULL;
28  if (uuid == IDD_MOFEMPrismInterface) {
29  *iface = const_cast<PrismInterface *>(this);
31  }
32  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
34 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static const MOFEMuuid IDD_MOFEMPrismInterface

◆ splitSides() [1/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 tetrahedral in children sets and add prism elements

The all new entities (prisms, tets) are added to refinement level given by bit

Parameters
meshsetmeshset to get entities from
BitRefLevelnew level where refinement would be stored
msIdmeshset ID imported from cubit
cubit_bc_typetype of meshset (NODESET, SIDESET or BLOCKSET and more)
add_interface_entitiesmeshset which contain the interface
recursiveif true parent meshset is searched recursively

Each inteface face has two tags, const int def_side[] = {0}; rval = moab.tag_get_handle("INTERFACE_SIDE",1,MB_TYPE_INTEGER, th_interface_side,MB_TAG_CREAT|MB_TAG_SPARSE,def_side); CHKERRQ_MOAB(rval);

            const EntityHandle def_node[] = {0};
            rval = moab.tag_get_handle("SIDE_INTFACE_ELEMENT",1,MB_TYPE_HANDLE,
              th_side_elem,MB_TAG_CREAT|MB_TAG_SPARSE,def_node); CHKERRQ_MOAB(rval);

First tag inform about inteface side, second tag inform about side adjacent inteface element.

Examples
forces_and_sources_testing_flat_prism_element.cpp, mesh_insert_interface_atom.cpp, and split_sideset.cpp.

Definition at line 543 of file PrismInterface.cpp.

548  {
549 
550  Interface &m_field = cOre;
551  MeshsetsManager *meshsets_manager_ptr;
553  CHKERR m_field.getInterface(meshsets_manager_ptr);
554  CubitMeshSet_multiIndex::index<
555  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator miit =
556  meshsets_manager_ptr->getMeshsetsMultindex()
557  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
558  .find(boost::make_tuple(msId, cubit_bc_type.to_ulong()));
559  if (miit != meshsets_manager_ptr->getMeshsetsMultindex()
560  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
561  .end()) {
562  CHKERR splitSides(meshset, bit, miit->meshset, add_interface_entities,
563  recursive, verb);
564  } else {
565  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "msId is not there");
566  }
568 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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 tetrahedral in children sets and add prism elements
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ 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 tetrahedral in children sets and add prism elements

The all new entities (prisms, tets) are added to refinement level given by bit

Definition at line 570 of file PrismInterface.cpp.

574  {
575 
577  CHKERR splitSides(meshset, bit, BitRefLevel(), BitRefLevel(), sideset,
578  add_interface_entities, recursive, verb);
580 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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 tetrahedral in children sets and add prism elements
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ splitSides() [3/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 tetrahedrons in children sets and add prism elements

The all new entities (prisms, tets) are added to refinement level given by bit

Parameters
meshset
Refinementbit level of new mesh
inhered_from_bit_levelinhered nodes and other entities form this bit level.
add_interface_entitiesadd prism elements at interface
recursivedo meshesets in the meshset
Note
inhered_from_bit_level is need to be specified to some meshsets with interfaces. Some nodes on some refinement levels dividing edges but not splitting faces. Inheriting those nodes will not split faces.

Definition at line 581 of file PrismInterface.cpp.

585  {
586 
587  Interface &m_field = cOre;
588  moab::Interface &moab = m_field.get_moab();
589  const RefEntity_multiIndex *refined_ents_ptr;
591 
592  CHKERR m_field.get_ref_ents(&refined_ents_ptr);
593 
594  std::vector<EntityHandle> children;
595  // get children meshsets
596  CHKERR moab.get_child_meshsets(sideset, children);
597  if (children.size() != 3)
598  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
599  "should be 3 child meshsets, each of them contains tets on two "
600  "sides of interface");
601 
602  // 3d ents on "father" side
603  Range side_ents3d;
604  CHKERR moab.get_entities_by_handle(children[0], side_ents3d, false);
605  // 3d ents on "mather" side
606  Range other_ents3d;
607  CHKERR moab.get_entities_by_handle(children[1], other_ents3d, false);
608  // faces of interface
609  Range triangles;
610  CHKERR moab.get_entities_by_type(sideset, MBTRI, triangles, recursive);
611  Range side_ents3d_tris;
612  CHKERR moab.get_adjacencies(side_ents3d, 2, true, side_ents3d_tris,
613  moab::Interface::UNION);
614  triangles = intersect(triangles, side_ents3d_tris);
615  // nodes on interface but not on crack front (those should not be splitted)
616  Range nodes;
617  CHKERR moab.get_entities_by_type(children[2], MBVERTEX, nodes, false);
618  Range meshset_3d_ents, meshset_2d_ents;
619  CHKERR moab.get_entities_by_dimension(meshset, 3, meshset_3d_ents, true);
620  Range meshset_tets = meshset_3d_ents.subset_by_type(MBTET);
621  CHKERR moab.get_adjacencies(meshset_tets, 2, false, meshset_2d_ents,
622  moab::Interface::UNION);
623  side_ents3d = intersect(meshset_3d_ents, side_ents3d);
624  other_ents3d = intersect(meshset_3d_ents, other_ents3d);
625  triangles = intersect(meshset_2d_ents, triangles);
626  if (verb >= VERY_VERBOSE) {
627  PetscPrintf(m_field.get_comm(), "split sides triangles %u\n",
628  triangles.size());
629  PetscPrintf(m_field.get_comm(), "split sides side_ents3d %u\n",
630  side_ents3d.size());
631  PetscPrintf(m_field.get_comm(), "split sides nodes %u\n", nodes.size());
632  }
633 
634  typedef std::map<EntityHandle, /*node on "mather" side*/
635  EntityHandle /*node on "father" side*/
636  >
637  MapNodes;
638  MapNodes map_nodes;
639 
640  // Map nodes on sides, set parent node and set bit ref level
641  {
642  struct CreateSideNodes {
643  MoFEM::Core &cOre;
644  MoFEM::Interface &m_field;
645  std::vector<EntityHandle> splitNodes;
646  std::vector<double> splitCoords[3];
647 
649 
650  CreateSideNodes(MoFEM::Core &core, int split_size = 0)
651  : cOre(core), m_field(core) {
652  splitNodes.reserve(split_size);
653  for (auto dd : {0, 1, 2})
654  splitCoords[dd].reserve(split_size);
655  }
656  MoFEMErrorCode operator()(const double coords[], const EntityHandle n) {
658  splitNodes.emplace_back(n);
659  for (auto dd : {0, 1, 2})
660  splitCoords[dd].emplace_back(coords[dd]);
662  }
663 
664  MoFEMErrorCode operator()(const BitRefLevel &bit, MapNodes &map_nodes) {
665  ReadUtilIface *iface;
667  int num_nodes = splitNodes.size();
668  std::vector<double *> arrays_coord;
669  EntityHandle startv;
670  CHKERR m_field.get_moab().query_interface(iface);
671  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays_coord);
672  Range verts(startv, startv + num_nodes - 1);
673  for (int dd = 0; dd != 3; ++dd)
674  std::copy(splitCoords[dd].begin(), splitCoords[dd].end(),
675  arrays_coord[dd]);
676  for (int nn = 0; nn != num_nodes; ++nn)
677  map_nodes[splitNodes[nn]] = verts[nn];
678  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
679  verts, &*splitNodes.begin());
680  CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
681  verts, bit, QUIET);
683  }
684  };
685  CreateSideNodes create_side_nodes(cOre, nodes.size());
686 
688  struct CreateParentEntView {
690  operator()(const BitRefLevel &bit, const BitRefLevel &mask,
691  const RefEntity_multiIndex *refined_ents_ptr,
693  &ref_parent_ents_view) const {
695  auto &ref_ents =
696  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
697  // view by parent type (VERTEX)
698  auto miit = ref_ents.lower_bound(boost::make_tuple(MBVERTEX, MBVERTEX));
699  auto hi_miit =
700  ref_ents.upper_bound(boost::make_tuple(MBVERTEX, MBVERTEX));
701  for (; miit != hi_miit; miit++) {
702  const auto &ent_bit = (*miit)->getBitRefLevel();
703  if ((ent_bit & bit).any() && (ent_bit & mask) == ent_bit) {
704  auto p_ref_ent_view = ref_parent_ents_view.insert(*miit);
705  if (!p_ref_ent_view.second)
706  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
707  "non unique insertion");
708  }
709  }
711  }
712  };
713  if (inhered_from_bit_level.any() && inhered_from_bit_level_mask.any())
714  CHKERR CreateParentEntView()(inhered_from_bit_level,
715  inhered_from_bit_level_mask,
716  refined_ents_ptr, ref_parent_ents_view);
717 
718  // add new nodes on interface and create map
719  Range add_bit_nodes;
720  for (auto pnit = nodes.pair_begin(); pnit != nodes.pair_end(); ++pnit) {
721  auto lo = refined_ents_ptr->lower_bound(pnit->first);
722  auto hi = refined_ents_ptr->upper_bound(pnit->second);
723  if (std::distance(lo, hi) != (pnit->second - pnit->first + 1))
724  SETERRQ(
725  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
726  "Can not find some nodes in database that are split on interface");
727  Range nodes_in_range;
728  insertOrdered(nodes_in_range, RefEntExtractor(), lo, hi);
729  std::vector<double> coords_range(nodes_in_range.size() * 3);
730  CHKERR moab.get_coords(nodes_in_range, &*coords_range.begin());
731  int pos = 0;
732  for (; lo != hi; ++lo, pos += 3) {
733  const EntityHandle node = (*lo)->getRefEnt();
734  EntityHandle child_entity = 0;
735  auto child_it = ref_parent_ents_view.find(node);
736  if (child_it != ref_parent_ents_view.end())
737  child_entity = (*child_it)->getRefEnt();
738  if (child_entity == 0) {
739  CHKERR create_side_nodes(&coords_range[pos], node);
740  } else {
741  map_nodes[node] = child_entity;
742  add_bit_nodes.insert(child_entity);
743  }
744  }
745  }
746  add_bit_nodes.merge(nodes);
747  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit_nodes,
748  bit);
749  CHKERR create_side_nodes(bit, map_nodes);
750  }
751 
752  // crete meshset for new mesh bit level
753  EntityHandle meshset_for_bit_level;
754  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
755  meshset_for_bit_level);
756  // subtract those elements which will be refined, i.e. disconnected form other
757  // side elements, and connected to new prisms, if they area created
758  meshset_3d_ents = subtract(meshset_3d_ents, side_ents3d);
759  CHKERR moab.add_entities(meshset_for_bit_level, meshset_3d_ents);
760 
761  // create new 3d ents on "father" side
762  Range new_3d_ents;
763  for (Range::iterator eit3d = side_ents3d.begin(); eit3d != side_ents3d.end();
764  eit3d++) {
765  auto miit_ref_ent = refined_ents_ptr->find(*eit3d);
766  if (miit_ref_ent == refined_ents_ptr->end())
767  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
768  "Tetrahedron not in database");
769 
770  int num_nodes;
771  const EntityHandle *conn;
772  CHKERR moab.get_connectivity(*eit3d, conn, num_nodes, true);
773  EntityHandle new_conn[num_nodes];
774  int nb_new_conn = 0;
775  for (int ii = 0; ii < num_nodes; ii++) {
776  std::map<EntityHandle, EntityHandle>::iterator mit =
777  map_nodes.find(conn[ii]);
778  if (mit != map_nodes.end()) {
779  new_conn[ii] = mit->second;
780  nb_new_conn++;
781  if (verb >= VERY_NOISY) {
782  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
783  new_conn[ii]);
784  }
785  } else {
786  new_conn[ii] = conn[ii];
787  }
788  }
789  if (nb_new_conn == 0) {
790  // Add this tet to bit ref level
791  CHKERR moab.add_entities(meshset_for_bit_level, &*eit3d, 1);
792  continue;
793  }
794 
795  // here is created new or prism is on interface
796  EntityHandle existing_ent = 0;
797  /* check if tet element with new connectivity is in database*/
798  auto child_iit =
799  refined_ents_ptr->get<Ent_Ent_mi_tag>().lower_bound(*eit3d);
800  auto hi_child_iit =
801  refined_ents_ptr->get<Ent_Ent_mi_tag>().upper_bound(*eit3d);
802 
803  // Check if child entity has the same connectivity
804  for (; child_iit != hi_child_iit; child_iit++) {
805  const EntityHandle *conn_ref_tet;
806  CHKERR moab.get_connectivity(child_iit->get()->getRefEnt(), conn_ref_tet,
807  num_nodes, true);
808  int nn = 0;
809  for (; nn < num_nodes; nn++) {
810  if (conn_ref_tet[nn] != new_conn[nn]) {
811  break;
812  }
813  }
814  if (nn == num_nodes) {
815  if (existing_ent != 0)
816  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
817  "Should be only one child entity with the same connectivity");
818  existing_ent = child_iit->get()->getRefEnt();
819  }
820  }
821 
822  switch (moab.type_from_handle(*eit3d)) {
823  case MBTET: {
824 
825  RefEntity_multiIndex::iterator child_it;
826  EntityHandle tet;
827  if (existing_ent == 0) {
828  Range new_conn_tet;
829  CHKERR moab.get_adjacencies(new_conn, 4, 3, false, new_conn_tet);
830  if (new_conn_tet.empty()) {
831  CHKERR moab.create_element(MBTET, new_conn, 4, tet);
832  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &tet, 1,
833  &*eit3d);
834  } else {
835  // FIXME: That takes firs element form the list. Should throw error
836  // if is more than one or handle it properly.
837  RefEntity_multiIndex::index<Ent_mi_tag>::type::iterator new_rit =
838  refined_ents_ptr->get<Ent_mi_tag>().find(*new_conn_tet.begin());
839  if (new_rit == refined_ents_ptr->get<Ent_mi_tag>().end())
840  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
841  "Can't find entity in database");
842  tet = *new_conn_tet.begin();
843  }
844  } else {
845  tet = existing_ent;
846  }
847 
848  CHKERR moab.add_entities(meshset_for_bit_level, &tet, 1);
849  new_3d_ents.insert(tet);
850 
851  } break;
852  case MBPRISM: {
853  EntityHandle prism;
854  if (existing_ent == 0) {
855  Range new_conn_prism;
856  CHKERR moab.get_adjacencies(new_conn, 6, 3, false, new_conn_prism);
857  if (new_conn_prism.empty()) {
858  CHKERR moab.create_element(MBPRISM, new_conn, 6, prism);
859  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &prism, 1,
860  &*eit3d);
861  } else {
862  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
863  "It is prism with such connectivity, that case has to be "
864  "handled but this is not implemented");
865  }
866  } else {
867  prism = existing_ent;
868  }
869  CHKERR moab.add_entities(meshset_for_bit_level, &prism, 1);
870  new_3d_ents.insert(prism);
871  } break;
872  default:
873  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
874  "Not implemented element");
875  }
876  }
877 
878  auto get_adj_ents = [&](const bool create) {
879  Range adj;
880  // create new entities by adjacencies form new tets
881  CHKERR moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET), 2, create,
882  adj, moab::Interface::UNION);
883  CHKERR moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET), 1, create,
884  adj, moab::Interface::UNION);
885  return adj;
886  };
887 
888  Range new_ents_existing = get_adj_ents(false);
889  Range new_ents = subtract(get_adj_ents(true), new_ents_existing);
890 
891  // Tags for setting side
892  Tag th_interface_side;
893  const int def_side[] = {0};
894  CHKERR moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
895  th_interface_side, MB_TAG_CREAT | MB_TAG_SPARSE,
896  def_side);
897 
898  struct SetParent {
899  map<EntityHandle, EntityHandle> parentsToChange;
900  MoFEMErrorCode operator()(const EntityHandle ent, const EntityHandle parent,
901  const RefEntity_multiIndex *ref_ents_ptr,
902  MoFEM::Core &cOre) {
903  MoFEM::Interface &m_field = cOre;
905  auto it = ref_ents_ptr->find(ent);
906  if (it != ref_ents_ptr->end()) {
907  if (it->get()->getParentEnt() != parent && ent != parent)
908  parentsToChange[ent] = parent;
909  } else {
910  if (ent != parent)
911  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
912  &ent, 1, &parent);
913  }
915  }
916  MoFEMErrorCode operator()(const RefEntity_multiIndex *ref_ents_ptr) {
918  for (auto mit = parentsToChange.begin(); mit != parentsToChange.end();
919  ++mit) {
920  auto it = ref_ents_ptr->find(mit->first);
921  bool success = const_cast<RefEntity_multiIndex *>(ref_ents_ptr)
922  ->modify(it, RefEntity_change_parent(mit->second));
923  if (!success)
924  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
925  "Impossible to set parent");
926  }
928  }
929  };
930  SetParent set_parent;
931 
932  // add new edges and triangles to mofem database
933  Range ents;
934  CHKERR moab.get_adjacencies(triangles, 1, false, ents,
935  moab::Interface::UNION);
936  ents.insert(triangles.begin(), triangles.end());
937  for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
938  int num_nodes;
939  const EntityHandle *conn;
940  CHKERR moab.get_connectivity(*eit, conn, num_nodes, true);
941  EntityHandle new_conn[num_nodes];
942  int nb_new_conn = 0;
943  for (int ii = 0; ii != num_nodes; ++ii) {
944  std::map<EntityHandle, EntityHandle>::iterator mit =
945  map_nodes.find(conn[ii]);
946  if (mit != map_nodes.end()) {
947  new_conn[ii] = mit->second;
948  nb_new_conn++;
949  if (verb >= VERY_NOISY) {
950  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
951  new_conn[ii]);
952  }
953  } else {
954  new_conn[ii] = conn[ii];
955  }
956  }
957  if (nb_new_conn == 0)
958  continue;
959  RefEntity_multiIndex::iterator miit_ref_ent = refined_ents_ptr->find(*eit);
960  if (miit_ref_ent == refined_ents_ptr->end())
961  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
962  "this entity (edge or tri) should be already in database");
963 
964  Range new_ent; // contains all entities (edges or triangles) added to mofem
965  // database
966  switch (moab.type_from_handle(*eit)) {
967  case MBTRI: {
968  // get entity based on its connectivity
969  CHKERR moab.get_adjacencies(new_conn, 3, 2, false, new_ent);
970  if (new_ent.size() != 1)
971  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
972  "this tri should be in moab database");
973  int new_side = 1;
974  CHKERR moab.tag_set_data(th_interface_side, &*new_ent.begin(), 1,
975  &new_side);
976  if (verb >= VERY_VERBOSE)
977  PetscPrintf(m_field.get_comm(), "new_ent %u\n", new_ent.size());
978  // add prism element
979  if (add_interface_entities) {
980  if (inhered_from_bit_level.any()) {
981  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
982  "not implemented for inhered_from_bit_level");
983  }
984  // set prism connectivity
985  EntityHandle prism_conn[6] = {conn[0], conn[1], conn[2],
986  new_conn[0], new_conn[1], new_conn[2]};
987  EntityHandle prism;
988  CHKERR moab.create_element(MBPRISM, prism_conn, 6, prism);
989  CHKERR moab.add_entities(meshset_for_bit_level, &prism, 1);
990  }
991  } break;
992  case MBEDGE: {
993  CHKERR moab.get_adjacencies(new_conn, 2, 1, false, new_ent);
994  if (new_ent.size() != 1) {
995  SETERRQ2(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
996  "this edge should be in moab database "
997  "new_ent.size() = %u nb_new_conn = %d",
998  new_ent.size(), nb_new_conn);
999  }
1000  } break;
1001  default:
1002  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1003  "Houston we have a problem !!!");
1004  }
1005  if (new_ent.size() != 1) {
1006  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1007  "new_ent.size() = %u, size always should be 1", new_ent.size());
1008  }
1009  CHKERR set_parent(new_ent[0], *eit, refined_ents_ptr, cOre);
1010  }
1011 
1012  // all other entities, some ents like triangles and faces on the side of tets
1013  auto all_others_adj_entities = [&](const bool create) {
1014  Range adj;
1015  for (auto d : {1, 2})
1016  CHKERR moab.get_adjacencies(side_ents3d.subset_by_type(MBTET), d, create,
1017  adj, moab::Interface::UNION);
1018  return adj;
1019  };
1020  Range side_adj_faces_and_edges = all_others_adj_entities(true);
1021 
1022  for (Range::iterator eit = side_adj_faces_and_edges.begin();
1023  eit != side_adj_faces_and_edges.end(); ++eit) {
1024  int num_nodes;
1025  const EntityHandle *conn;
1026  CHKERR moab.get_connectivity(*eit, conn, num_nodes, true);
1027  EntityHandle new_conn[num_nodes];
1028  int nb_new_conn = 0;
1029  for (int ii = 0; ii < num_nodes; ii++) {
1030  std::map<EntityHandle, EntityHandle>::iterator mit =
1031  map_nodes.find(conn[ii]);
1032  if (mit != map_nodes.end()) {
1033  new_conn[ii] = mit->second;
1034  nb_new_conn++;
1035  if (verb >= VERY_NOISY) {
1036  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
1037  new_conn[ii]);
1038  }
1039  } else {
1040  new_conn[ii] = conn[ii];
1041  }
1042  }
1043  if (nb_new_conn == 0)
1044  continue;
1045  RefEntity_multiIndex::iterator miit_ref_ent = refined_ents_ptr->find(*eit);
1046  if (miit_ref_ent == refined_ents_ptr->end()) {
1047  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1048  "entity should be in MoFem database, num_nodes = %d", num_nodes);
1049  }
1050  Range new_ent;
1051  switch (moab.type_from_handle(*eit)) {
1052  case MBEDGE:
1053  CHKERR moab.get_adjacencies(new_conn, 2, 1, false, new_ent);
1054  break;
1055  case MBTRI:
1056  CHKERR moab.get_adjacencies(new_conn, 3, 2, false, new_ent);
1057  break;
1058  default:
1059  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1060  "Houston we have a problem");
1061  }
1062  if (new_ent.size() != 1) {
1063  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1064  "database inconsistency, new_ent.size() = %u", new_ent.size());
1065  }
1066  CHKERR set_parent(new_ent[0], *eit, refined_ents_ptr, cOre);
1067  }
1068 
1069  // add new prisms which parents are part of other intefaces
1070  Range new_3d_prims = new_3d_ents.subset_by_type(MBPRISM);
1071  for (Range::iterator pit = new_3d_prims.begin(); pit != new_3d_prims.end();
1072  ++pit) {
1073  // get parent entity
1074  EntityHandle parent_prism;
1075  CHKERR moab.tag_get_data(cOre.get_th_RefParentHandle(), &*pit, 1,
1076  &parent_prism);
1077  if (moab.type_from_handle(parent_prism) != MBPRISM) {
1078  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1079  "this prism should have parent which is prism as well");
1080  }
1081  int num_nodes;
1082  // parent prism
1083  const EntityHandle *conn_parent;
1084  CHKERR moab.get_connectivity(parent_prism, conn_parent, num_nodes, true);
1085  Range face_side3_parent, face_side4_parent;
1086  CHKERR moab.get_adjacencies(conn_parent, 3, 2, false, face_side3_parent);
1087  CHKERR moab.get_adjacencies(&conn_parent[3], 3, 2, false,
1088  face_side4_parent);
1089  if (face_side3_parent.size() != 1) {
1090  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1091  "parent face3.size() = %u", face_side3_parent.size());
1092  }
1093  if (face_side4_parent.size() != 1) {
1094  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1095  "parent face4.size() = %u", face_side4_parent.size());
1096  }
1097  // new prism
1098  const EntityHandle *conn;
1099  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
1100  Range face_side3, face_side4;
1101  CHKERR moab.get_adjacencies(conn, 3, 2, false, face_side3);
1102  CHKERR moab.get_adjacencies(&conn[3], 3, 2, false, face_side4);
1103  if (face_side3.size() != 1) {
1104  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "face3 is missing");
1105  }
1106  if (face_side4.size() != 1) {
1107  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "face4 is missing");
1108  }
1109  std::vector<EntityHandle> face(2), parent_face(2);
1110  face[0] = *face_side3.begin();
1111  face[1] = *face_side4.begin();
1112  parent_face[0] = *face_side3_parent.begin();
1113  parent_face[1] = *face_side4_parent.begin();
1114  for (int ff = 0; ff != 2; ++ff) {
1115  if (parent_face[ff] == face[ff])
1116  continue;
1117  int interface_side;
1118  CHKERR moab.tag_get_data(th_interface_side, &parent_face[ff], 1,
1119  &interface_side);
1120  CHKERR moab.tag_set_data(th_interface_side, &face[ff], 1,
1121  &interface_side);
1122  EntityHandle parent_tri;
1123  CHKERR moab.tag_get_data(cOre.get_th_RefParentHandle(), &face[ff], 1,
1124  &parent_tri);
1125  if (parent_tri != parent_face[ff]) {
1126  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1127  "wrong parent %lu", parent_tri);
1128  }
1129  }
1130  }
1131 
1132  // finalise by adding new tets and prism ti bit level
1133  // FIXME: This is switch of, you can not change parent.
1134  // CHKERR set_parent(refined_ents_ptr);
1135 
1136  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
1137  meshset_for_bit_level, 3, bit);
1138  CHKERR moab.delete_entities(&meshset_for_bit_level, 1);
1139  CHKERR moab.clear_meshset(&children[0], 3);
1141 }
Tag get_th_RefParentHandle() const
Definition: Core.hpp:150
virtual moab::Interface & get_moab()=0
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
Core (interface) class.
Definition: Core.hpp:50
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:590
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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::getRefEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_hashed_parent_entity
multi-index view of RefEntity by parent entity
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, 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::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0

Member Data Documentation

◆ cOre

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

Definition at line 42 of file PrismInterface.hpp.


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