v0.5.86
Public Member Functions | Public Attributes | List of all members
MoFEM::PrismInterface Struct Reference
MoFEM » Prism interface

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

PetscErrorCode queryInterface (const MOFEMuuid &uuid, UnknownInterface **iface)
 
 PrismInterface (const MoFEM::Core &core)
 
 ~PrismInterface ()
 destructor More...
 
virtual PetscErrorCode getSides (const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=-1)
 create two children meshsets in the meshset containing tetrahedral on two sides of faces More...
 
virtual PetscErrorCode getSides (const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, int verb=-1)
 create two children meshsets in the meshset containing tetrahedral on two sides of faces More...
 
virtual PetscErrorCode splitSides (const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_iterfece_entities, const bool recursive=false, int verb=-1)
 split nodes and other entities of tetrahedral in children sets and add prism elements More...
 
virtual PetscErrorCode splitSides (const EntityHandle meshset, const BitRefLevel &bit, const EntityHandle sideset, const bool add_iterfece_entities, const bool recursive=false, int verb=-1)
 split nodes and other entities of tetrahedral in children sets and add prism elements More...
 
virtual PetscErrorCode splitSides (const EntityHandle meshset, const BitRefLevel &bit, const BitRefLevel &inheret_from_bit_level, const BitRefLevel &inheret_from_bit_level_mask, const EntityHandle sideset, const bool add_iterfece_entities, const bool recursive=false, int verb=-1)
 split nodes and other entities of tetrahedrons in children sets and add prism elements More...
 
- Public Member Functions inherited from MoFEM::UnknownInterface
virtual ~UnknownInterface ()
 
virtual PetscErrorCode getLibVersion (Version &version) const
 Get library version. More...
 
virtual const PetscErrorCode getFileVersion (moab::Interface &moab, Version &version) const
 Get database major version. More...
 
virtual PetscErrorCode getInterfaceVersion (Version &version) const
 Get database major version. More...
 

Public Attributes

MoFEM::Core & cOre
 

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.

Definition at line 36 of file PrismInterface.hpp.

Constructor & Destructor Documentation

◆ PrismInterface()

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

Definition at line 70 of file PrismInterface.cpp.

70  :
71 cOre(const_cast<MoFEM::Core&>(core)) {
72 }

◆ ~PrismInterface()

MoFEM::PrismInterface::~PrismInterface ( )

destructor

Definition at line 44 of file PrismInterface.hpp.

44 {}

Member Function Documentation

◆ getSides() [1/2]

PetscErrorCode MoFEM::PrismInterface::getSides ( const int  msId,
const CubitBCType  cubit_bc_type,
const BitRefLevel  mesh_bit_level,
const bool  recursive,
int  verb = -1 
)
virtual

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

Definition at line 74 of file PrismInterface.cpp.

76  {
77 
78  MoFEM::Interface &m_field = cOre;
79  MeshsetsManager *meshsets_manager_ptr;
80  PetscFunctionBegin;
81  ierr = m_field.query_interface(meshsets_manager_ptr); CHKERRQ(ierr);
82  CubitMeshSet_multiIndex::index<Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator
83  miit = meshsets_manager_ptr->getMeshsetsMultindex().get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
84  .find(boost::make_tuple(msId,cubit_bc_type.to_ulong()));
85  if(miit!=meshsets_manager_ptr->getMeshsetsMultindex().get<Composite_Cubit_msId_And_MeshSetType_mi_tag>().end()) {
86  ierr = getSides(miit->meshset,mesh_bit_level,recursive,verb); CHKERRQ(ierr);
87  } else {
88  SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_FOUND,"msId is not there");
89  }
90  PetscFunctionReturn(0);
91 }
MoFEM::ierr
static PetscErrorCode ierr
Definition: Common.hpp:26
MoFEM::Interface::query_interface
PetscErrorCode query_interface(IFace *&ptr) const
Definition: Interface.hpp:47
CHKERRQ
CHKERRQ(ierr)
MoFEM::PrismInterface::getSides
virtual PetscErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=-1)
create two children meshsets in the meshset containing tetrahedral on two sides of faces ...
Definition: PrismInterface.cpp:74
MoFEM::Interface
InterfaceThis interface is used by user to:
Definition: Interface.hpp:38

◆ getSides() [2/2]

PetscErrorCode MoFEM::PrismInterface::getSides ( const EntityHandle  sideset,
const BitRefLevel  mesh_bit_level,
const bool  recursive,
int  verb = -1 
)
virtual

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 eges on the boundary of the face which is in the voulume 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 93 of file PrismInterface.cpp.

93  {
94 
95 
96  MoFEM::Interface &m_field = cOre;
97  moab::Interface &moab = m_field.get_moab();
98  PetscFunctionBegin;
99  Range mesh_level_ents3d,mesh_level_ents3d_tris;
100  Range mesh_level_tris;
101  Range mesh_level_edges;
102  Range mesh_level_nodes;
103  if(mesh_bit_level.any()) {
104  ierr = m_field.get_entities_by_type_and_ref_level(mesh_bit_level,BitRefLevel().set(),MBTET,mesh_level_ents3d); CHKERRQ(ierr);
105  ierr = m_field.get_entities_by_type_and_ref_level(mesh_bit_level,BitRefLevel().set(),MBTRI,mesh_level_tris); CHKERRQ(ierr);
106  ierr = m_field.get_entities_by_type_and_ref_level(mesh_bit_level,BitRefLevel().set(),MBEDGE,mesh_level_edges); CHKERRQ(ierr);
107  ierr = m_field.get_entities_by_type_and_ref_level(mesh_bit_level,BitRefLevel().set(),MBVERTEX,mesh_level_nodes); CHKERRQ(ierr);
108  rval = moab.get_adjacencies(mesh_level_ents3d,2,false,mesh_level_ents3d_tris,moab::Interface::UNION); CHKERRQ_MOAB(rval);
109 
110  }
111  Range mesh_level_prisms;
112  if(mesh_bit_level.any()) {
113  ierr = m_field.get_entities_by_type_and_ref_level(mesh_bit_level,BitRefLevel().set(),MBPRISM,mesh_level_prisms); CHKERRQ(ierr);
114  mesh_level_ents3d.merge(mesh_level_prisms);
115  }
116  Skinner skin(&moab);
117  //get interface triangles from side set
118  Range triangles;
119  rval = moab.get_entities_by_type(sideset,MBTRI,triangles,recursive); CHKERRQ_MOAB(rval);
120  if(mesh_bit_level.any()) {
121  triangles = intersect(triangles,mesh_level_ents3d_tris);
122  }
123  if(verb>1) {
124  PetscPrintf(m_field.get_comm(),"Nb. of triangles in set %u\n",triangles.size());
125  }
126  //get nodes, edges and 3d ents (i.e. tets and prisms)
127  Range nodes; // nodes from triangles
128  rval = moab.get_connectivity(triangles,nodes,true); CHKERRQ_MOAB(rval);
129  Range ents3d,ents3d_with_prisms; // 3d ents form nodes
130  rval = moab.get_adjacencies(nodes,3,false,ents3d_with_prisms,moab::Interface::UNION); CHKERRQ_MOAB(rval);
131  if(mesh_bit_level.any()) {
132  ents3d_with_prisms = intersect(ents3d_with_prisms,mesh_level_ents3d);
133  }
134  ents3d = ents3d_with_prisms.subset_by_type(MBTET); // take only tets, add prism later
135  //take skin faces
136  Range skin_faces; // skin faces from 3d ents
137  rval = skin.find_skin(0,ents3d,false,skin_faces); CHKERRQ_MOAB(rval);
138  //take skin edges (boundary of surface if there is any)
139  Range skin_edges_boundary; //skin edges from triangles
140  rval = skin.find_skin(0,triangles,false,skin_edges_boundary); CHKERRQ_MOAB(rval);
141  if(verb>3) PetscPrintf(m_field.get_comm(),"skin_edges_boundary %u\n",skin_edges_boundary.size());
142  //take all edges on skin faces (i.e. skin surface)
143  Range skin_faces_edges; //edges from skin faces of 3d ents
144  rval = moab.get_adjacencies(skin_faces,1,false,skin_faces_edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
145  if(mesh_bit_level.any()) {
146  skin_faces_edges = intersect(skin_faces_edges,mesh_level_edges);
147  }
148  if(verb>3) PetscPrintf(m_field.get_comm(),"skin_faces_edges %u\n",skin_faces_edges.size());
149  //note: that skin faces edges do not contain internal boundary
150  //note: that prisms are not included in ents3d, so if ents3d have border with other inteface is like external boundary
151  //skin edges boundary are internal edge <- skin_faces_edges contains edges which are on the body boundary <- that is the trick
152  skin_edges_boundary = subtract(skin_edges_boundary,skin_faces_edges); // from skin edges subtract edges from skin faces of 3d ents (only internal edges)
153  if(verb>3) {
154  EntityHandle out_meshset;
155  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
156  rval = moab.add_entities(out_meshset,triangles); CHKERRQ_MOAB(rval);
157  rval = moab.write_file("triangles.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
158  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
159  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
160  rval = moab.add_entities(out_meshset,ents3d); CHKERRQ_MOAB(rval);
161  rval = moab.write_file("ents3d.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
162  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
163  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
164  rval = moab.add_entities(out_meshset,skin_edges_boundary); CHKERRQ_MOAB(rval);
165  rval = moab.write_file("skin_edges_boundary.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
166  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
167  }
168  if(verb>3) PetscPrintf(m_field.get_comm(),"subtract skin_edges_boundary %u\n",skin_edges_boundary.size());
169  //Get nodes on boundary edge
170  Range skin_nodes_boundary;
171  rval = moab.get_connectivity(skin_edges_boundary,skin_nodes_boundary,true); CHKERRQ_MOAB(rval);
172  //Remove noded which are boundary with other existing interface
173  Range prisms_nodes;
174  rval = moab.get_connectivity(ents3d_with_prisms.subset_by_type(MBPRISM),prisms_nodes,true); CHKERRQ_MOAB(rval);
175  skin_nodes_boundary = subtract(skin_nodes_boundary,prisms_nodes);
176  if(verb>3) PetscPrintf(m_field.get_comm(),"subtract skin_nodes_boundary %u\n",skin_nodes_boundary.size());
177  //use nodes on body boundary and interface (without internal boundary) to find adjacent tets
178  Range nodes_without_front = subtract(nodes,skin_nodes_boundary); // nodes_without_front adjacent to all split face edges except those on internal edge
179  if(verb>3) PetscPrintf(m_field.get_comm(),"adj. node if ents3d but not on the internal edge %u\n",nodes_without_front.size());
180  //ents3 that are adjacent to front nodes on split faces but not those which are on the front nodes on internal edge
181  ents3d.clear();
182  ents3d_with_prisms.clear();
183  rval = moab.get_adjacencies(
184  nodes_without_front,3,false,ents3d_with_prisms,moab::Interface::UNION
185  ); CHKERRQ_MOAB(rval);
186  if(mesh_bit_level.any()) {
187  ents3d_with_prisms = intersect(ents3d_with_prisms,mesh_level_ents3d);
188  }
189  ents3d = ents3d_with_prisms.subset_by_type(MBTET);
190  if(verb>3) PetscPrintf(m_field.get_comm(),"adj. ents3d to fornt nodes %u\n",ents3d.size());
191  Range side_ents3d;
192  unsigned int nb_side_ents3d = side_ents3d.size();
193  side_ents3d.insert(*ents3d.begin());
194  Range side_ents3d_tris_on_surface;
195 
196  do {
197  do {
198  Range adj_tris,adj_ents3d;
199  nb_side_ents3d = side_ents3d.size();
200  if(verb>2) PetscPrintf(m_field.get_comm(),"nb_side_ents3d %u\n",nb_side_ents3d);
201  //get faces
202  rval = moab.get_adjacencies(
203  side_ents3d.subset_by_type(MBTET),2,false,adj_tris,moab::Interface::UNION
204  ); CHKERRQ_MOAB(rval);
205  if(mesh_bit_level.any()) {
206  adj_tris = intersect(adj_tris,mesh_level_tris);
207  }
208  //subtrace from faces interface
209  adj_tris = subtract(adj_tris,triangles);
210  if(verb>2) PetscPrintf(m_field.get_comm(),"adj_tris %u\n",adj_tris.size());
211  //get tets adjacent to faces
212  rval = moab.get_adjacencies(
213  adj_tris,3,true,adj_ents3d,moab::Interface::UNION
214  ); CHKERRQ_MOAB(rval);
215  //intersect tets with tets adjacent to inetface
216  adj_ents3d = intersect(adj_ents3d,ents3d_with_prisms);
217  if(verb>2) PetscPrintf(m_field.get_comm(),"adj_ents3d %u\n",adj_ents3d.size());
218  //add tets to side
219  side_ents3d.insert(adj_ents3d.begin(),adj_ents3d.end());
220  } while (nb_side_ents3d != side_ents3d.size());
221 
222  Range side_ents3d_tris;
223  rval = moab.get_adjacencies(
224  side_ents3d,2,false,side_ents3d_tris,moab::Interface::UNION
225  ); CHKERRQ_MOAB(rval);
226  side_ents3d_tris_on_surface = intersect(side_ents3d_tris,triangles);
227  if(side_ents3d_tris_on_surface.size()!=triangles.size()) {
228  Range left_triangles = subtract(triangles,side_ents3d_tris_on_surface);
229  Range tets;
230  rval = moab.get_adjacencies(&*left_triangles.begin(),1,3,false,tets); CHKERRQ_MOAB(rval);
231  tets = intersect(tets,ents3d_with_prisms);
232  if(tets.empty()) {
233  SETERRQ(
234  m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,
235  "not all faces on surface going to be split"
236  );
237  }
238  side_ents3d.insert(*tets.begin());
239  }
240 
241  } while (side_ents3d_tris_on_surface.size()!=triangles.size());
242 
243  if(ents3d_with_prisms.size() == side_ents3d.size()) {
244  SETERRQ(
245  m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,
246  "all tets on one side, no-interface"
247  );
248  }
249  //other side ents
250  Range other_side = subtract(ents3d_with_prisms,side_ents3d);
251  //side nodes
252  Range side_nodes;
253  rval = moab.get_connectivity(side_ents3d.subset_by_type(MBTET),side_nodes,true); CHKERRQ_MOAB(rval);
254  //nodes on crack surface without front
255  nodes_without_front = intersect(nodes_without_front,side_nodes);
256  Range side_edges;
257  rval = moab.get_adjacencies(side_ents3d.subset_by_type(MBTET),1,false,side_edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
258  skin_edges_boundary = intersect(skin_edges_boundary,side_edges);
259  //make child meshsets
260  std::vector<EntityHandle> children;
261  rval = moab.get_child_meshsets(sideset,children); CHKERRQ_MOAB(rval);
262  if(children.empty()) {
263  children.resize(3);
264  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,children[0]); CHKERRQ_MOAB(rval);
265  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,children[1]); CHKERRQ_MOAB(rval);
266  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,children[2]); CHKERRQ_MOAB(rval);
267  rval = moab.add_child_meshset(sideset,children[0]); CHKERRQ_MOAB(rval);
268  rval = moab.add_child_meshset(sideset,children[1]); CHKERRQ_MOAB(rval);
269  rval = moab.add_child_meshset(sideset,children[2]); CHKERRQ_MOAB(rval);
270  } else {
271  if(children.size()!=3) {
272  SETERRQ(PETSC_COMM_SELF,1,"this meshset shuld have 3 children meshsets");
273  }
274  children.resize(3);
275  ierr = moab.clear_meshset(&children[0],3); CHKERRQ(ierr);
276  }
277  EntityHandle &child_side = children[0];
278  EntityHandle &child_other_side = children[1];
279  EntityHandle &child_nodes_and_skin_edges = children[2];
280  rval = moab.add_entities(child_side,side_ents3d); CHKERRQ_MOAB(rval);
281  rval = moab.add_entities(child_other_side,other_side); CHKERRQ_MOAB(rval);
282  rval = moab.add_entities(child_nodes_and_skin_edges,nodes_without_front); CHKERRQ_MOAB(rval);
283  rval = moab.add_entities(child_nodes_and_skin_edges,skin_edges_boundary); CHKERRQ_MOAB(rval);
284  if(verb>1) {
285  PetscPrintf(m_field.get_comm(),"Nb. of side ents3d in set %u\n",side_ents3d.size());
286  PetscPrintf(m_field.get_comm(),"Nb. of other side ents3d in set %u\n",other_side.size());
287  PetscPrintf(m_field.get_comm(),"Nb. of boundary edges %u\n",skin_edges_boundary.size());
288  }
289  if(verb>3) {
290  ierr = moab.write_file("side.vtk","VTK","",&children[0],1); CHKERRQ(ierr);
291  ierr = moab.write_file("other_side.vtk","VTK","",&children[1],1); CHKERRQ(ierr);
292  }
293  PetscFunctionReturn(0);
294 }
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:362
MoFEM::ierr
static PetscErrorCode ierr
Definition: Common.hpp:26
MoFEM::rval
static MoABErrorCode rval
Definition: Common.hpp:25
MoFEM::Interface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::Interface::get_entities_by_type_and_ref_level
virtual PetscErrorCode get_entities_by_type_and_ref_level(const BitRefLevel &bit, const BitRefLevel &mask, const EntityType type, const EntityHandle meshset, int verb=-1)=0
add all ents from ref level given by bit to meshset
MoFEM::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:48
CHKERRQ
CHKERRQ(ierr)
MoFEM::Interface
InterfaceThis interface is used by user to:
Definition: Interface.hpp:38
MoFEM::Interface::get_comm
virtual MPI_Comm & get_comm() const =0

◆ queryInterface()

PetscErrorCode MoFEM::PrismInterface::queryInterface ( const MOFEMuuid &  uuid,
UnknownInterface **  iface 
)
virtual

Implements MoFEM::UnknownInterface.

Definition at line 55 of file PrismInterface.cpp.

55  {
56  PetscFunctionBegin;
57  *iface = NULL;
58  if(uuid == IDD_MOFEMPrismInterface) {
59  *iface = dynamic_cast<PrismInterface*>(this);
60  PetscFunctionReturn(0);
61  }
62  if(uuid == IDD_MOFEMUnknown) {
63  *iface = dynamic_cast<UnknownInterface*>(this);
64  PetscFunctionReturn(0);
65  }
66  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"unknown interface");
67  PetscFunctionReturn(0);
68 }
MoFEM::IDD_MOFEMUnknown
static const MOFEMuuid IDD_MOFEMUnknown
Definition: UnknownInterface.hpp:46
MoFEM::PrismInterface::PrismInterface
PrismInterface(const MoFEM::Core &core)
Definition: PrismInterface.cpp:70
MoFEM::IDD_MOFEMPrismInterface
static const MOFEMuuid IDD_MOFEMPrismInterface
Definition: PrismInterface.hpp:26

◆ splitSides() [1/3]

PetscErrorCode MoFEM::PrismInterface::splitSides ( const EntityHandle  meshset,
const BitRefLevel &  bit,
const int  msId,
const CubitBCType  cubit_bc_type,
const bool  add_iterfece_entities,
const bool  recursive = false,
int  verb = -1 
)
virtual

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_intefece_entitiesmeshset which contain the interface
recursiveif true parent meshset is searched recursively

Each inteface face has two tages, 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 abot inteface side, second tag inform about side adjacent inteface element.

Definition at line 295 of file PrismInterface.cpp.

298  {
299 
300  MoFEM::Interface &m_field = cOre;
301  MeshsetsManager *meshsets_manager_ptr;
302  PetscFunctionBegin;
303  ierr = m_field.query_interface(meshsets_manager_ptr); CHKERRQ(ierr);
304  CubitMeshSet_multiIndex::index<Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator
305  miit = meshsets_manager_ptr->getMeshsetsMultindex().get<Composite_Cubit_msId_And_MeshSetType_mi_tag>().find(boost::make_tuple(msId,cubit_bc_type.to_ulong()));
306  if(miit!=meshsets_manager_ptr->getMeshsetsMultindex().get<Composite_Cubit_msId_And_MeshSetType_mi_tag>().end()) {
307  ierr = splitSides(
308  meshset,bit,miit->meshset,add_iterfece_entities,recursive,verb); CHKERRQ(ierr);
309  } else {
310  SETERRQ(PETSC_COMM_SELF,1,"msId is not there");
311  }
312  PetscFunctionReturn(0);
313 }
MoFEM::ierr
static PetscErrorCode ierr
Definition: Common.hpp:26
MoFEM::Interface::query_interface
PetscErrorCode query_interface(IFace *&ptr) const
Definition: Interface.hpp:47
CHKERRQ
CHKERRQ(ierr)
MoFEM::Interface
InterfaceThis interface is used by user to:
Definition: Interface.hpp:38
MoFEM::PrismInterface::splitSides
virtual PetscErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_iterfece_entities, const bool recursive=false, int verb=-1)
split nodes and other entities of tetrahedral in children sets and add prism elements ...
Definition: PrismInterface.cpp:295

◆ splitSides() [2/3]

PetscErrorCode MoFEM::PrismInterface::splitSides ( const EntityHandle  meshset,
const BitRefLevel &  bit,
const EntityHandle  sideset,
const bool  add_iterfece_entities,
const bool  recursive = false,
int  verb = -1 
)
virtual

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 314 of file PrismInterface.cpp.

318  {
319 
320  PetscFunctionBegin;
321  ierr = splitSides(meshset,bit,
322  BitRefLevel(),BitRefLevel(),sideset,add_iterfece_entities,recursive,verb
323  ); CHKERRQ(ierr);
324  PetscFunctionReturn(0);
325 }
MoFEM::ierr
static PetscErrorCode ierr
Definition: Common.hpp:26
MoFEM::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:48
CHKERRQ
CHKERRQ(ierr)
MoFEM::PrismInterface::splitSides
virtual PetscErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_iterfece_entities, const bool recursive=false, int verb=-1)
split nodes and other entities of tetrahedral in children sets and add prism elements ...
Definition: PrismInterface.cpp:295

◆ splitSides() [3/3]

PetscErrorCode MoFEM::PrismInterface::splitSides ( const EntityHandle  meshset,
const BitRefLevel &  bit,
const BitRefLevel &  inheret_from_bit_level,
const BitRefLevel &  inheret_from_bit_level_mask,
const EntityHandle  sideset,
const bool  add_iterfece_entities,
const bool  recursive = false,
int  verb = -1 
)
virtual

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
inheret_from_bit_levelinheret nodes and other entities form this bit level.
add_iterfece_entitiesadd prism elements at interface
recuslsivedo meshesets in the meshset

note inheret_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 326 of file PrismInterface.cpp.

331  {
332 
333 
334  MoFEM::Interface &m_field = cOre;
335  moab::Interface &moab = m_field.get_moab();
336  PetscFunctionBegin;
337  std::vector<EntityHandle> children;
338  //get children meshsets
339  rval = moab.get_child_meshsets(sideset,children); CHKERRQ_MOAB(rval);
340  if(children.size()!=3) {
341  SETERRQ(
342  m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,
343  "should be 3 child meshsets, each of them contains tets on two sides of interface"
344  );
345  }
346  //3d ents on "father" side
347  Range side_ents3d;
348  rval = moab.get_entities_by_handle(children[0],side_ents3d,false); CHKERRQ_MOAB(rval);
349  //3d ents on "mather" side
350  Range other_ents3d;
351  rval = moab.get_entities_by_handle(children[1],other_ents3d,false); CHKERRQ_MOAB(rval);
352  //faces of interface
353  Range triangles;
354  rval = moab.get_entities_by_type(sideset,MBTRI,triangles,recursive); CHKERRQ_MOAB(rval);
355  Range side_ents3d_tris;
356  rval = moab.get_adjacencies(side_ents3d,2,true,side_ents3d_tris,moab::Interface::UNION); CHKERRQ_MOAB(rval);
357  triangles = intersect(triangles,side_ents3d_tris);
358  //nodes on interface but not on crack front (those should not be splitted)
359  Range nodes;
360  rval = moab.get_entities_by_type(children[2],MBVERTEX,nodes,false); CHKERRQ_MOAB(rval);
361  Range meshset_3d_ents,meshset_2d_ents;
362  rval = moab.get_entities_by_dimension(meshset,3,meshset_3d_ents,true); CHKERRQ_MOAB(rval);
363  Range meshset_tets = meshset_3d_ents.subset_by_type(MBTET);
364  rval = moab.get_adjacencies(meshset_tets,2,false,meshset_2d_ents,moab::Interface::UNION); CHKERRQ_MOAB(rval);
365  side_ents3d = intersect(meshset_3d_ents,side_ents3d);
366  other_ents3d = intersect(meshset_3d_ents,other_ents3d);
367  triangles = intersect(meshset_2d_ents,triangles);
368  if(verb>3) {
369  PetscPrintf(m_field.get_comm(),"triangles %u\n",triangles.size());
370  PetscPrintf(m_field.get_comm(),"side_ents3d %u\n",side_ents3d.size());
371  PetscPrintf(m_field.get_comm(),"nodes %u\n",nodes.size());
372  }
373 
374  const RefEntity_multiIndex *refined_ents_ptr;
375  ierr = m_field.get_ref_ents(&refined_ents_ptr); CHKERRQ(ierr);
376  typedef RefEntity_multiIndex::index<Ent_mi_tag>::type RefEntsByEntType;
377  const RefEntsByEntType &ref_ents_by_type = refined_ents_ptr->get<Ent_mi_tag>();
378  RefEntity_multiIndex_view_by_parent_entity ref_parent_ents_view;
379  //create view index by parent entity
380  {
381  typedef RefEntity_multiIndex::index<Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
382  const RefEntsByComposite &ref_ents =
383  refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
384 
385  RefEntsByComposite::iterator miit;
386  RefEntsByComposite::iterator hi_miit;
387  //view by parent type (VERTEX)
388  miit = ref_ents.lower_bound(boost::make_tuple(MBVERTEX,MBVERTEX));
389  hi_miit = ref_ents.upper_bound(boost::make_tuple(MBVERTEX,MBVERTEX));
390  for(;miit!=hi_miit;miit++) {
391  if(((*miit)->getBitRefLevel()&inheret_from_bit_level_mask) == (*miit)->getBitRefLevel()) {
392  if(((*miit)->getBitRefLevel()&inheret_from_bit_level).any()) {
393  std::pair<RefEntity_multiIndex_view_by_parent_entity::iterator,bool> p_ref_ent_view;
394  p_ref_ent_view = ref_parent_ents_view.insert(*miit);
395  if(!p_ref_ent_view.second) {
396  SETERRQ(PETSC_COMM_SELF,1,"non uniqe insertion");
397  }
398  }
399  }
400  }
401  }
402  //maps nodes on "father" and "mather" side
403  std::map<
404  EntityHandle, /*node on "mather" side*/
405  EntityHandle /*node on "father" side*/
406  > map_nodes;
407  //add new nodes on interface and create map
408  Range::iterator nit = nodes.begin();
409  double coord[3];
410  for(;nit!=nodes.end();nit++) {
411  //find ref enet
412  RefEntsByEntType::iterator miit_ref_ent = ref_ents_by_type.find(*nit);
413  if(miit_ref_ent == ref_ents_by_type.end()) {
414  SETERRQ(PETSC_COMM_SELF,1,"can not find node in MoFEM database");
415  }
416  EntityHandle child_entity = 0;
417  RefEntity_multiIndex::iterator child_it;
418  RefEntity_multiIndex_view_by_parent_entity::iterator child_iit;
419  child_iit = ref_parent_ents_view.find(*nit);
420  if(child_iit != ref_parent_ents_view.end()) {
421  child_it = refined_ents_ptr->find((*child_iit)->getRefEnt());
422  BitRefLevel bit_child = (*child_it)->getBitRefLevel();
423  if( (inheret_from_bit_level&bit_child).none() ) {
424  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"data inconsistency");
425  }
426  child_entity = (*child_it)->getRefEnt();
427  }
428  //
429  bool success;
430  if(child_entity == 0) {
431  rval = moab.get_coords(&*nit,1,coord); CHKERRQ_MOAB(rval);
432  EntityHandle new_node;
433  rval = moab.create_vertex(coord,new_node); CHKERRQ_MOAB(rval);
434  map_nodes[*nit] = new_node;
435  //create new node on "father" side
436  //parent is node on "mather" side
437  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(),&new_node,1,&*nit); CHKERRQ_MOAB(rval);
438  std::pair<RefEntity_multiIndex::iterator,bool> p_ref_ent =
439  const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->insert(
440  boost::shared_ptr<RefEntity>
441  (new RefEntity(m_field.get_basic_entity_data_ptr(),new_node))
442  );
443  //set ref bit level to node on "father" side
444  success = const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->
445  modify(p_ref_ent.first,RefEntity_change_add_bit(bit));
446  if(!success) SETERRQ(m_field.get_comm(),MOFEM_OPERATION_UNSUCCESSFUL,"modification unsuccessful");
447  } else {
448  map_nodes[*nit] = child_entity;
449  //set ref bit level to node on "father" side
450  success = const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->
451  modify(child_it,RefEntity_change_add_bit(bit));
452  if(!success) SETERRQ(m_field.get_comm(),MOFEM_OPERATION_UNSUCCESSFUL,"modification unsuccessful");
453  }
454  //set ref bit level to node on "mather" side
455  success = const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->
456  modify(miit_ref_ent,RefEntity_change_add_bit(bit));
457  if(!success) SETERRQ(m_field.get_comm(),MOFEM_OPERATION_UNSUCCESSFUL,"modification unsuccessful");
458  }
459  //crete meshset for new mesh bit level
460  EntityHandle meshset_for_bit_level;
461  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,meshset_for_bit_level); CHKERRQ_MOAB(rval);
462  //subtract those elements which will be refined, i.e. disconnected form other side elements, and connected to new prisms, if they area created
463  meshset_3d_ents = subtract(meshset_3d_ents,side_ents3d);
464  rval = moab.add_entities(meshset_for_bit_level,meshset_3d_ents); CHKERRQ_MOAB(rval);
465  for(int dd = 0;dd<3;dd++) {
466  Range ents_dd;
467  rval = moab.get_adjacencies(meshset_3d_ents,dd,false,ents_dd,moab::Interface::UNION); CHKERRQ_MOAB(rval);
468  rval = moab.add_entities(meshset_for_bit_level,ents_dd); CHKERRQ_MOAB(rval);
469  }
470  //
471  //typedef RefEntity_multiIndex::index<Composite_ParentEnt_And_EntType_mi_tag>::type ref_ent_by_composite;
472  //ref_ent_by_composite &by_composite = refinedEntities.get<Composite_ParentEnt_And_EntType_mi_tag>();
473  //create new 3d ents on "father" side
474  Range new_3d_ents;
475  Range::iterator eit3d = side_ents3d.begin();
476  for(;eit3d!=side_ents3d.end();eit3d++) {
477  RefEntsByEntType::iterator miit_ref_ent = ref_ents_by_type.find(*eit3d);
478  if(miit_ref_ent==ref_ents_by_type.end()) {
479  SETERRQ(m_field.get_comm(),MOFEM_OPERATION_UNSUCCESSFUL,"tet not in database");
480  }
481  int num_nodes;
482  const EntityHandle* conn;
483  rval = moab.get_connectivity(*eit3d,conn,num_nodes,true); CHKERRQ_MOAB(rval);
484  EntityHandle new_conn[num_nodes];
485  int nb_new_conn = 0;
486  int ii = 0;
487  for(; ii<num_nodes; ii++) {
488  std::map<EntityHandle,EntityHandle>::iterator mit = map_nodes.find(conn[ii]);
489  if(mit != map_nodes.end()) {
490  new_conn[ii] = mit->second;
491  nb_new_conn++;
492  if(verb>6) {
493  PetscPrintf(m_field.get_comm(),"nodes %u -> %d\n",conn[ii],new_conn[ii]);
494  }
495  } else {
496  new_conn[ii] = conn[ii];
497  }
498  }
499  if(nb_new_conn==0) {
500  if(verb>3) {
501  EntityHandle meshset_error_out;
502  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,meshset_error_out); CHKERRQ_MOAB(rval);
503  rval = moab.add_entities(meshset_error_out,&*eit3d,1); CHKERRQ_MOAB(rval);
504  ierr = moab.write_file("error_out.vtk","VTK","",&meshset_error_out,1); CHKERRQ(ierr);
505  rval = moab.delete_entities(&meshset_error_out,1); CHKERRQ_MOAB(rval);
506 
507  }
508  SETERRQ1(PETSC_COMM_SELF,1,"database inconsistency, in side_ent3 is a tet which has no common node with interface, num_nodes = %d",num_nodes);
509  }
510 
511  const RefElement_multiIndex *refined_finite_elements_ptr;
512  ierr = m_field.get_ref_finite_elements(&refined_finite_elements_ptr); CHKERRQ(ierr);
513 
514  //here is created new or prism is on interface
515  EntityHandle existing_ent = 0;
516  /* check if tet element whith new connectivity is in database*/
517  RefElement_multiIndex::index<Ent_Ent_mi_tag>::type::iterator child_iit,hi_child_iit;
518  child_iit = refined_finite_elements_ptr->get<Ent_Ent_mi_tag>().lower_bound(*eit3d);
519  hi_child_iit = refined_finite_elements_ptr->get<Ent_Ent_mi_tag>().upper_bound(*eit3d);
520  for(;child_iit!=hi_child_iit;child_iit++) {
521  const EntityHandle* conn_ref_tet;
522  rval = moab.get_connectivity(child_iit->getRefEnt(),conn_ref_tet,num_nodes,true); CHKERRQ_MOAB(rval);
523  int nn = 0;
524  for(;nn<num_nodes;nn++) {
525  if(conn_ref_tet[nn]!=new_conn[nn]) {
526  break;
527  }
528  }
529  if(nn == num_nodes) {
530  if(existing_ent != 0) {
531  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"database inconsistency");
532  }
533  existing_ent = child_iit->getRefEnt();
534  }
535  }
536  switch (moab.type_from_handle(*eit3d)) {
537  case MBTET: {
538  RefEntsByEntType::iterator child_it;
539  EntityHandle tet;
540  if(existing_ent == 0) {
541  Range new_conn_tet;
542  rval = moab.get_adjacencies(new_conn,4,3,false,new_conn_tet); CHKERRQ_MOAB(rval);
543  if(new_conn_tet.empty()) {
544  rval = moab.create_element(MBTET,new_conn,4,tet); CHKERRQ_MOAB(rval);
545  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(),&tet,1,&*eit3d); CHKERRQ_MOAB(rval);
546  } else {
547  RefElement_multiIndex::index<Ent_mi_tag>::type::iterator rit,new_rit;
548  rit = refined_finite_elements_ptr->get<Ent_mi_tag>().find(*eit3d);
549  if(rit==refined_finite_elements_ptr->get<Ent_mi_tag>().end()) {
550  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"can't find this in database");
551  }
552  new_rit = refined_finite_elements_ptr->get<Ent_mi_tag>().find(*new_conn_tet.begin());
553  if(new_rit==refined_finite_elements_ptr->get<Ent_mi_tag>().end()) {
554  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"can't find this in database");
555  }
556  tet = *new_conn_tet.begin();
557  /*std::ostringstream ss;
558  ss << "nb new conns: " << nb_new_conn << std::endl;
559  ss << "new_conn_tets.size() " << new_conn_tet.size() << std::endl;
560  ss << "data inconsistency\n";
561  ss << "this ent:\n";
562  ss << *rit->ref_ptr << std::endl;
563  ss << "found this ent:\n";
564  ss << *new_rit->ref_ptr << std::endl;
565  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());*/
566  }
567  } else {
568  tet = existing_ent;
569  }
570  rval = moab.add_entities(meshset_for_bit_level,&tet,1); CHKERRQ_MOAB(rval);
571  rval = moab.add_entities(meshset_for_bit_level,new_conn,4); CHKERRQ_MOAB(rval);
572  new_3d_ents.insert(tet);
573  } break;
574  case MBPRISM: {
575  EntityHandle prism;
576  if(verb>3) {
577  PetscPrintf(m_field.get_comm(),"prims nb_new_nodes %d\n",nb_new_conn);
578  }
579  if(existing_ent == 0) {
580  Range new_conn_prism;
581  rval = moab.get_adjacencies(new_conn,6,3,false,new_conn_prism); CHKERRQ_MOAB(rval);
582  if(new_conn_prism.empty()) {
583  rval = moab.create_element(MBPRISM,new_conn,6,prism); CHKERRQ_MOAB(rval);
584  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(),&prism,1,&*eit3d); CHKERRQ_MOAB(rval);
585  } else {
586  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"data inconsistency");
587  }
588  } else {
589  prism = existing_ent;
590  }
591  rval = moab.add_entities(meshset_for_bit_level,&prism,1); CHKERRQ_MOAB(rval);
592  rval = moab.add_entities(meshset_for_bit_level,new_conn,4); CHKERRQ_MOAB(rval);
593  new_3d_ents.insert(prism);
594  } break;
595  default:
596  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"not implemented");
597  }
598  }
599  Range new_ents;
600  //create new entities by adjecies form new tets
601  rval = moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET),2,true,new_ents,moab::Interface::UNION); CHKERRQ_MOAB(rval);
602  rval = moab.get_adjacencies(new_3d_ents.subset_by_type(MBTET),1,true,new_ents,moab::Interface::UNION); CHKERRQ_MOAB(rval);
603  //Tags for setting side
604  Tag th_interface_side;
605  const int def_side[] = {0};
606  rval = moab.tag_get_handle("INTERFACE_SIDE",1,MB_TYPE_INTEGER,
607  th_interface_side,MB_TAG_CREAT|MB_TAG_SPARSE,def_side); CHKERRQ_MOAB(rval);
608  //add new edges and triangles to mofem database
609  Range ents;
610  rval = moab.get_adjacencies(triangles,1,false,ents,moab::Interface::UNION); CHKERRQ_MOAB(rval);
611  ents.insert(triangles.begin(),triangles.end());
612  Range new_ents_in_database; //this range contains all new entities
613  Range::iterator eit = ents.begin();
614  for(;eit!=ents.end();eit++) {
615  int num_nodes;
616  const EntityHandle* conn;
617  rval = moab.get_connectivity(*eit,conn,num_nodes,true); CHKERRQ_MOAB(rval);
618  EntityHandle new_conn[num_nodes];
619  int nb_new_conn = 0;
620  int ii = 0;
621  for(;ii<num_nodes; ii++) {
622  std::map<EntityHandle,EntityHandle>::iterator mit = map_nodes.find(conn[ii]);
623  if(mit != map_nodes.end()) {
624  new_conn[ii] = mit->second;
625  nb_new_conn++;
626  if(verb>6) {
627  PetscPrintf(m_field.get_comm(),"nodes %u -> %d\n",conn[ii],new_conn[ii]);
628  }
629  } else {
630  new_conn[ii] = conn[ii];
631  }
632  }
633  if(nb_new_conn==0) continue;
634  RefEntsByEntType::iterator miit_ref_ent = ref_ents_by_type.find(*eit);
635  if(miit_ref_ent == ref_ents_by_type.end()) {
636  SETERRQ(PETSC_COMM_SELF,1,"this entity (edge or tri) should be already in database");
637  }
638  Range new_ent; //contains all entities (edges or triangles) added to mofem database
639  switch (moab.type_from_handle(*eit)) {
640  case MBTRI: {
641  //get entity based on its connectivity
642  rval = moab.get_adjacencies(new_conn,3,2,false,new_ent); CHKERRQ_MOAB(rval);
643  if(new_ent.size() != 1) SETERRQ(PETSC_COMM_SELF,1,"this tri should be in moab database");
644  int new_side = 1;
645  rval = moab.tag_set_data(th_interface_side,&*new_ent.begin(),1,&new_side); CHKERRQ_MOAB(rval);
646  if(verb>3) PetscPrintf(m_field.get_comm(),"new_ent %u\n",new_ent.size());
647  //add prism element
648  if(add_iterfece_entities) {
649  if(inheret_from_bit_level.any()) {
650  SETERRQ(PETSC_COMM_SELF,1,"not implemented for inheret_from_bit_level");
651  }
652  //set prism connectivity
653  EntityHandle prism_conn[6] = {
654  conn[0],conn[1],conn[2],
655  new_conn[0],new_conn[1],new_conn[2]
656  };
657  EntityHandle prism;
658  rval = moab.create_element(MBPRISM,prism_conn,6,prism); CHKERRQ_MOAB(rval);
659  ierr = cOre.addPrismToDatabase(prism,verb); CHKERRQ(ierr);
660  rval = moab.add_entities(meshset_for_bit_level,&prism,1); CHKERRQ_MOAB(rval);
661  }
662  } break;
663  case MBEDGE: {
664  rval = moab.get_adjacencies(new_conn,2,1,false,new_ent); CHKERRQ_MOAB(rval);
665  if(new_ent.size()!=1) {
666  if(m_field.get_comm_rank()==0) {
667  EntityHandle out_meshset;
668  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
669  rval = moab.add_entities(out_meshset,&*eit,1); CHKERRQ_MOAB(rval);
670  rval = moab.write_file("debug_splitSides.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
671  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
672  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
673  rval = moab.add_entities(out_meshset,side_ents3d); CHKERRQ_MOAB(rval);
674  rval = moab.write_file("debug_splitSides_side_ents3d.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
675  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
676  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
677  rval = moab.add_entities(out_meshset,other_ents3d); CHKERRQ_MOAB(rval);
678  rval = moab.write_file("debug_splitSides_other_ents3d.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
679  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
680  rval = moab.create_meshset(MESHSET_SET|MESHSET_TRACK_OWNER,out_meshset); CHKERRQ_MOAB(rval);
681  rval = moab.add_entities(out_meshset,triangles); CHKERRQ_MOAB(rval);
682  rval = moab.write_file("debug_get_msId_3dENTS_split_triangles.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
683  rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
684  }
685  SETERRQ2(PETSC_COMM_SELF,1,
686  "this edge should be in moab database new_ent.size() = %u nb_new_conn = %d",
687  new_ent.size(),nb_new_conn);
688  }
689  } break;
690  default:
691  SETERRQ(PETSC_COMM_SELF,1,"Houston we have a problem !!!");
692  }
693  if(new_ent.size()!=1) {
694  SETERRQ1(PETSC_COMM_SELF,1,"new_ent.size() = %u, size always should be 1",new_ent.size());
695  }
696  //set parent
697  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(),&*new_ent.begin(),1,&*eit); CHKERRQ_MOAB(rval);
698  //add to database
699  std::pair<RefEntity_multiIndex::iterator,bool> p_ref_ent =
700  const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->
701  insert(boost::shared_ptr<RefEntity>(new RefEntity(m_field.get_basic_entity_data_ptr(),new_ent[0])));
702  const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->
703  modify(p_ref_ent.first,RefEntity_change_add_bit(bit));
704  new_ents_in_database.insert(new_ent.begin(),new_ent.end());
705  }
706  //all other entities, some ents like triangles and faces on the side of tets
707  Range side_adj_faces_and_edges;
708  rval = moab.get_adjacencies(side_ents3d.subset_by_type(MBTET),1,true,side_adj_faces_and_edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
709  rval = moab.get_adjacencies(side_ents3d.subset_by_type(MBTET),2,true,side_adj_faces_and_edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
710  //subtract entities already added to mofem database
711  side_adj_faces_and_edges = subtract(side_adj_faces_and_edges,new_ents_in_database);
712  eit = side_adj_faces_and_edges.begin();
713  for(;eit!=side_adj_faces_and_edges.end();eit++) {
714  int num_nodes;
715  const EntityHandle* conn;
716  rval = moab.get_connectivity(*eit,conn,num_nodes,true); CHKERRQ_MOAB(rval);
717  EntityHandle new_conn[num_nodes];
718  int nb_new_conn = 0;
719  int ii = 0;
720  for(;ii<num_nodes; ii++) {
721  std::map<EntityHandle,EntityHandle>::iterator mit = map_nodes.find(conn[ii]);
722  if(mit != map_nodes.end()) {
723  new_conn[ii] = mit->second;
724  nb_new_conn++;
725  if(verb>6) {
726  PetscPrintf(m_field.get_comm(),"nodes %u -> %d\n",conn[ii],new_conn[ii]);
727  }
728  } else {
729  new_conn[ii] = conn[ii];
730  }
731  }
732  if(nb_new_conn==0) continue;
733  RefEntsByEntType::iterator miit_ref_ent = ref_ents_by_type.find(*eit);
734  if(miit_ref_ent == ref_ents_by_type.end()) {
735  SETERRQ1(PETSC_COMM_SELF,1,"entity should be in MoFem database, num_nodes = %d",num_nodes);
736  }
737  Range new_ent;
738  switch (moab.type_from_handle(*eit)) {
739  case MBTRI: {
740  rval = moab.get_adjacencies(new_conn,3,2,false,new_ent); CHKERRQ_MOAB(rval);
741  }
742  break;
743  case MBEDGE: {
744  rval = moab.get_adjacencies(new_conn,2,1,false,new_ent); CHKERRQ_MOAB(rval);
745  }
746  break;
747  default:
748  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"Houston we have a problem");
749  }
750  if(new_ent.size()!=1) {
751  SETERRQ1(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"database inconsistency, new_ent.size() = %u",new_ent.size());
752  }
753  //add entity to mofem database
754  rval = moab.tag_set_data(cOre.get_th_RefParentHandle(),&*new_ent.begin(),1,&*eit); CHKERRQ_MOAB(rval);
755  std::pair<RefEntity_multiIndex::iterator,bool> p_ref_ent =
756  const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->
757  insert(
758  boost::shared_ptr<RefEntity>(new RefEntity(m_field.get_basic_entity_data_ptr(),new_ent[0]))
759  );
760  const_cast<RefEntity_multiIndex*>(refined_ents_ptr)->
761  modify(p_ref_ent.first,RefEntity_change_add_bit(bit));
762  if(verb>3) PetscPrintf(m_field.get_comm(),"new_ent %u\n",new_ent.size());
763  new_ents_in_database.insert(new_ent.begin(),new_ent.end());
764  }
765  //add new prisms which parents are part of other intefaces
766  Range new_3d_prims = new_3d_ents.subset_by_type(MBPRISM);
767  for(Range::iterator pit = new_3d_prims.begin();pit!=new_3d_prims.end();pit++) {
768  ierr = cOre.addPrismToDatabase(*pit,verb); CHKERRQ(ierr);
769  //get parent entity
770  EntityHandle parent_prism;
771  rval = moab.tag_get_data(cOre.get_th_RefParentHandle(),&*pit,1,&parent_prism); CHKERRQ_MOAB(rval);
772  const EntityHandle root_meshset = moab.get_root_set();
773  if(parent_prism == root_meshset) {
774  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"this prism should have parent");
775  }
776  if(moab.type_from_handle(parent_prism)!=MBPRISM) {
777  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"this prism should have parent which is prism as well");
778  }
779  int num_nodes;
780  //parent prism
781  const EntityHandle* conn_parent;
782  rval = moab.get_connectivity(parent_prism,conn_parent,num_nodes,true); MOAB_THROW(rval);
783  Range face_side3_parent,face_side4_parent;
784  rval = moab.get_adjacencies(conn_parent,3,2,false,face_side3_parent); CHKERRQ_MOAB(rval);
785  rval = moab.get_adjacencies(&conn_parent[3],3,2,false,face_side4_parent); CHKERRQ_MOAB(rval);
786  if(face_side3_parent.size()!=1) {
787  SETERRQ1(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"parent face3.size() = %u",face_side3_parent.size());
788  }
789  if(face_side4_parent.size()!=1) {
790  SETERRQ1(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"parent face4.size() = %u",face_side4_parent.size());
791  }
792  //new prism
793  const EntityHandle* conn;
794  rval = moab.get_connectivity(*pit,conn,num_nodes,true); MOAB_THROW(rval);
795  Range face_side3,face_side4;
796  rval = moab.get_adjacencies(conn,3,2,false,face_side3); CHKERRQ_MOAB(rval);
797  rval = moab.get_adjacencies(&conn[3],3,2,false,face_side4); CHKERRQ_MOAB(rval);
798  if(face_side3.size()!=1) {
799  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"face3 is missing");
800  }
801  if(face_side4.size()!=1) {
802  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"face4 is missing");
803  }
804  //
805  std::vector<EntityHandle> face(2),parent_face(2);
806  face[0] = *face_side3.begin();
807  face[1] = *face_side4.begin();
808  parent_face[0] = *face_side3_parent.begin();
809  parent_face[1] = *face_side4_parent.begin();
810  for(int ff = 0;ff<2;ff++) {
811  if(parent_face[ff] == face[ff]) continue;
812  int interface_side;
813  rval = moab.tag_get_data(th_interface_side,&parent_face[ff],1,&interface_side); CHKERRQ_MOAB(rval);
814  rval = moab.tag_set_data(th_interface_side,&face[ff],1,&interface_side); CHKERRQ_MOAB(rval);
815  EntityHandle parent_tri;
816  rval = moab.tag_get_data(cOre.get_th_RefParentHandle(),&face[ff],1,&parent_tri); CHKERRQ_MOAB(rval);
817  if(parent_tri != parent_face[ff]) {
818  SETERRQ1(PETSC_COMM_SELF,1,"wrong parent %lu",parent_tri);
819  }
820  if(new_ents_in_database.find(face[ff])==new_ents_in_database.end()) {
821  RefEntity_multiIndex::index<Ent_mi_tag>::type::iterator miit_ref_ent
822  = refined_ents_ptr->get<Ent_mi_tag>().find(face[ff]);
823  if(miit_ref_ent==refined_ents_ptr->get<Ent_mi_tag>().end()) {
824  SETERRQ(m_field.get_comm(),MOFEM_DATA_INCONSISTENCY,"this is not in database, but should not be");
825  }
826  }
827  }
828  }
829  //finalise by adding new tets and prism ti bitlelvel
830  ierr = m_field.query_interface<BitRefManager>()->setBitRefLevelByDim(meshset_for_bit_level,3,bit); CHKERRQ(ierr);
831  rval = moab.delete_entities(&meshset_for_bit_level,1); CHKERRQ_MOAB(rval);
832  ierr = moab.clear_meshset(&children[0],3); CHKERRQ(ierr);
833  PetscFunctionReturn(0);
834 }
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:362
MoFEM::Core::get_th_RefParentHandle
Tag get_th_RefParentHandle() const
Definition: Core.hpp:61
MoFEM::ierr
static PetscErrorCode ierr
Definition: Common.hpp:26
MoFEM::RefElement_multiIndex
multi_index_container< ptrWrapperRefElement, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityHandle, &ptrWrapperRefElement::getRefEnt > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityHandle, &ptrWrapperRefElement::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityType, &ptrWrapperRefElement::getEntType > >, ordered_non_unique< tag< Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag >, composite_key< ptrWrapperRefElement, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityHandle, &ptrWrapperRefElement::getParentEnt >, const_mem_fun< ptrWrapperRefElement::interface_type_RefElement, int, &ptrWrapperRefElement::getBitRefEdgesUlong > > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< ptrWrapperRefElement, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityHandle, &ptrWrapperRefElement::getRefEnt >, const_mem_fun< ptrWrapperRefElement::interface_type_RefEntity, EntityHandle, &ptrWrapperRefElement::getParentEnt > > > >> RefElement_multiIndex
Definition: FEMultiIndices.hpp:258
MoFEM::rval
static MoABErrorCode rval
Definition: Common.hpp:25
MoFEM::RefEntity_multiIndex_view_by_parent_entity
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< hashed_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_parent_entity
multi-index view of RefEntity by parent entity
Definition: EntsMultiIndices.hpp:456
MoFEM::Interface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::Core::addPrismToDatabase
PetscErrorCode addPrismToDatabase(const EntityHandle prism, int verb=-1)
Definition: Core.cpp:449
MOAB_THROW
#define MOAB_THROW(a)
Definition: definitions.h:376
MoFEM::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:48
MoFEM::Interface::query_interface
PetscErrorCode query_interface(IFace *&ptr) const
Definition: Interface.hpp:47
MoFEM::Interface::get_basic_entity_data_ptr
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
CHKERRQ
CHKERRQ(ierr)
dd
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:28
MoFEM::Interface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::Interface::get_ref_finite_elements
virtual PetscErrorCode get_ref_finite_elements(const RefElement_multiIndex **refined_finite_elements_ptr) const =0
Get ref finite elements multi-index form database.
MoFEM::RefEntity_multiIndex
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, EntityHandle,&RefEntity::getParentEnt >, const_mem_fun< RefEntity::BasicEntity, EntityType,&RefEntity::getEntType > > > >> RefEntity_multiIndex
Definition: EntsMultiIndices.hpp:436
MoFEM::Interface
InterfaceThis interface is used by user to:
Definition: Interface.hpp:38
MoFEM::Interface::get_ref_ents
virtual PetscErrorCode get_ref_ents(const RefEntity_multiIndex **refined_ents_ptr) const =0
Get ref entities multi-index from database.
MoFEM::Interface::get_comm
virtual MPI_Comm & get_comm() const =0

Member Data Documentation

◆ cOre

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

Definition at line 40 of file PrismInterface.hpp.

The documentation for this struct was generated from the following files:
Generated on Tue Sep 12 2017 13:32:37 by Doxygen 1.8.13 and hosted at University of Glasgow