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:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
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:406

◆ 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:476
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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:476
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
#define CHKERR
Inline error check.
Definition: definitions.h:595
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:406

◆ 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:500
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
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:476
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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:476
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:595
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:406

◆ 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 {
689  typedef RefEntity_multiIndex::index<
690  Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
692  operator()(const BitRefLevel &bit, const BitRefLevel &mask,
693  const RefEntity_multiIndex *refined_ents_ptr,
695  &ref_parent_ents_view) const {
697  auto &ref_ents =
698  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
699  // view by parent type (VERTEX)
700  auto miit = ref_ents.lower_bound(boost::make_tuple(MBVERTEX, MBVERTEX));
701  auto hi_miit =
702  ref_ents.upper_bound(boost::make_tuple(MBVERTEX, MBVERTEX));
703  for (; miit != hi_miit; miit++) {
704  const auto &ent_bit = (*miit)->getBitRefLevel();
705  if ((ent_bit & bit).any() && (ent_bit & mask) == ent_bit) {
706  auto p_ref_ent_view = ref_parent_ents_view.insert(*miit);
707  if (!p_ref_ent_view.second)
708  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
709  "non unique insertion");
710  }
711  }
713  }
714  };
715  if (inhered_from_bit_level.any() && inhered_from_bit_level_mask.any())
716  CHKERR CreateParentEntView()(inhered_from_bit_level,
717  inhered_from_bit_level_mask,
718  refined_ents_ptr, ref_parent_ents_view);
719 
720  // add new nodes on interface and create map
721  Range add_bit_nodes;
722  for (auto pnit = nodes.pair_begin(); pnit != nodes.pair_end(); ++pnit) {
723  auto lo = refined_ents_ptr->lower_bound(pnit->first);
724  auto hi = refined_ents_ptr->upper_bound(pnit->second);
725  if (std::distance(lo, hi) != (pnit->second - pnit->first + 1))
726  SETERRQ(
727  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
728  "Can not find some nodes in database that are split on interface");
729  Range nodes_in_range;
730  insertOrdered(nodes_in_range, RefEntExtractor(), lo, hi);
731  std::vector<double> coords_range(nodes_in_range.size() * 3);
732  CHKERR moab.get_coords(nodes_in_range, &*coords_range.begin());
733  int pos = 0;
734  for (; lo != hi; ++lo, pos += 3) {
735  const EntityHandle node = (*lo)->getRefEnt();
736  EntityHandle child_entity = 0;
737  auto child_it = ref_parent_ents_view.find(node);
738  if (child_it != ref_parent_ents_view.end())
739  child_entity = (*child_it)->getRefEnt();
740  if (child_entity == 0) {
741  CHKERR create_side_nodes(&coords_range[pos], node);
742  } else {
743  map_nodes[node] = child_entity;
744  add_bit_nodes.insert(child_entity);
745  }
746  }
747  }
748  add_bit_nodes.merge(nodes);
749  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit_nodes,
750  bit);
751  CHKERR create_side_nodes(bit, map_nodes);
752  }
753 
754  // crete meshset for new mesh bit level
755  EntityHandle meshset_for_bit_level;
756  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
757  meshset_for_bit_level);
758  // subtract those elements which will be refined, i.e. disconnected form other
759  // side elements, and connected to new prisms, if they area created
760  meshset_3d_ents = subtract(meshset_3d_ents, side_ents3d);
761  CHKERR moab.add_entities(meshset_for_bit_level, meshset_3d_ents);
762 
763  // create new 3d ents on "father" side
764  Range new_3d_ents;
765  for (Range::iterator eit3d = side_ents3d.begin(); eit3d != side_ents3d.end();
766  eit3d++) {
767  auto miit_ref_ent = refined_ents_ptr->find(*eit3d);
768  if (miit_ref_ent == refined_ents_ptr->end())
769  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
770  "Tetrahedron not in database");
771 
772  int num_nodes;
773  const EntityHandle *conn;
774  CHKERR moab.get_connectivity(*eit3d, conn, num_nodes, true);
775  EntityHandle new_conn[num_nodes];
776  int nb_new_conn = 0;
777  for (int ii = 0; ii < num_nodes; ii++) {
778  std::map<EntityHandle, EntityHandle>::iterator mit =
779  map_nodes.find(conn[ii]);
780  if (mit != map_nodes.end()) {
781  new_conn[ii] = mit->second;
782  nb_new_conn++;
783  if (verb >= VERY_NOISY) {
784  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
785  new_conn[ii]);
786  }
787  } else {
788  new_conn[ii] = conn[ii];
789  }
790  }
791  if (nb_new_conn == 0) {
792  // Add this tet to bit ref level
793  CHKERR moab.add_entities(meshset_for_bit_level, &*eit3d, 1);
794  continue;
795  }
796 
797  // here is created new or prism is on interface
798  EntityHandle existing_ent = 0;
799  /* check if tet element with new connectivity is in database*/
800  auto child_iit =
801  refined_ents_ptr->get<Ent_Ent_mi_tag>().lower_bound(*eit3d);
802  auto hi_child_iit =
803  refined_ents_ptr->get<Ent_Ent_mi_tag>().upper_bound(*eit3d);
804 
805  // Check if child entity has the same connectivity
806  for (; child_iit != hi_child_iit; child_iit++) {
807  const EntityHandle *conn_ref_tet;
808  CHKERR moab.get_connectivity(child_iit->get()->getRefEnt(), conn_ref_tet,
809  num_nodes, true);
810  int nn = 0;
811  for (; nn < num_nodes; nn++) {
812  if (conn_ref_tet[nn] != new_conn[nn]) {
813  break;
814  }
815  }
816  if (nn == num_nodes) {
817  if (existing_ent != 0)
818  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
819  "Should be only one child entity with the same connectivity");
820  existing_ent = child_iit->get()->getRefEnt();
821  }
822  }
823 
824  switch (moab.type_from_handle(*eit3d)) {
825  case MBTET: {
826 
827  RefEntity_multiIndex::iterator child_it;
828  EntityHandle tet;
829  if (existing_ent == 0) {
830  Range new_conn_tet;
831  CHKERR moab.get_adjacencies(new_conn, 4, 3, false, new_conn_tet);
832  if (new_conn_tet.empty()) {
833  CHKERR moab.create_element(MBTET, new_conn, 4, tet);
834  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &tet, 1,
835  &*eit3d);
836  } else {
837  // FIXME: That takes firs element form the list. Should throw error
838  // if is more than one or handle it properly.
839  RefEntity_multiIndex::index<Ent_mi_tag>::type::iterator new_rit =
840  refined_ents_ptr->get<Ent_mi_tag>().find(*new_conn_tet.begin());
841  if (new_rit == refined_ents_ptr->get<Ent_mi_tag>().end())
842  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
843  "Can't find entity in database");
844  tet = *new_conn_tet.begin();
845  }
846  } else {
847  tet = existing_ent;
848  }
849 
850  CHKERR moab.add_entities(meshset_for_bit_level, &tet, 1);
851  new_3d_ents.insert(tet);
852 
853  } break;
854  case MBPRISM: {
855  EntityHandle prism;
856  if (existing_ent == 0) {
857  Range new_conn_prism;
858  CHKERR moab.get_adjacencies(new_conn, 6, 3, false, new_conn_prism);
859  if (new_conn_prism.empty()) {
860  CHKERR moab.create_element(MBPRISM, new_conn, 6, prism);
861  CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), &prism, 1,
862  &*eit3d);
863  } else {
864  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
865  "It is prism with such connectivity, that case has to be "
866  "handled but this is not implemented");
867  }
868  } else {
869  prism = existing_ent;
870  }
871  CHKERR moab.add_entities(meshset_for_bit_level, &prism, 1);
872  new_3d_ents.insert(prism);
873  } break;
874  default:
875  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
876  "Not implemented element");
877  }
878  }
879 
880  auto get_adj_ents = [&](const bool create) {
881  Range adj;
882  // create new entities by adjacencies form new tets
883  CHKERR moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET), 2, create,
884  adj, moab::Interface::UNION);
885  CHKERR moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET), 1, create,
886  adj, moab::Interface::UNION);
887  return adj;
888  };
889 
890  Range new_ents_existing = get_adj_ents(false);
891  Range new_ents = subtract(get_adj_ents(true), new_ents_existing);
892 
893  // Tags for setting side
894  Tag th_interface_side;
895  const int def_side[] = {0};
896  CHKERR moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
897  th_interface_side, MB_TAG_CREAT | MB_TAG_SPARSE,
898  def_side);
899 
900  struct SetParent {
901  map<EntityHandle, EntityHandle> parentsToChange;
902  MoFEMErrorCode operator()(const EntityHandle ent, const EntityHandle parent,
903  const RefEntity_multiIndex *ref_ents_ptr,
904  MoFEM::Core &cOre) {
905  MoFEM::Interface &m_field = cOre;
907  auto it = ref_ents_ptr->find(ent);
908  if (it != ref_ents_ptr->end()) {
909  if (it->get()->getParentEnt() != parent && ent != parent)
910  parentsToChange[ent] = parent;
911  } else {
912  if (ent != parent)
913  CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
914  &ent, 1, &parent);
915  }
917  }
918  MoFEMErrorCode operator()(const RefEntity_multiIndex *ref_ents_ptr) {
920  for (auto mit = parentsToChange.begin(); mit != parentsToChange.end();
921  ++mit) {
922  auto it = ref_ents_ptr->find(mit->first);
923  bool success = const_cast<RefEntity_multiIndex *>(ref_ents_ptr)
924  ->modify(it, RefEntity_change_parent(mit->second));
925  if (!success)
926  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
927  "Impossible to set parent");
928  }
930  }
931  };
932  SetParent set_parent;
933 
934  // add new edges and triangles to mofem database
935  Range ents;
936  CHKERR moab.get_adjacencies(triangles, 1, false, ents,
937  moab::Interface::UNION);
938  ents.insert(triangles.begin(), triangles.end());
939  for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
940  int num_nodes;
941  const EntityHandle *conn;
942  CHKERR moab.get_connectivity(*eit, conn, num_nodes, true);
943  EntityHandle new_conn[num_nodes];
944  int nb_new_conn = 0;
945  for (int ii = 0; ii != num_nodes; ++ii) {
946  std::map<EntityHandle, EntityHandle>::iterator mit =
947  map_nodes.find(conn[ii]);
948  if (mit != map_nodes.end()) {
949  new_conn[ii] = mit->second;
950  nb_new_conn++;
951  if (verb >= VERY_NOISY) {
952  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
953  new_conn[ii]);
954  }
955  } else {
956  new_conn[ii] = conn[ii];
957  }
958  }
959  if (nb_new_conn == 0)
960  continue;
961  RefEntity_multiIndex::iterator miit_ref_ent = refined_ents_ptr->find(*eit);
962  if (miit_ref_ent == refined_ents_ptr->end())
963  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
964  "this entity (edge or tri) should be already in database");
965 
966  Range new_ent; // contains all entities (edges or triangles) added to mofem
967  // database
968  switch (moab.type_from_handle(*eit)) {
969  case MBTRI: {
970  // get entity based on its connectivity
971  CHKERR moab.get_adjacencies(new_conn, 3, 2, false, new_ent);
972  if (new_ent.size() != 1)
973  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
974  "this tri should be in moab database");
975  int new_side = 1;
976  CHKERR moab.tag_set_data(th_interface_side, &*new_ent.begin(), 1,
977  &new_side);
978  if (verb >= VERY_VERBOSE)
979  PetscPrintf(m_field.get_comm(), "new_ent %u\n", new_ent.size());
980  // add prism element
981  if (add_interface_entities) {
982  if (inhered_from_bit_level.any()) {
983  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
984  "not implemented for inhered_from_bit_level");
985  }
986  // set prism connectivity
987  EntityHandle prism_conn[6] = {conn[0], conn[1], conn[2],
988  new_conn[0], new_conn[1], new_conn[2]};
989  EntityHandle prism;
990  CHKERR moab.create_element(MBPRISM, prism_conn, 6, prism);
991  CHKERR moab.add_entities(meshset_for_bit_level, &prism, 1);
992  }
993  } break;
994  case MBEDGE: {
995  CHKERR moab.get_adjacencies(new_conn, 2, 1, false, new_ent);
996  if (new_ent.size() != 1) {
997  SETERRQ2(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
998  "this edge should be in moab database "
999  "new_ent.size() = %u nb_new_conn = %d",
1000  new_ent.size(), nb_new_conn);
1001  }
1002  } break;
1003  default:
1004  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1005  "Houston we have a problem !!!");
1006  }
1007  if (new_ent.size() != 1) {
1008  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1009  "new_ent.size() = %u, size always should be 1", new_ent.size());
1010  }
1011  CHKERR set_parent(new_ent[0], *eit, refined_ents_ptr, cOre);
1012  }
1013 
1014  // all other entities, some ents like triangles and faces on the side of tets
1015  auto all_others_adj_entities = [&](const bool create) {
1016  Range adj;
1017  for (auto d : {1, 2})
1018  CHKERR moab.get_adjacencies(side_ents3d.subset_by_type(MBTET), d, create,
1019  adj, moab::Interface::UNION);
1020  return adj;
1021  };
1022  Range side_adj_faces_and_edges = all_others_adj_entities(true);
1023 
1024  for (Range::iterator eit = side_adj_faces_and_edges.begin();
1025  eit != side_adj_faces_and_edges.end(); ++eit) {
1026  int num_nodes;
1027  const EntityHandle *conn;
1028  CHKERR moab.get_connectivity(*eit, conn, num_nodes, true);
1029  EntityHandle new_conn[num_nodes];
1030  int nb_new_conn = 0;
1031  for (int ii = 0; ii < num_nodes; ii++) {
1032  std::map<EntityHandle, EntityHandle>::iterator mit =
1033  map_nodes.find(conn[ii]);
1034  if (mit != map_nodes.end()) {
1035  new_conn[ii] = mit->second;
1036  nb_new_conn++;
1037  if (verb >= VERY_NOISY) {
1038  PetscPrintf(m_field.get_comm(), "nodes %u -> %d\n", conn[ii],
1039  new_conn[ii]);
1040  }
1041  } else {
1042  new_conn[ii] = conn[ii];
1043  }
1044  }
1045  if (nb_new_conn == 0)
1046  continue;
1047  RefEntity_multiIndex::iterator miit_ref_ent = refined_ents_ptr->find(*eit);
1048  if (miit_ref_ent == refined_ents_ptr->end()) {
1049  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1050  "entity should be in MoFem database, num_nodes = %d", num_nodes);
1051  }
1052  Range new_ent;
1053  switch (moab.type_from_handle(*eit)) {
1054  case MBEDGE:
1055  CHKERR moab.get_adjacencies(new_conn, 2, 1, false, new_ent);
1056  break;
1057  case MBTRI:
1058  CHKERR moab.get_adjacencies(new_conn, 3, 2, false, new_ent);
1059  break;
1060  default:
1061  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1062  "Houston we have a problem");
1063  }
1064  if (new_ent.size() != 1) {
1065  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1066  "database inconsistency, new_ent.size() = %u", new_ent.size());
1067  }
1068  CHKERR set_parent(new_ent[0], *eit, refined_ents_ptr, cOre);
1069  }
1070 
1071  // add new prisms which parents are part of other intefaces
1072  Range new_3d_prims = new_3d_ents.subset_by_type(MBPRISM);
1073  for (Range::iterator pit = new_3d_prims.begin(); pit != new_3d_prims.end();
1074  ++pit) {
1075  // get parent entity
1076  EntityHandle parent_prism;
1077  CHKERR moab.tag_get_data(cOre.get_th_RefParentHandle(), &*pit, 1,
1078  &parent_prism);
1079  if (moab.type_from_handle(parent_prism) != MBPRISM) {
1080  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1081  "this prism should have parent which is prism as well");
1082  }
1083  int num_nodes;
1084  // parent prism
1085  const EntityHandle *conn_parent;
1086  CHKERR moab.get_connectivity(parent_prism, conn_parent, num_nodes, true);
1087  Range face_side3_parent, face_side4_parent;
1088  CHKERR moab.get_adjacencies(conn_parent, 3, 2, false, face_side3_parent);
1089  CHKERR moab.get_adjacencies(&conn_parent[3], 3, 2, false,
1090  face_side4_parent);
1091  if (face_side3_parent.size() != 1) {
1092  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1093  "parent face3.size() = %u", face_side3_parent.size());
1094  }
1095  if (face_side4_parent.size() != 1) {
1096  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1097  "parent face4.size() = %u", face_side4_parent.size());
1098  }
1099  // new prism
1100  const EntityHandle *conn;
1101  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
1102  Range face_side3, face_side4;
1103  CHKERR moab.get_adjacencies(conn, 3, 2, false, face_side3);
1104  CHKERR moab.get_adjacencies(&conn[3], 3, 2, false, face_side4);
1105  if (face_side3.size() != 1) {
1106  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "face3 is missing");
1107  }
1108  if (face_side4.size() != 1) {
1109  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "face4 is missing");
1110  }
1111  std::vector<EntityHandle> face(2), parent_face(2);
1112  face[0] = *face_side3.begin();
1113  face[1] = *face_side4.begin();
1114  parent_face[0] = *face_side3_parent.begin();
1115  parent_face[1] = *face_side4_parent.begin();
1116  for (int ff = 0; ff != 2; ++ff) {
1117  if (parent_face[ff] == face[ff])
1118  continue;
1119  int interface_side;
1120  CHKERR moab.tag_get_data(th_interface_side, &parent_face[ff], 1,
1121  &interface_side);
1122  CHKERR moab.tag_set_data(th_interface_side, &face[ff], 1,
1123  &interface_side);
1124  EntityHandle parent_tri;
1125  CHKERR moab.tag_get_data(cOre.get_th_RefParentHandle(), &face[ff], 1,
1126  &parent_tri);
1127  if (parent_tri != parent_face[ff]) {
1128  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1129  "wrong parent %lu", parent_tri);
1130  }
1131  }
1132  }
1133 
1134  // finalise by adding new tets and prism ti bit level
1135  // FIXME: This is switch of, you can not change parent.
1136  // CHKERR set_parent(refined_ents_ptr);
1137 
1138  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
1139  meshset_for_bit_level, 3, bit);
1140  CHKERR moab.delete_entities(&meshset_for_bit_level, 1);
1141  CHKERR moab.clear_meshset(&children[0], 3);
1143 }
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:476
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:536
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:595
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:406
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: