v0.15.0
Loading...
Searching...
No Matches
CutMeshInterface.hpp
Go to the documentation of this file.
1/** \file CutMeshInterface.hpp
2 * \brief Cut mesh interface
3 *
4 * \ingroup mesh_cut
5 */
6
7#ifndef __CUTMESHINTERFACE_HPP__
8 #define __CUTMESHINTERFACE_HPP__
9
10 #include "UnknownInterface.hpp"
11
12namespace MoFEM {
13
14/**
15 * \brief Interface to cut meshes
16 *
17 * \ingroup mesh_cut
18 */
20
21 MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
22 UnknownInterface **iface) const;
23
25 CutMeshInterface(const MoFEM::Core &core);
27
31
33
34 /**
35 * \brief Get options from command line
36 * @return error code
37 */
40 PetscOptionsBegin(PETSC_COMM_WORLD, "cut_", "MOFEM Cut mesh options",
41 "none");
42 CHKERR PetscOptionsInt("-linesearch_steps",
43 "number of bisection steps which line search do to "
44 "find optimal merged nodes position",
46 PETSC_NULLPTR);
47
48 CHKERR PetscOptionsInt(
49 "-max_merging_cycles", "number of maximal merging cycles", "",
50 nbMaxMergingCycles, &nbMaxMergingCycles, PETSC_NULLPTR);
51
52 CHKERR PetscOptionsScalar("-project_entities_quality_trashold",
53 "project entities quality trashold", "",
55 &projectEntitiesQualityTrashold, PETSC_NULLPTR);
56
57 PetscOptionsEnd();
59 }
60
62
63 /**
64 * \brief set surface entities
65 * @param surface entities which going to be added
66 * @return error code
67 */
69
70 /**
71 * \brief copy surface entities
72 * @param surface entities which going to be added
73 * @return error code
74 */
76 double *shift = NULL, double *origin = NULL,
77 double *transform = NULL,
78 const std::string save_mesh = "");
79
80 /**
81 * \brief set volume entities
82 * @param volume entities which going to be added
83 * @return error code
84 */
85 MoFEMErrorCode setVolume(const Range volume);
86
87 /**
88 * @brief Set the constrain surface object
89 *
90 * Add surfaces which are restricted by mesh cut. Example of surface which
91 * needs to be respected is an interface, the boundary between two materials,
92 * or crack surface.
93 *
94 * @param surf
95 * @return MoFEMErrorCode
96 */
98
99 /**
100 * \brief merge surface entities
101 * @param surface entities which going to be added
102 * @return error code
103 */
105
106 /**
107 * \brief merge volume entities
108 * @param volume entities which going to be added
109 * @return error code
110 */
111 MoFEMErrorCode mergeVolumes(const Range volume);
112
114 const double rel_tol,
115 const double abs_tol, Tag th = nullptr,
116 const bool debug = false);
117
118 MoFEMErrorCode snapSurfaceToEdges(const Range surface_edges,
119 const Range fixed_edges,
120 const double rel_tol, const double abs_tol,
121 Tag th = nullptr, const bool debug = false);
122
123 /**
124 * \brief build tree
125 * @return error code
126 */
128
129 /**
130 * @brief Cut mesh only
131 *
132 * @param vol
133 * @param cut_bit
134 * @param th
135 * @param tol_geometry
136 * @param tol_cut_close
137 * @param fixed_edges
138 * @param corner_nodes
139 * @param update_meshsets
140 * @param debug
141 * @return MoFEMErrorCode
142 */
143 MoFEMErrorCode cutOnly(Range vol, const BitRefLevel cut_bit, Tag th,
144 const double tol_geometry, const double tol_cut_close,
145 Range *fixed_edges = NULL, Range *corner_nodes = NULL,
146 const bool update_meshsets = false,
147 const bool debug = false);
148
149 /**
150 * @brief Trim mesh only
151 *
152 * @param trim_bit
153 * @param th
154 * @param tol_cut_close
155 * @param fixed_edges
156 * @param corner_nodes
157 * @param update_meshsets
158 * @param debug
159 * @return MoFEMErrorCode
160 */
161 MoFEMErrorCode trimOnly(const BitRefLevel trim_bit, Tag th,
162 const double tol_trim_close,
163 Range *fixed_edges = NULL, Range *corner_nodes = NULL,
164 const bool update_meshsets = false,
165 const bool debug = false);
166
167 /**
168 * @brief Cut and trim
169 *
170 * @param first_bit
171 * @param th
172 * @param tol_geometry
173 * @param tol_cut_close
174 * @param tol_trim_close
175 * @param fixed_edges
176 * @param corner_nodes
177 * @param update_meshsets
178 * @param debug
179 * @return MoFEMErrorCode
180 */
182 cutAndTrim(int &first_bit, Tag th, const double tol_geometry,
183 const double tol_cut_close, const double tol_trim_close,
184 Range *fixed_edges = NULL, Range *corner_nodes = NULL,
185 const bool update_meshsets = false, const bool debug = false);
186
187 /**
188 * @brief Cut, trim and merge
189 *
190 * @param first_bit first bit of bit revel, subsequent set bits are for trim
191 * and merge
192 * @param fraction_level fraction of edges merged at each merge step
193 * @param th tag storring mesh node positions
194 * @param tol_geometry tolerance how mesh node should be close to cut surface
195 * (mesh node is moved), should be small
196 * @param tol_cut_close how crack node should be close to mesh (cut surface
197 * node is moved), can be big
198 * @param tol_trim_close how front node should be close to mesh, can be big
199 * @param fixed_edges edges on which nodes can not be moved
200 * @param corner_nodes nodes which can not be moved
201 * @param update_meshsets update meshsets by parents
202 * @param debug switch on debugging
203 * @return MoFEMErrorCode
204 */
205 MoFEMErrorCode cutTrimAndMerge(int &first_bit, const int fraction_level,
206 Tag th, const double tol_geometry,
207 const double tol_cut_close,
208 const double tol_trim_close,
209 Range &fixed_edges, Range &corner_nodes,
210 const bool update_meshsets = false,
211 const bool debug = false);
212
213 /**
214 * @brief Create front from the surface
215 *
216 * @param debug
217 * @return MoFEMErrorCode
218 */
219 MoFEMErrorCode makeFront(const bool debug = false);
220
221 /**
222 * @brief Calculate distance from mesh nodes to cut surface
223 *
224 * @param intersect_vol
225 * @param verb
226 * @param debug
227 * @return MoFEMErrorCode
228 */
230 const bool debug = false);
231
232 /**
233 * @brief Calculate distance from mesh nodes to surface front
234 *
235 * @param vol
236 * @param th if not null take coordinates from tag
237 * @param verb
238 * @param debug
239 * @return MoFEMErrorCode
240 */
242 int verb = QUIET,
243 const bool debug = false);
244
245 /**
246 * @brief Find level set volumes
247 *
248 * @param th
249 * @param vol_edges
250 * @param remove_adj_prims_edges
251 * @param verb
252 * @param debug
253 * @param edges_file_name
254 * @return MoFEMErrorCode
255 */
257 findLevelSetVolumes(Tag th, Range &vol_edges,
258 const bool remove_adj_prims_edges, int verb = QUIET,
259 const bool debug = false,
260 const std::string edges_file_name = string());
261
262 /**
263 * @brief Find level set volumes
264 *
265 * @param update_front
266 * @param verb
267 * @param debug
268 * @return MoFEMErrorCode
269 */
271 const bool debug = false);
272
273 /**
274 * @brief Refine and set level sets
275 *
276 * \note Should be run befor cutting
277 *
278 * @param refine_front refine nodes at front
279 * @param update_front update level set at front
280 * @param init_bit_level inital bit ref level to store refined meshes
281 * @param surf_levels number of mesh surface refinement
282 * @param front_levels
283 * @param fixed_edges
284 * @param verb
285 * @param debug
286 * @return MoFEMErrorCode
287 */
288 MoFEMErrorCode refineMesh(const int init_bit_level, const int surf_levels,
289 const int front_levels,
290 Range *fixed_edges = nullptr, int verb = QUIET,
291 const bool debug = false);
292
293 /**
294 * @brief find edges to cut
295 *
296 * @param vol is tetrahedrons search to cut
297 * @param verb verbosity level
298 * @param debug debugging
299 * @return MoFEMErrorCode
300 */
302 const bool debug = false);
303
304 /**
305 * @brief Find entities on cut surface which can be projected
306 *
307 * @param fixed_edges pointer to fix edges
308 * @param corner_nodes pointer to corner nodes
309 * @param close tolerance is tolerance how close entities has to be
310 * @param verb verbosity level
311 * @param debug true for debuging purposes
312 *
313 */
315 Range *corner_nodes,
316 const double geometry_tol = 0,
317 const double close_tol = 0,
318 const int verb = QUIET,
319 const bool debug = false);
320
321 /**
322 * \brief cut edges
323 *
324 * For edges to cut (calculated by findEdgesToCut), edges are split in the
325 * middle and then using MoFEM::MeshRefinement interface, tetrahedra mesh
326 * are cut.
327 *
328 * @param bit BitRefLevel of new mesh created by cutting edges
329 * @return error code
330 */
332 Range &cut_surf, Range &cut_verts,
333 const bool debug = false);
334
335 /**
336 * \brief projecting of mid edge nodes on new mesh on surface
337 * @return error code
338 */
340
341 /**
342 * \brief Find edges to trimEdges
343
344 * To make this work, you need to find edges to cut (findEdgesToCut), then
345 * cut edges in the middle (cutEdgesInMiddle) and finally project edges on
346 * the surface (moveMidNodesOnCutEdges)
347
348 * @param verb verbosity level
349 * @return error code
350 */
351 MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes,
352 Tag th = NULL, const double tol = 1e-4,
353 const bool debug = false);
354
355 /**
356 * \brief trim edges
357 * @param bit bit level of the trimmed mesh
358 * @return error code
359 */
361 const bool debug = false);
362
363 /**
364 * \brief move trimmed edges mid nodes
365 * @return error code
366 */
368
369 /**
370 * @brief Trim surface from faces beyond front
371 *
372 * @param fixed_edge
373 * @param corner_nodes
374 * @param debug
375 * @return MoFEMErrorCode
376 */
377 MoFEMErrorCode trimSurface(Range *fixed_edge, Range *corner_nodes,
378 const double tol = 1e-4, Tag th = NULL,
379 const bool debug = false);
380
381 /**
382 * \brief Remove pathological elements on surface internal front
383 *
384 * Internal surface skin is a set of edges in iterior of the body on boundary
385 * of surface. This set of edges is called surface front. If surface face has
386 * three nodes on surface front, non of the face nodes is split and should be
387 * removed from surface if it is going to be split.
388 *
389 * @param split_bit split bit level
390 * @param bit bit level of split mesh
391 * @param ents ents on the surface which is going to be split
392 * @return error code
393 */
395 Range &ents);
396
397 /**
398 * \brief split sides
399 * @param split_bit split bit level
400 * @param bit bit level of split mesh
401 * @param ents ents on the surface which is going to be split
402 * @return error code
403 */
404 MoFEMErrorCode splitSides(const BitRefLevel split_bit, const BitRefLevel bit,
405 const Range &ents, Tag th = NULL);
406
407 /**
408 * @brief Merge edges
409 *
410 * Sort all edges, where sorting is by quality calculated as edge length times
411 * quality of tets adjacent to the edge. Edge is merged if quality if the mesh
412 * is improved.
413 *
414 * @param fraction_level Fraction of edges attemt to be merged at iteration
415 * @param tets Tets of the mesh which edges are merged
416 * @param surface Surface created by edge spliting
417 * @param fixed_edges edges which are geometrical corners of the body
418 * @param corner_nodes vertices on the corners
419 * @param merged_nodes merged nodes
420 * @param out_tets returned test after merge
421 * @param new_surf new surface without merged edges
422 * @param th tag with nodal positons
423 * @param bit_ptr set bit ref level to mesh without merged edges
424 * @param debug
425 * @return MoFEMErrorCode
426 */
427 MoFEMErrorCode mergeBadEdges(const int fraction_level, const Range &tets,
428 const Range &surface, const Range &fixed_edges,
429 const Range &corner_nodes, Range &merged_nodes,
430 Range &out_tets, Range &new_surf, Tag th,
431 const bool update_meshsets = false,
432 const BitRefLevel *bit_ptr = NULL,
433 const bool debug = false);
434
435 /**
436 * @brief Merge edges
437 *
438 * Sort all edges, where sorting is by quality calculated as edge length times
439 * quality of tets adjacent to the edge. Edge is merged if quality if the mesh
440 * is improved.
441 */
443 mergeBadEdges(const int fraction_level, const BitRefLevel cut_bit,
444 const BitRefLevel trim_bit, const BitRefLevel bit,
445 const Range &surface, const Range &fixed_edges,
446 const Range &corner_nodes, Tag th = NULL,
447 const bool update_meshsets = false, const bool debug = false);
448
449 /**
450 * \brief set coords to tag
451 * @param th tag handle
452 * @return error code
453 */
455
456 /**
457 * \brief set coords from tag
458 * @param th tag handle
459 * @return error code
460 */
462 const BitRefLevel mask = BitRefLevel().set());
463
464 inline const Range &getVolume() const { return vOlume; }
465 inline const Range &getSurface() const { return sUrface; }
466 inline const Range &getFront() const { return fRont; }
467
468 inline const Range &getCutEdges() const { return cutEdges; }
469 inline const Range &getNewCutVolumes() const { return cutNewVolumes; }
470 inline const Range &getNewCutSurfaces() const { return cutNewSurfaces; }
471 inline const Range &getNewCutVertices() const { return cutNewVertices; }
472 inline const Range &projectZeroDistanceEnts() const {
473 return zeroDistanceEnts;
474 }
475
476 inline const Range &getTrimEdges() const { return trimEdges; }
477 inline const Range &getNewTrimVolumes() const { return trimNewVolumes; }
478 inline const Range &getNewTrimSurfaces() const { return trimNewSurfaces; }
479 inline const Range &getNewTrimVertices() const { return trimNewVertices; }
480
481 inline const Range &getMergedVolumes() const { return mergedVolumes; }
482 inline const Range &getMergedSurfaces() const { return mergedSurfaces; }
483
484 inline const Range &getCutSurfaceVolumes() const { return cutSurfaceVolumes; }
485 inline const Range &getCutFrontVolumes() const { return cutFrontVolumes; }
486
487 MoFEMErrorCode saveCutEdges(const std::string prefix = "");
488
490
491 inline boost::shared_ptr<OrientedBoxTreeTool> &getTreeSurfPtr() {
492 return treeSurfPtr;
493 }
494
496
497 struct SaveData {
498 moab::Interface &moab;
499 SaveData(moab::Interface &moab) : moab(moab) {}
500 MoFEMErrorCode operator()(const std::string name, const Range &ents) {
502 EntityHandle meshset;
503 CHKERR moab.create_meshset(MESHSET_SET, meshset);
504 CHKERR moab.add_entities(meshset, ents);
505 CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1);
506 CHKERR moab.delete_entities(&meshset, 1);
508 }
509 };
510
512 double lEngth;
513 double qUality;
515 bool skip;
516 LengthMapData(const double l, double q, const EntityHandle e)
517 : lEngth(l), qUality(q), eDge(e), skip(false) {}
518 };
519
520 typedef multi_index_container<
521 LengthMapData,
522 indexed_by<
523
524 ordered_non_unique<
525 member<LengthMapData, double, &LengthMapData::lEngth>>,
526
527 hashed_unique<
528 member<LengthMapData, EntityHandle, &LengthMapData::eDge>>,
529
530 ordered_non_unique<
531 member<LengthMapData, double, &LengthMapData::qUality>>
532
533 >>
535
536private:
540
541 boost::shared_ptr<OrientedBoxTreeTool> treeSurfPtr;
543
550
555
558
565
566 map<EntityHandle, TreeData> edgesToCut;
567 map<EntityHandle, TreeData> verticesOnCutEdges;
568 map<EntityHandle, TreeData> edgesToTrim;
569 map<EntityHandle, TreeData> verticesOnTrimEdges;
570
571 double aveLength; ///< Average edge length
572 double maxLength; ///< Maximal edge length
573
577};
578} // namespace MoFEM
579
580#endif //__CUTMESHINTERFACE_HPP__
581
582/**
583 * \defgroup mesh_cut CutMeshInterface
584 * \brief Interface to mesh cut mesh
585 *
586 * \ingroup mofem
587 */
MoFEM interface.
@ QUIET
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto bit
set bit
FTensor::Index< 'l', 3 > l
double tol
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
Definition Types.hpp:92
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
static const bool debug
Core (interface) class.
Definition Core.hpp:82
LengthMapData(const double l, double q, const EntityHandle e)
MoFEMErrorCode operator()(const std::string name, const Range &ents)
Interface to cut meshes.
MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, const bool debug=false)
Find edges to trimEdges.
MoFEMErrorCode copySurface(const Range surface, Tag th=NULL, double *shift=NULL, double *origin=NULL, double *transform=NULL, const std::string save_mesh="")
copy surface entities
const Range & getNewCutSurfaces() const
MoFEMErrorCode setFront(const Range surface)
MoFEMErrorCode findLevelSetVolumes(Tag th, Range &vol_edges, const bool remove_adj_prims_edges, int verb=QUIET, const bool debug=false, const std::string edges_file_name=string())
Find level set volumes.
MoFEMErrorCode moveMidNodesOnTrimmedEdges(Tag th=NULL)
move trimmed edges mid nodes
const Range & getFront() const
MoFEMErrorCode trimOnly(const BitRefLevel trim_bit, Tag th, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Trim mesh only.
const Range & getCutEdges() const
MoFEMErrorCode getOptions()
Get options from command line.
MoFEMErrorCode trimSurface(Range *fixed_edge, Range *corner_nodes, const double tol=1e-4, Tag th=NULL, const bool debug=false)
Trim surface from faces beyond front.
const Range & getVolume() const
double aveLength
Average edge length.
const Range & getTrimEdges() const
MoFEMErrorCode trimEdgesInTheMiddle(const BitRefLevel bit, const bool debug=false)
trim edges
const Range & getNewCutVolumes() const
MoFEMErrorCode setCoords(Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
set coords from tag
const Range & projectZeroDistanceEnts() const
const Range & getCutFrontVolumes() const
MoFEMErrorCode findEdgesToCut(Range vol, int verb=QUIET, const bool debug=false)
find edges to cut
const Range & getMergedSurfaces() const
const Range & getNewTrimVolumes() const
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
MoFEMErrorCode cutEdgesInMiddle(const BitRefLevel bit, Range &cut_vols, Range &cut_surf, Range &cut_verts, const bool debug=false)
cut edges
MoFEMErrorCode setSurface(const Range surface)
set surface entities
map< EntityHandle, TreeData > edgesToTrim
map< EntityHandle, TreeData > verticesOnTrimEdges
MoFEMErrorCode buildTree()
build tree
double maxLength
Maximal edge length.
MoFEMErrorCode splitSides(const BitRefLevel split_bit, const BitRefLevel bit, const Range &ents, Tag th=NULL)
split sides
MoFEMErrorCode mergeBadEdges(const int fraction_level, const Range &tets, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Range &merged_nodes, Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets=false, const BitRefLevel *bit_ptr=NULL, const bool debug=false)
Merge edges.
MoFEMErrorCode mergeVolumes(const Range volume)
merge volume entities
MoFEMErrorCode getSubInterfaceOptions()
boost::shared_ptr< OrientedBoxTreeTool > & getTreeSurfPtr()
MoFEMErrorCode moveMidNodesOnCutEdges(Tag th=NULL)
projecting of mid edge nodes on new mesh on surface
map< EntityHandle, TreeData > verticesOnCutEdges
MoFEMErrorCode saveCutEdges(const std::string prefix="")
MoFEMErrorCode snapSurfaceSkinToEdges(const Range fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
const Range & getNewTrimVertices() const
const Range & getSurface() const
MoFEMErrorCode removePathologicalFrontTris(const BitRefLevel split_bit, Range &ents)
Remove pathological elements on surface internal front.
MoFEMErrorCode setConstrainSurface(const Range surf)
Set the constrain surface object.
MoFEMErrorCode cutAndTrim(int &first_bit, Tag th, const double tol_geometry, const double tol_cut_close, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Cut and trim.
const Range & getNewTrimSurfaces() const
MoFEMErrorCode cutOnly(Range vol, const BitRefLevel cut_bit, Tag th, const double tol_geometry, const double tol_cut_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Cut mesh only.
MoFEMErrorCode setVolume(const Range volume)
set volume entities
const Range & getMergedVolumes() const
MoFEMErrorCode createFrontLevelSets(Range vol, Tag th=nullptr, int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to surface front.
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
CutMeshInterface(const MoFEM::Core &core)
MoFEMErrorCode cutTrimAndMerge(int &first_bit, const int fraction_level, Tag th, const double tol_geometry, const double tol_cut_close, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)
Cut, trim and merge.
map< EntityHandle, TreeData > edgesToCut
MoFEMErrorCode makeFront(const bool debug=false)
Create front from the surface.
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
MoFEMErrorCode mergeSurface(const Range surface)
merge surface entities
MoFEMErrorCode refineMesh(const int init_bit_level, const int surf_levels, const int front_levels, Range *fixed_edges=nullptr, int verb=QUIET, const bool debug=false)
Refine and set level sets.
const Range & getNewCutVertices() const
MoFEMErrorCode snapSurfaceToEdges(const Range surface_edges, const Range fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
multi_index_container< LengthMapData, indexed_by< ordered_non_unique< member< LengthMapData, double, &LengthMapData::lEngth > >, hashed_unique< member< LengthMapData, EntityHandle, &LengthMapData::eDge > >, ordered_non_unique< member< LengthMapData, double, &LengthMapData::qUality > > > > LengthMapData_multi_index
const Range & getCutSurfaceVolumes() const
MoFEMErrorCode createSurfaceLevelSets(int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to cut surface.
base class for all interface classes