13 {
15
16 auto core_log = logging::core::get();
17 core_log->add_sink(
21
22 try {
23
24
26 PetscBool flg_file = PETSC_FALSE;
27 PetscBool squash_bit_levels = 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",
41 CHKERR PetscOptionsString(
"-file_name",
"mesh file name",
"",
"mesh.h5m",
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",
"",
56 PetscOptionsEnd();
57
58 moab::Core mb_instance;
59 moab::Interface &moab = mb_instance;
60
61 const char *option;
62 option = "";
64
65
68
69 if (flg_file != PETSC_TRUE) {
71 "*** ERROR -my_file (mesh file needed)");
72 }
73 if (flg_list_of_sidesets != PETSC_TRUE && flg_block != PETSC_TRUE) {
75 "Block name or kist of sidesets not given ...");
76 }
77
78
80
82
84
86
87 for (
auto m : m_mng->getCubitMeshsetPtr(
SIDESET)) {
89 <<
"Sideset on the mesh id = " <<
m->getMeshsetId();
90 }
91
93 std::vector<BitRefLevel> bit_levels;
95
96 auto update_meshsets = [&]() {
100 CHKERR bit_mng->updateMeshsetByEntitiesChildren(
101 cubit_meshset, bit_levels.back(), cubit_meshset, MBMAXTYPE, true);
102 }
104 };
105
106
107 if (split_corner_edges) {
110
111 auto refine = [&](auto mit, auto meshset_of_edges_to_refine_ptr) {
114 CHKERR moab.get_entities_by_type(mit->getMeshset(), MBTRI, faces,
true);
116 CHKERR moab.get_adjacencies(faces, 1,
true, faces_edges,
117 moab::Interface::UNION);
118
120 CHKERR skin.find_skin(0, faces,
false, skin_edges);
122 CHKERR moab.get_connectivity(skin_edges, skin_verts,
true);
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);
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
161 CHKERR moab.write_file(
"out_edges_to_refine.vtk",
"VTK",
"",
162 meshset_of_edges_to_refine_ptr->get_ptr(), 1);
163
165 CHKERR refine_mng->addVerticesInTheMiddleOfEdges(
166 *meshset_of_edges_to_refine_ptr, bit_levels.back());
167 CHKERR refine_mng->refineTets(0, bit_levels.back());
169 }
170
171 auto split = [&](auto mit) {
173 const auto cubit_meshset = mit->getMeshset();
174 {
175
177 CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit_levels.back(),
179 *ref_level_meshset_ptr);
180 CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
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
188 CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(),
true,
189 0);
190
192 << "Add bit level " << bit_levels.size();
193 bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
194
195 CHKERR interface_ptr->splitSides(*ref_level_meshset_ptr,
196 bit_levels.back(), cubit_meshset,
197 add_prisms, true, 0);
198 }
199
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 ) {
216 }
217
218 if (squash_bit_levels == PETSC_TRUE) {
219 for (unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
221 true);
222 }
223 CHKERR bit_mng->shiftRightBitRef(bit_levels.size() - 1);
224 }
225
226
228 auto &list = m_mng->getMeshsetsMultindex();
229
230 for (
auto &
m : list) {
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++) {
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
250 if (squash_bit_levels)
252 else
253 bit = bit_levels.back();
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)
262 "Communicator should be allocated");
263
264 CHKERR pcomm->assign_global_ids(0, 3, 1,
false);
265 CHKERR moab.write_file(
"out.h5m");
266 }
268
270
271 return 0;
272}
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
FTensor::Index< 'm', 3 > m
virtual moab::Interface & get_moab()=0
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
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Mesh refinement interface.
Interface for managing meshsets containing materials and boundary conditions.
Create interface from given surface and insert flat prisms in-between.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.