v0.9.1
Functions | Variables
split_sideset.cpp File Reference

Split sidesets. More...

#include <MoFEM.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Detailed Description

Split sidesets.

Definition in file split_sideset.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
split_sideset.cpp.

Definition at line 26 of file split_sideset.cpp.

26  {
27  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
28 
29  try {
30 
31  // global variables
32  char mesh_file_name[255];
33  PetscBool flg_file = PETSC_FALSE;
34  PetscBool squash_bit_levels = PETSC_TRUE;
35  PetscBool flg_list_of_sidesets = PETSC_FALSE;
36  PetscBool output_vtk = PETSC_TRUE;
37 
38  int nb_sidesets = 10;
39  int sidesets[nb_sidesets];
40 
41  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Split sides options", "none");
42  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
43  mesh_file_name, 255, &flg_file);
44  CHKERR PetscOptionsBool("-squash_bit_levels", "squash bit levels", "",
45  squash_bit_levels, &squash_bit_levels, NULL);
46  CHKERR PetscOptionsIntArray("-side_sets", "get list of sidesets", "",
47  sidesets, &nb_sidesets, &flg_list_of_sidesets);
48  CHKERR PetscOptionsBool("-output_vtk", "if true outout vtk file", "",
49  output_vtk, &output_vtk, PETSC_NULL);
50  ierr = PetscOptionsEnd(); CHKERRG(ierr);
51 
52  moab::Core mb_instance;
53  moab::Interface &moab = mb_instance;
54  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
55  if (pcomm == NULL)
56  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
57  const char *option;
58  option = "";
59  CHKERR moab.load_file(mesh_file_name, 0, option);
60 
61  // Create MoFEM database
62  MoFEM::Core core(moab);
63  MoFEM::Interface &m_field = core;
64 
65  if (flg_file != PETSC_TRUE) {
66  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
67  "*** ERROR -my_file (MESH FILE NEEDED)");
68  }
69  if (flg_list_of_sidesets != PETSC_TRUE) {
70  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
71  "List of sidesets not given -my_side_sets ...");
72  }
73 
74  // Get interface to meshsets manager
75  MeshsetsManager *m_mng;
76  CHKERR m_field.getInterface(m_mng);
77  // Get interface for splitting manager
78  PrismInterface *interface_ptr;
79  CHKERR m_field.getInterface(interface_ptr);
80  BitRefManager *bit_mng;
81  CHKERR m_field.getInterface(bit_mng);
82 
83  // Seed mesh with bit levels
84  CHKERR bit_mng->setBitRefLevelByDim(0, 3, BitRefLevel().set(0));
85  std::vector<BitRefLevel> bit_levels;
86  bit_levels.push_back(BitRefLevel().set(0));
87 
88  typedef CubitMeshSet_multiIndex::index<
90  typedef CMeshsetByType::iterator CMIteratorByType;
91  typedef CubitMeshSet_multiIndex::index<
93  typedef CMeshsetByIdType::iterator CMIteratorByIdType;
94 
95  CubitMeshSet_multiIndex &meshsets_index = m_mng->getMeshsetsMultindex();
96  CMeshsetByType &m_by_type =
97  meshsets_index.get<CubitMeshSets_mask_meshset_mi_tag>();
98  CMeshsetByIdType &m_by_id_and_type =
99  meshsets_index.get<Composite_Cubit_msId_And_MeshSetType_mi_tag>();
100 
101  for (CMIteratorByType mit = m_by_type.lower_bound(SIDESET);
102  mit != m_by_type.upper_bound(SIDESET); mit++) {
103 
104  std::cout << "Sideset on the mesh id = " << mit->getMeshsetId()
105  << std::endl;
106  }
107 
108  // iterate sideset and split
109  for (int mm = 0; mm != nb_sidesets; mm++) {
110 
111  // find side set
112 
113  CMIteratorByIdType mit;
114  mit = m_by_id_and_type.find(boost::make_tuple(sidesets[mm], SIDESET));
115  if (mit == m_by_id_and_type.end()) {
116  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
117  "No sideset in database id = %d", sidesets[mm]);
118  }
119  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Split sideset %d\n",
120  mit->getMeshsetId());
121 
122  EntityHandle cubit_meshset = mit->getMeshset();
123  {
124  // get tet entities form back bit_level
125  EntityHandle ref_level_meshset = 0;
126  CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
128  bit_levels.back(), BitRefLevel().set(), MBTET, ref_level_meshset);
130  bit_levels.back(), BitRefLevel().set(), MBPRISM, ref_level_meshset);
131  Range ref_level_tets;
132  CHKERR moab.get_entities_by_handle(ref_level_meshset, ref_level_tets,
133  true);
134 
135  // get faces and test to split
136  CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(), true,
137  0);
138  // set new bit level
139  bit_levels.push_back(BitRefLevel().set(mm + 1));
140  // split faces and
141  CHKERR interface_ptr->splitSides(ref_level_meshset, bit_levels.back(),
142  cubit_meshset, false, true, 0);
143 
144  // clean meshsets
145  CHKERR moab.delete_entities(&ref_level_meshset, 1);
146  }
147  // Update cubit meshsets
149  (core.getInterface<MeshsetsManager &, 0>()),ciit) ) {
150 
151  EntityHandle cubit_meshset = ciit->meshset;
153  cubit_meshset, bit_levels.back(), cubit_meshset,
154  MBVERTEX, true);
156  cubit_meshset, bit_levels.back(), cubit_meshset, MBEDGE,
157  true);
159  cubit_meshset, bit_levels.back(), cubit_meshset, MBTRI,
160  true);
162  cubit_meshset, bit_levels.back(), cubit_meshset, MBTET,
163  true);
164  }
165  }
166 
167  if (squash_bit_levels == PETSC_TRUE) {
168  for (unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
169  CHKERR m_field.delete_ents_by_bit_ref(bit_levels[ll], bit_levels[ll],
170  true);
171  }
172  CHKERR bit_mng->shiftRightBitRef(bit_levels.size() - 1);
173  }
174 
175  if (output_vtk) {
176  EntityHandle meshset;
177  CHKERR moab.create_meshset(MESHSET_SET, meshset);
178  BitRefLevel bit;
179  if (squash_bit_levels)
180  bit = bit_levels[0];
181  else
182  bit = bit_levels.back();
184  bit, BitRefLevel().set(), MBTET, meshset);
185  CHKERR moab.write_file("out.vtk", "VTK", "", &meshset, 1);
186  CHKERR moab.delete_entities(&meshset, 1);
187  }
188 
189  CHKERR moab.write_file("out.h5m");
190  }
191  CATCH_ERRORS;
192 
194 
195  return 0;
196 }
Deprecated interface functions.
MoFEMErrorCode getEntitiesByTypeAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const EntityType type, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
PetscBool output_vtk
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:550
virtual MoFEMErrorCode delete_ents_by_bit_ref(const BitRefLevel bit, const BitRefLevel mask, const bool remove_parent=false, int verb=DEFAULT_VERBOSITY)=0
delete entities form mofem and moab database
Interface for managing meshsets containing materials and boundary conditions.
Core (interface) class.
Definition: Core.hpp:50
static char help[]
CubitMeshSet_multiIndex & getMeshsetsMultindex()
char mesh_file_name[255]
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
multi_index_container< CubitMeshSets, indexed_by< hashed_unique< tag< Meshset_mi_tag >, member< CubitMeshSets, EntityHandle, &CubitMeshSets::meshset > >, ordered_non_unique< tag< CubitMeshSets_mi_tag >, const_mem_fun< CubitMeshSets, unsigned long int, &CubitMeshSets::getBcTypeULong > >, ordered_non_unique< tag< CubitMeshSets_mask_meshset_mi_tag >, const_mem_fun< CubitMeshSets, unsigned long int, &CubitMeshSets::getMaksedBcTypeULong > >, ordered_non_unique< tag< CubitMeshSets_name >, const_mem_fun< CubitMeshSets, std::string, &CubitMeshSets::getName > >, hashed_unique< tag< Composite_Cubit_msId_And_MeshSetType_mi_tag >, composite_key< CubitMeshSets, const_mem_fun< CubitMeshSets, int, &CubitMeshSets::getMeshsetId >, const_mem_fun< CubitMeshSets, unsigned long int, &CubitMeshSets::getMaksedBcTypeULong > > > > > CubitMeshSet_multiIndex
Stores data about meshsets (see CubitMeshSets) storing data about boundary conditions,...
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode updateMeshsetByEntitiesChildren(const EntityHandle parent, const BitRefLevel &parent_bit, const BitRefLevel &parent_mask, const BitRefLevel &child_bit, const BitRefLevel &child_mask, const EntityHandle child, EntityType child_type, const bool recursive=false, int verb=0)
Get child entities form meshset containing parent entitiesSearch for refined entities of given type w...
Managing BitRefLevels.
Create interface from given surface and insert flat prisms in-between.
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
#define CHKERR
Inline error check.
Definition: definitions.h:602
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:291
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1788
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
MoFEMErrorCode shiftRightBitRef(const int shift, const BitRefLevel mask=BitRefLevel().set(), int verb=-1) const
right shift bit ref level
MoFEMErrorCode setBitRefLevelByDim(const EntityHandle meshset, const int dim, const BitRefLevel bit, int verb=QUIET) const
Set the Bit Ref Level By Dim object.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:77

Variable Documentation

◆ help

char help[] = "...\n\n"
static
Examples
split_sideset.cpp.

Definition at line 24 of file split_sideset.cpp.