v0.15.0
Loading...
Searching...
No Matches
Functions | Variables
split_sideset.cpp File Reference

Split sidesets/blockset and add internal surface. 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"
 
constexpr bool debug = false
 

Detailed Description

Split sidesets/blockset and add internal surface.

Definition in file split_sideset.cpp.

Function Documentation

◆ main()

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

Definition at line 13 of file split_sideset.cpp.

13 {
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 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_NULLPTR);
51 CHKERR PetscOptionsBool("-split_corner_edges", "if true output vtk file",
52 "", split_corner_edges, &split_corner_edges,
53 PETSC_NULLPTR);
54 CHKERR PetscOptionsBool("-output_vtk", "if true output vtk file", "",
55 output_vtk, &output_vtk, PETSC_NULLPTR);
56 PetscOptionsEnd();
57
58 moab::Core mb_instance;
59 moab::Interface &moab = mb_instance;
60
61 const char *option;
62 option = "";
63 CHKERR moab.load_file(mesh_file_name, 0, option);
64
65 // Create MoFEM database
66 MoFEM::Core core(moab);
67 MoFEM::Interface &m_field = core;
68
69 if (flg_file != PETSC_TRUE) {
70 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
71 "*** ERROR -my_file (mesh file needed)");
72 }
73 if (flg_list_of_sidesets != PETSC_TRUE && flg_block != PETSC_TRUE) {
74 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
75 "Block name or kist of sidesets not given ...");
76 }
77
78 // Get interface to meshsets manager
79 auto m_mng = m_field.getInterface<MeshsetsManager>();
80 // Get interface for splitting manager
81 auto interface_ptr = m_field.getInterface<PrismInterface>();
82 // Managing bits
83 auto bit_mng = m_field.getInterface<BitRefManager>();
84 // Refine mesh
85 auto refine_mng = m_field.getInterface<MeshRefinement>();
86
87 for (auto m : m_mng->getCubitMeshsetPtr(SIDESET)) {
88 MOFEM_LOG("SPLIT", Sev::verbose)
89 << "Sideset on the mesh id = " << m->getMeshsetId();
90 }
91
92 CHKERR bit_mng->setBitRefLevelByDim(0, 3, BitRefLevel().set(0));
93 std::vector<BitRefLevel> bit_levels;
94 bit_levels.push_back(BitRefLevel().set(0));
95
96 auto update_meshsets = [&]() {
97 MoFEMFunctionBegin;
98 for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
99 EntityHandle cubit_meshset = ciit->meshset;
100 CHKERR bit_mng->updateMeshsetByEntitiesChildren(
101 cubit_meshset, bit_levels.back(), cubit_meshset, MBMAXTYPE, true);
102 }
103 MoFEMFunctionReturn(0);
104 };
105
106 // refine corner mesh
107 if (split_corner_edges) {
108 Skinner skin(&m_field.get_moab());
109 auto meshset_of_edges_to_refine_ptr = get_temp_meshset_ptr(moab);
110
111 auto refine = [&](auto mit, auto meshset_of_edges_to_refine_ptr) {
112 MoFEMFunctionBegin;
113 Range faces;
114 CHKERR moab.get_entities_by_type(mit->getMeshset(), MBTRI, faces, true);
115 Range faces_edges;
116 CHKERR moab.get_adjacencies(faces, 1, true, faces_edges,
117 moab::Interface::UNION);
118
119 Range skin_edges;
120 CHKERR skin.find_skin(0, faces, false, skin_edges);
121 Range skin_verts;
122 CHKERR moab.get_connectivity(skin_edges, skin_verts, true);
123 Range vertex_edges;
124 CHKERR moab.get_adjacencies(skin_verts, 1, true, vertex_edges,
125 moab::Interface::UNION);
126 Range vertex_edges_verts;
127 CHKERR moab.get_connectivity(vertex_edges, vertex_edges_verts, true);
128 vertex_edges_verts = subtract(vertex_edges_verts, skin_verts);
129 Range vertex_edges_verts_edges;
130 CHKERR moab.get_adjacencies(vertex_edges_verts, 1, true,
131 vertex_edges_verts_edges,
132 moab::Interface::UNION);
133 vertex_edges = subtract(vertex_edges, vertex_edges_verts_edges);
134 vertex_edges = subtract(vertex_edges, skin_edges);
135 vertex_edges = intersect(vertex_edges, faces_edges);
136 CHKERR moab.add_entities(*meshset_of_edges_to_refine_ptr, vertex_edges);
137 MoFEMFunctionReturn(0);
138 };
139
140 for (int mm = 0; mm != nb_sidesets; mm++) {
141 CHKERR refine(m_mng->getCubitMeshsetPtr(sidesets[mm], SIDESET),
142 meshset_of_edges_to_refine_ptr);
143 }
144
145 for (auto m : m_mng->getCubitMeshsetPtr(std::regex(
146
147 (boost::format("%s") % block_name).str()
148
149 ))
150
151 ) {
152 CHKERR refine(m, meshset_of_edges_to_refine_ptr);
153 }
154
155 int nb_tris;
156 CHKERR moab.get_number_entities_by_type(*meshset_of_edges_to_refine_ptr,
157 MBEDGE, nb_tris, true);
158 MOFEM_LOG("SPLIT", Sev::inform) << "Refine corner edges " << nb_tris;
159
160 if (debug)
161 CHKERR moab.write_file("out_edges_to_refine.vtk", "VTK", "",
162 meshset_of_edges_to_refine_ptr->get_ptr(), 1);
163
164 bit_levels.push_back(BitRefLevel().set(1));
165 CHKERR refine_mng->addVerticesInTheMiddleOfEdges(
166 *meshset_of_edges_to_refine_ptr, bit_levels.back());
167 CHKERR refine_mng->refineTets(0, bit_levels.back());
168 CHKERR update_meshsets();
169 }
170
171 auto split = [&](auto mit) {
172 MoFEMFunctionBegin;
173 const auto cubit_meshset = mit->getMeshset();
174 {
175 // get tet entities form back bit_level
176 auto ref_level_meshset_ptr = get_temp_meshset_ptr(moab);
177 CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit_levels.back(),
178 BitRefLevel().set(), MBTET,
179 *ref_level_meshset_ptr);
180 CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
181 bit_levels.back(), BitRefLevel().set(), MBPRISM,
182 *ref_level_meshset_ptr);
183 Range ref_level_tets;
184 CHKERR moab.get_entities_by_handle(*ref_level_meshset_ptr,
185 ref_level_tets, true);
186
187 // get faces and test to split
188 CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(), true,
189 0);
190 // set new bit level
191 MOFEM_LOG("SPLIT", Sev::verbose)
192 << "Add bit level " << bit_levels.size();
193 bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
194 // split faces and
195 CHKERR interface_ptr->splitSides(*ref_level_meshset_ptr,
196 bit_levels.back(), cubit_meshset,
197 add_prisms, true, 0);
198 }
199 // Update cubit meshsets
200 CHKERR update_meshsets();
201 MoFEMFunctionReturn(0);
202 };
203
204 for (int mm = 0; mm != nb_sidesets; mm++) {
205 CHKERR split(m_mng->getCubitMeshsetPtr(sidesets[mm], SIDESET));
206 }
207
208 for (auto m : m_mng->getCubitMeshsetPtr(std::regex(
209
210 (boost::format("%s") % block_name).str()
211
212 ))
213
214 ) {
215 CHKERR split(m);
216 }
217
218 if (squash_bit_levels == PETSC_TRUE) {
219 for (unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
220 CHKERR m_field.delete_ents_by_bit_ref(bit_levels[ll], bit_levels[ll],
221 true);
222 }
223 CHKERR bit_mng->shiftRightBitRef(bit_levels.size() - 1);
224 }
225
226 // extract filter and delete extra entities not in meshsets
227 Range meshset_ents;
228 auto &list = m_mng->getMeshsetsMultindex();
229
230 for (auto &m : list) {
231 Range ents;
232 CHKERR moab.get_entities_by_handle(m.getMeshset(), ents, false);
233 meshset_ents.merge(ents);
234 }
235
236 Range not_in_meshset_ents;
237 for (EntityType type = MBEDGE; type != MBENTITYSET; type++) {
238 Range ents;
239 CHKERR moab.get_entities_by_type(0, type, ents, false);
240 not_in_meshset_ents.merge(ents);
241 }
242
243 not_in_meshset_ents = subtract(not_in_meshset_ents, meshset_ents);
244
245 CHKERR moab.delete_entities(not_in_meshset_ents);
246
247 if (output_vtk) {
248 auto meshset_ptr = get_temp_meshset_ptr(moab);
249 BitRefLevel bit;
250 if (squash_bit_levels)
251 bit = bit_levels[0];
252 else
253 bit = bit_levels.back();
254 CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
255 MBTET, *meshset_ptr);
256 CHKERR moab.write_file("out.vtk", "VTK", "", meshset_ptr->get_ptr(), 1);
257 }
258
259 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
260 if (pcomm == NULL)
261 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
262 "Communicator should be allocated");
263
264 CHKERR pcomm->assign_global_ids(0, 3, 1, false);
265 CHKERR moab.write_file("out.h5m");
266 }
267 CATCH_ERRORS;
268
269 MoFEM::Core::Finalize();
270
271 return 0;
272}
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
SIDESET
@ SIDESET
Definition definitions.h:160
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition definitions.h:36
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition LogManager.cpp:392
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition LogManager.hpp:347
_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:57
bit
auto bit
set bit
Definition hanging_node_approx.cpp:75
output_vtk
PetscBool output_vtk
Definition mesh_smoothing.cpp:25
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
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::debug
static const bool debug
Definition ConstrainMatrixCtx.cpp:9
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition Templates.hpp:1980
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
help
static char help[]
Definition split_sideset.cpp:10
MoFEM::BitRefManager
Managing BitRefLevels.
Definition BitRefManager.hpp:21
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
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 from MoFEM and MOAB database.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
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::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
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:293
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition LogManager.cpp:339
MoFEM::MeshRefinement
Mesh refinement interface.
Definition MeshRefinement.hpp:26
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition MeshsetsManager.hpp:208
MoFEM::PrismInterface
Create interface from given surface and insert flat prisms in-between.
Definition PrismInterface.hpp:23
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93

Variable Documentation

◆ debug

constexpr bool debug = false
constexpr

Definition at line 11 of file split_sideset.cpp.

◆ help

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

Definition at line 10 of file split_sideset.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow