v0.14.0
split_sideset.cpp
Go to the documentation of this file.
1 /** \file split_sideset.cpp
2  \brief Split sidesets/blockset and add internal surface
3  \example split_sideset.cpp
4 
5 */
6 
7 #include <MoFEM.hpp>
8 using namespace MoFEM;
9 
10 static char help[] = "...\n\n";
11 constexpr bool debug = false;
12 
13 int main(int argc, char *argv[]) {
14  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
15 
16  auto core_log = logging::core::get();
17  core_log->add_sink(
18  LogManager::createSink(LogManager::getStrmWorld(), "SPLIT"));
19  LogManager::setLog("SPLIT");
20  MOFEM_LOG_TAG("SPLIT", "split");
21 
22  try {
23 
24  // global variables
25  char mesh_file_name[255] = "mesh.h5m";
26  PetscBool flg_file = PETSC_FALSE;
27  PetscBool squash_bit_levels = PETSC_TRUE;
28  PetscBool output_vtk = PETSC_TRUE;
29  PetscBool add_prisms = PETSC_FALSE;
30  PetscBool split_corner_edges = PETSC_FALSE;
31  char block_name[255] = "CRACK";
32  PetscBool flg_block = PETSC_FALSE;
33  int nb_sidesets = 10;
34  int sidesets[nb_sidesets];
35  PetscBool flg_list_of_sidesets = PETSC_FALSE;
36 
37  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Split sides options",
38  "none");
39  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
40  mesh_file_name, 255, &flg_file);
41  CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
42  mesh_file_name, 255, &flg_file);
43  CHKERR PetscOptionsBool("-squash_bit_levels", "squash bit levels", "",
44  squash_bit_levels, &squash_bit_levels, NULL);
45  CHKERR PetscOptionsIntArray("-side_sets", "get list of sidesets", "",
46  sidesets, &nb_sidesets, &flg_list_of_sidesets);
47  CHKERR PetscOptionsString("-block_name", "block_name", "", "mesh.h5m",
48  block_name, 255, &flg_block);
49  CHKERR PetscOptionsBool("-add_prisms", "add prisms", "", add_prisms,
50  &add_prisms, PETSC_NULL);
51  CHKERR PetscOptionsBool("-split_corner_edges", "if true output vtk file",
52  "", split_corner_edges, &split_corner_edges,
53  PETSC_NULL);
54  CHKERR PetscOptionsBool("-output_vtk", "if true output vtk file", "",
55  output_vtk, &output_vtk, PETSC_NULL);
56  ierr = PetscOptionsEnd();
57  CHKERRG(ierr);
58 
59  moab::Core mb_instance;
60  moab::Interface &moab = mb_instance;
61 
62  const char *option;
63  option = "";
64  CHKERR moab.load_file(mesh_file_name, 0, option);
65 
66  // Create MoFEM database
67  MoFEM::Core core(moab);
68  MoFEM::Interface &m_field = core;
69 
70  if (flg_file != PETSC_TRUE) {
71  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
72  "*** ERROR -my_file (mesh file needed)");
73  }
74  if (flg_list_of_sidesets != PETSC_TRUE && flg_block != PETSC_TRUE) {
75  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
76  "Block name or kist of sidesets not given ...");
77  }
78 
79  // Get interface to meshsets manager
80  auto m_mng = m_field.getInterface<MeshsetsManager>();
81  // Get interface for splitting manager
82  auto interface_ptr = m_field.getInterface<PrismInterface>();
83  // Managing bits
84  auto bit_mng = m_field.getInterface<BitRefManager>();
85  // Refine mesh
86  auto refine_mng = m_field.getInterface<MeshRefinement>();
87 
88  for (auto m : m_mng->getCubitMeshsetPtr(SIDESET)) {
89  MOFEM_LOG("SPLIT", Sev::verbose)
90  << "Sideset on the mesh id = " << m->getMeshsetId();
91  }
92 
93  CHKERR bit_mng->setBitRefLevelByDim(0, 3, BitRefLevel().set(0));
94  std::vector<BitRefLevel> bit_levels;
95  bit_levels.push_back(BitRefLevel().set(0));
96 
97  auto update_meshsets = [&]() {
98  MoFEMFunctionBegin;
99  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
100  EntityHandle cubit_meshset = ciit->meshset;
101  CHKERR bit_mng->updateMeshsetByEntitiesChildren(
102  cubit_meshset, bit_levels.back(), cubit_meshset, MBMAXTYPE, true);
103  }
104  MoFEMFunctionReturn(0);
105  };
106 
107  // refine corner mesh
108  if (split_corner_edges) {
109  Skinner skin(&m_field.get_moab());
110  auto meshset_of_edges_to_refine_ptr = get_temp_meshset_ptr(moab);
111 
112  auto refine = [&](auto mit, auto meshset_of_edges_to_refine_ptr) {
113  MoFEMFunctionBegin;
114  Range faces;
115  CHKERR moab.get_entities_by_type(mit->getMeshset(), MBTRI, faces, true);
116  Range faces_edges;
117  CHKERR moab.get_adjacencies(faces, 1, true, faces_edges,
118  moab::Interface::UNION);
119 
120  Range skin_edges;
121  CHKERR skin.find_skin(0, faces, false, skin_edges);
122  Range skin_verts;
123  CHKERR moab.get_connectivity(skin_edges, skin_verts, true);
124  Range vertex_edges;
125  CHKERR moab.get_adjacencies(skin_verts, 1, true, vertex_edges,
126  moab::Interface::UNION);
127  Range vertex_edges_verts;
128  CHKERR moab.get_connectivity(vertex_edges, vertex_edges_verts, true);
129  vertex_edges_verts = subtract(vertex_edges_verts, skin_verts);
130  Range vertex_edges_verts_edges;
131  CHKERR moab.get_adjacencies(vertex_edges_verts, 1, true,
132  vertex_edges_verts_edges,
133  moab::Interface::UNION);
134  vertex_edges = subtract(vertex_edges, vertex_edges_verts_edges);
135  vertex_edges = subtract(vertex_edges, skin_edges);
136  vertex_edges = intersect(vertex_edges, faces_edges);
137  CHKERR moab.add_entities(*meshset_of_edges_to_refine_ptr, vertex_edges);
138  MoFEMFunctionReturn(0);
139  };
140 
141  for (int mm = 0; mm != nb_sidesets; mm++) {
142  CHKERR refine(m_mng->getCubitMeshsetPtr(sidesets[mm], SIDESET),
143  meshset_of_edges_to_refine_ptr);
144  }
145 
146  for (auto m : m_mng->getCubitMeshsetPtr(std::regex(
147 
148  (boost::format("%s(.*)") % block_name).str()
149 
150  ))
151 
152  ) {
153  CHKERR refine(m, meshset_of_edges_to_refine_ptr);
154  }
155 
156  int nb_tris;
157  CHKERR moab.get_number_entities_by_type(*meshset_of_edges_to_refine_ptr,
158  MBEDGE, nb_tris, true);
159  MOFEM_LOG("SPLIT", Sev::inform) << "Refine corner edges " << nb_tris;
160 
161  if (debug)
162  CHKERR moab.write_file("out_edges_to_refine.vtk", "VTK", "",
163  meshset_of_edges_to_refine_ptr->get_ptr(), 1);
164 
165  bit_levels.push_back(BitRefLevel().set(1));
166  CHKERR refine_mng->addVerticesInTheMiddleOfEdges(
167  *meshset_of_edges_to_refine_ptr, bit_levels.back());
168  CHKERR refine_mng->refineTets(0, bit_levels.back());
169  CHKERR update_meshsets();
170  }
171 
172  auto split = [&](auto mit) {
173  MoFEMFunctionBegin;
174  const auto cubit_meshset = mit->getMeshset();
175  {
176  // get tet entities form back bit_level
177  auto ref_level_meshset_ptr = get_temp_meshset_ptr(moab);
178  CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit_levels.back(),
179  BitRefLevel().set(), MBTET,
180  *ref_level_meshset_ptr);
181  CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
182  bit_levels.back(), BitRefLevel().set(), MBPRISM,
183  *ref_level_meshset_ptr);
184  Range ref_level_tets;
185  CHKERR moab.get_entities_by_handle(*ref_level_meshset_ptr,
186  ref_level_tets, true);
187 
188  // get faces and test to split
189  CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(), true,
190  0);
191  // set new bit level
192  MOFEM_LOG("SPLIT", Sev::verbose)
193  << "Add bit level " << bit_levels.size();
194  bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
195  // split faces and
196  CHKERR interface_ptr->splitSides(*ref_level_meshset_ptr,
197  bit_levels.back(), cubit_meshset,
198  add_prisms, true, 0);
199  }
200  // Update cubit meshsets
201  CHKERR update_meshsets();
202  MoFEMFunctionReturn(0);
203  };
204 
205  for (int mm = 0; mm != nb_sidesets; mm++) {
206  CHKERR split(m_mng->getCubitMeshsetPtr(sidesets[mm], SIDESET));
207  }
208 
209  for (auto m : m_mng->getCubitMeshsetPtr(std::regex(
210 
211  (boost::format("%s(.*)") % block_name).str()
212 
213  ))
214 
215  ) {
216  CHKERR split(m);
217  }
218 
219  if (squash_bit_levels == PETSC_TRUE) {
220  for (unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
221  CHKERR m_field.delete_ents_by_bit_ref(bit_levels[ll], bit_levels[ll],
222  true);
223  }
224  CHKERR bit_mng->shiftRightBitRef(bit_levels.size() - 1);
225  }
226 
227  if (output_vtk) {
228  auto meshset_ptr = get_temp_meshset_ptr(moab);
229  BitRefLevel bit;
230  if (squash_bit_levels)
231  bit = bit_levels[0];
232  else
233  bit = bit_levels.back();
234  CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
235  MBTET, *meshset_ptr);
236  CHKERR moab.write_file("out.vtk", "VTK", "", meshset_ptr->get_ptr(), 1);
237  }
238 
239  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
240  if (pcomm == NULL)
241  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
242  "Communicator should be allocated");
243 
244  CHKERR pcomm->assign_global_ids(0, 3, 1, false);
245  CHKERR moab.write_file("out.h5m");
246  }
247  CATCH_ERRORS;
248 
249  MoFEM::Core::Finalize();
250 
251  return 0;
252 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
SIDESET
@ SIDESET
Definition: definitions.h:160
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::debug
const static bool debug
Definition: ConstrainMatrixCtx.cpp:9
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
EntityHandle
MoFEM::PrismInterface
Create interface from given surface and insert flat prisms in-between.
Definition: PrismInterface.hpp:23
main
int main(int argc, char *argv[])
Definition: split_sideset.cpp:13
MoFEM.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
output_vtk
PetscBool output_vtk
Definition: mesh_smoothing.cpp:25
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1886
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
help
static char help[]
Definition: split_sideset.cpp:10
MoFEM::CoreInterface::delete_ents_by_bit_ref
virtual MoFEMErrorCode delete_ents_by_bit_ref(const BitRefLevel bit, const BitRefLevel mask, const bool remove_parent=false, int verb=DEFAULT_VERBOSITY, MoFEMTypes mf=MF_ZERO)=0
delete entities form mofem and moab database
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
Range
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
_IT_CUBITMESHSETS_FOR_LOOP_
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
Definition: MeshsetsManager.hpp:34
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
Generated by Doxygen 1.8.17 and hosted at University of Glasgow