v0.9.0
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<
89  CubitMeshSets_mask_meshset_mi_tag>::type CMeshsetByType;
90  typedef CMeshsetByType::iterator CMIteratorByType;
91  typedef CubitMeshSet_multiIndex::index<
92  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type CMeshsetByIdType;
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 }
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 tetrahedral in children sets and add prism elements
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
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.
Make interface on given faces and create flat prism in that space.
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
create two children meshsets in the meshset containing tetrahedral on two sides of faces
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:285
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:433
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:61

Variable Documentation

◆ help

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

Definition at line 24 of file split_sideset.cpp.