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