v0.14.0
split_sideset.cpp
Go to the documentation of this file.
1 /** \file split_sideset.cpp
2  \brief Split sidesets
3  \example split_sideset.cpp
4 
5 */
6 
7 
8 
9 #include <MoFEM.hpp>
10 using namespace MoFEM;
11 
12 static char help[] = "...\n\n";
13 constexpr bool debug = false;
14 
15 int main(int argc, char *argv[]) {
16  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
17 
18  auto core_log = logging::core::get();
19  core_log->add_sink(
20  LogManager::createSink(LogManager::getStrmWorld(), "SPLIT"));
21  LogManager::setLog("SPLIT");
22  MOFEM_LOG_TAG("SPLIT", "split");
23 
24  try {
25 
26  // global variables
27  char mesh_file_name[255];
28  PetscBool flg_file = PETSC_FALSE;
29  PetscBool squash_bit_levels = PETSC_TRUE;
30  PetscBool flg_list_of_sidesets = PETSC_FALSE;
31  PetscBool output_vtk = PETSC_TRUE;
32  PetscBool add_prisms = PETSC_FALSE;
33  PetscBool split_corner_edges = PETSC_FALSE;
34  int nb_sidesets = 10;
35  int sidesets[nb_sidesets];
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 PetscOptionsBool("-squash_bit_levels", "squash bit levels", "",
42  squash_bit_levels, &squash_bit_levels, NULL);
43  CHKERR PetscOptionsIntArray("-side_sets", "get list of sidesets", "",
44  sidesets, &nb_sidesets, &flg_list_of_sidesets);
45  CHKERR PetscOptionsBool("-output_vtk", "if true outout vtk file", "",
46  output_vtk, &output_vtk, PETSC_NULL);
47  CHKERR PetscOptionsBool("-add_prisms", "if true outout vtk file", "",
48  add_prisms, &add_prisms, PETSC_NULL);
49  CHKERR PetscOptionsBool("-split_corner_edges", "if true outout vtk file",
50  "", split_corner_edges, &split_corner_edges,
51  PETSC_NULL);
52  CHKERR PetscOptionsBool("-output_vtk", "if true outout vtk file", "",
53  output_vtk, &output_vtk, PETSC_NULL);
54  ierr = PetscOptionsEnd();
55  CHKERRG(ierr);
56 
57  moab::Core mb_instance;
58  moab::Interface &moab = mb_instance;
59 
60  const char *option;
61  option = "";
62  CHKERR moab.load_file(mesh_file_name, 0, option);
63 
64  // Create MoFEM database
65  MoFEM::Core core(moab);
66  MoFEM::Interface &m_field = core;
67 
68  if (flg_file != PETSC_TRUE) {
69  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
70  "*** ERROR -my_file (mesh file needed)");
71  }
72  if (flg_list_of_sidesets != PETSC_TRUE) {
73  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
74  "List of sidesets not given -my_side_sets ...");
75  }
76 
77  // Get interface to meshsets manager
78  auto m_mng = m_field.getInterface<MeshsetsManager>();
79  // Get interface for splitting manager
80  auto interface_ptr = m_field.getInterface<PrismInterface>();
81  // Managing bits
82  auto bit_mng = m_field.getInterface<BitRefManager>();
83  // Refine mesh
84  auto refine_mng = m_field.getInterface<MeshRefinement>();
85 
86  auto &meshsets_index = m_mng->getMeshsetsMultindex();
87  auto &m_by_type = meshsets_index.get<CubitMeshsetMaskedType_mi_tag>();
88  auto &m_by_id_and_type =
89  meshsets_index.get<Composite_Cubit_msId_And_MeshsetType_mi_tag>();
90 
91  for (auto mit = m_by_type.lower_bound(SIDESET);
92  mit != m_by_type.upper_bound(SIDESET); mit++) {
93  MOFEM_LOG("SPLIT", Sev::verbose)
94  << "Sideset on the mesh id = " << mit->getMeshsetId();
95  }
96 
97  CHKERR bit_mng->setBitRefLevelByDim(0, 3, BitRefLevel().set(0));
98  std::vector<BitRefLevel> bit_levels;
99  bit_levels.push_back(BitRefLevel().set(0));
100 
101  auto update_meshsets = [&]() {
102  MoFEMFunctionBegin;
103  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
104  EntityHandle cubit_meshset = ciit->meshset;
105  CHKERR bit_mng->updateMeshsetByEntitiesChildren(
106  cubit_meshset, bit_levels.back(), cubit_meshset, MBMAXTYPE, true);
107  }
108  MoFEMFunctionReturn(0);
109  };
110 
111  // refine corner mesh
112  if (split_corner_edges) {
113  Skinner skin(&m_field.get_moab());
114  auto meshset_of_edges_to_refine_ptr = get_temp_meshset_ptr(moab);
115 
116  for (int mm = 0; mm != nb_sidesets; mm++) {
117 
118  // find side set
119  auto mit =
120  m_by_id_and_type.find(boost::make_tuple(sidesets[mm], SIDESET));
121  if (mit == m_by_id_and_type.end())
122  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
123  "No sideset in database id = %d", sidesets[mm]);
124 
125  Range faces;
126  CHKERR moab.get_entities_by_type(mit->getMeshset(), MBTRI, faces, true);
127  Range faces_edges;
128  CHKERR moab.get_adjacencies(faces, 1, true, faces_edges,
129  moab::Interface::UNION);
130 
131  Range skin_edges;
132  CHKERR skin.find_skin(0, faces, false, skin_edges);
133  Range skin_verts;
134  CHKERR moab.get_connectivity(skin_edges, skin_verts, true);
135  Range vertex_edges;
136  CHKERR moab.get_adjacencies(skin_verts, 1, true, vertex_edges,
137  moab::Interface::UNION);
138  Range vertex_edges_verts;
139  CHKERR moab.get_connectivity(vertex_edges, vertex_edges_verts, true);
140  vertex_edges_verts = subtract(vertex_edges_verts, skin_verts);
141  Range vertex_edges_verts_edges;
142  CHKERR moab.get_adjacencies(vertex_edges_verts, 1, true,
143  vertex_edges_verts_edges,
144  moab::Interface::UNION);
145  vertex_edges = subtract(vertex_edges, vertex_edges_verts_edges);
146  vertex_edges = subtract(vertex_edges, skin_edges);
147  vertex_edges = intersect(vertex_edges, faces_edges);
148  CHKERR moab.add_entities(*meshset_of_edges_to_refine_ptr, vertex_edges);
149  }
150 
151  int nb_tris;
152  CHKERR moab.get_number_entities_by_type(*meshset_of_edges_to_refine_ptr,
153  MBEDGE, nb_tris, true);
154  MOFEM_LOG("SPLIT", Sev::inform) << "Refine corner edges " << nb_tris;
155 
156  if (debug)
157  CHKERR moab.write_file("out_edges_to_refine.vtk", "VTK", "",
158  meshset_of_edges_to_refine_ptr->get_ptr(), 1);
159 
160  bit_levels.push_back(BitRefLevel().set(1));
161  CHKERR refine_mng->addVerticesInTheMiddleOfEdges(
162  *meshset_of_edges_to_refine_ptr, bit_levels.back());
163  CHKERR refine_mng->refineTets(0, bit_levels.back());
164  CHKERR update_meshsets();
165  }
166 
167  // iterate sideset and split
168  for (int mm = 0; mm != nb_sidesets; mm++) {
169 
170  // find side set
171  auto mit =
172  m_by_id_and_type.find(boost::make_tuple(sidesets[mm], SIDESET));
173  if (mit == m_by_id_and_type.end())
174  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
175  "No sideset in database id = %d", sidesets[mm]);
176 
177  MOFEM_LOG("SPLIT", Sev::inform)
178  << "Split sideset " << mit->getMeshsetId();
179 
180  const auto cubit_meshset = mit->getMeshset();
181  {
182  // get tet entities form back bit_level
183  auto ref_level_meshset_ptr = get_temp_meshset_ptr(moab);
184  CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit_levels.back(),
185  BitRefLevel().set(), MBTET,
186  *ref_level_meshset_ptr);
187  CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
188  bit_levels.back(), BitRefLevel().set(), MBPRISM,
189  *ref_level_meshset_ptr);
190  Range ref_level_tets;
191  CHKERR moab.get_entities_by_handle(*ref_level_meshset_ptr,
192  ref_level_tets, true);
193 
194  // get faces and test to split
195  CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(), true,
196  0);
197  // set new bit level
198  MOFEM_LOG("SPLIT", Sev::verbose)
199  << "Add bit level " << bit_levels.size();
200  bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
201  // split faces and
202  CHKERR interface_ptr->splitSides(*ref_level_meshset_ptr,
203  bit_levels.back(), cubit_meshset,
204  add_prisms, true, 0);
205  }
206  // Update cubit meshsets
207  CHKERR update_meshsets();
208  }
209 
210  if (squash_bit_levels == PETSC_TRUE) {
211  for (unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
212  CHKERR m_field.delete_ents_by_bit_ref(bit_levels[ll], bit_levels[ll],
213  true);
214  }
215  CHKERR bit_mng->shiftRightBitRef(bit_levels.size() - 1);
216  }
217 
218  if (output_vtk) {
219  auto meshset_ptr = get_temp_meshset_ptr(moab);
220  BitRefLevel bit;
221  if (squash_bit_levels)
222  bit = bit_levels[0];
223  else
224  bit = bit_levels.back();
225  CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
226  MBTET, *meshset_ptr);
227  CHKERR moab.write_file("out.vtk", "VTK", "", meshset_ptr->get_ptr(), 1);
228  }
229 
230  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
231  if (pcomm == NULL)
232  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
233  "Communicator should be allocated");
234 
235  CHKERR pcomm->assign_global_ids(0, 3, 1, false);
236  CHKERR moab.write_file("out.h5m");
237  }
238  CATCH_ERRORS;
239 
240  MoFEM::Core::Finalize();
241 
242  return 0;
243 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
SIDESET
@ SIDESET
Definition: definitions.h:147
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
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:15
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:1975
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
output_vtk
PetscBool output_vtk
Definition: mesh_smoothing.cpp:26
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:1857
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:12
debug
constexpr bool debug
Definition: split_sideset.cpp:13
MoFEM::Composite_Cubit_msId_And_MeshsetType_mi_tag
Definition: TagMultiIndices.hpp:15
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:23
MoFEM::CubitMeshsetMaskedType_mi_tag
Definition: TagMultiIndices.hpp:13
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:372
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
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
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:416
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
Generated by Doxygen 1.8.17 and hosted at University of Glasgow