v0.14.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, UnknownInterface **iface) const;
22
24 CutMeshInterface(const MoFEM::Core &core);
26
30
32
33 /**
34 * \brief Get options from command line
35 * @return error code
36 */
39 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "cut_", "MOFEM Cut mesh options",
40 "none");
41
42 CHKERR PetscOptionsInt("-linesearch_steps",
43 "number of bisection steps which line search do to "
44 "find optimal merged nodes position",
45 "", lineSearchSteps, &lineSearchSteps, PETSC_NULL);
46
47 CHKERR PetscOptionsInt("-max_merging_cycles",
48 "number of maximal merging cycles", "",
50
51 CHKERR PetscOptionsScalar("-project_entities_quality_trashold",
52 "project entities quality trashold", "",
55
56 ierr = 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 */
75 MoFEMErrorCode copySurface(const Range surface, Tag th = NULL,
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 /**
89 * @brief Set the constrain surface object
90 *
91 * Add surfaces which are restricted by mesh cut. Example of surface which
92 * needs to be respected is an interface, the boundary between two materials,
93 * or crack surface.
94 *
95 * @param surf
96 * @return MoFEMErrorCode
97 */
99
100 /**
101 * \brief merge surface entities
102 * @param surface entities which going to be added
103 * @return error code
104 */
106
107 /**
108 * \brief merge volume entities
109 * @param volume entities which going to be added
110 * @return error code
111 */
112 MoFEMErrorCode mergeVolumes(const Range volume);
113
115 const double rel_tol,
116 const double abs_tol, Tag th = nullptr,
117 const bool debug = false);
118
119 MoFEMErrorCode snapSurfaceToEdges(const Range surface_edges,
120 const Range fixed_edges,
121 const double rel_tol, const double abs_tol,
122 Tag th = nullptr, const bool debug = false);
123
124 /**
125 * \brief build tree
126 * @return error code
127 */
129
130 /**
131 * @brief Cut mesh only
132 *
133 * @param vol
134 * @param cut_bit
135 * @param th
136 * @param tol_geometry
137 * @param tol_cut_close
138 * @param fixed_edges
139 * @param corner_nodes
140 * @param update_meshsets
141 * @param debug
142 * @return MoFEMErrorCode
143 */
144 MoFEMErrorCode cutOnly(Range vol, const BitRefLevel cut_bit, Tag th,
145 const double tol_geometry, const double tol_cut_close,
146 Range *fixed_edges = NULL, Range *corner_nodes = NULL,
147 const bool update_meshsets = false,
148 const bool debug = false);
149
150 /**
151 * @brief Trim mesh only
152 *
153 * @param trim_bit
154 * @param th
155 * @param tol_cut_close
156 * @param fixed_edges
157 * @param corner_nodes
158 * @param update_meshsets
159 * @param debug
160 * @return MoFEMErrorCode
161 */
162 MoFEMErrorCode trimOnly(const BitRefLevel trim_bit, Tag th,
163 const double tol_trim_close,
164 Range *fixed_edges = NULL, Range *corner_nodes = NULL,
165 const bool update_meshsets = false,
166 const bool debug = false);
167
168 /**
169 * @brief Cut and trim
170 *
171 * @param first_bit
172 * @param th
173 * @param tol_geometry
174 * @param tol_cut_close
175 * @param tol_trim_close
176 * @param fixed_edges
177 * @param corner_nodes
178 * @param update_meshsets
179 * @param debug
180 * @return MoFEMErrorCode
181 */
183 cutAndTrim(int &first_bit, Tag th, const double tol_geometry,
184 const double tol_cut_close, const double tol_trim_close,
185 Range *fixed_edges = NULL, Range *corner_nodes = NULL,
186 const bool update_meshsets = false, const bool debug = false);
187
188 /**
189 * @brief Cut, trim and merge
190 *
191 * @param first_bit first bit of bit revel, subsequent set bits are for trim
192 * and merge
193 * @param fraction_level fraction of edges merged at each merge step
194 * @param th tag storring mesh node positions
195 * @param tol_geometry tolerance how mesh node should be close to cut surface (mesh
196 * node is moved), should be small
197 * @param tol_cut_close how crack node should be close to mesh (cut surface
198 * node is moved), can be big
199 * @param tol_trim_close how front node should be close to mesh, can be big
200 * @param fixed_edges edges on which nodes can not be moved
201 * @param corner_nodes nodes which can not be moved
202 * @param update_meshsets update meshsets by parents
203 * @param debug swich on debugging
204 * @return MoFEMErrorCode
205 */
206 MoFEMErrorCode cutTrimAndMerge(int &first_bit, const int fraction_level,
207 Tag th, const double tol_geometry,
208 const double tol_cut_close,
209 const double tol_trim_close,
210 Range &fixed_edges, Range &corner_nodes,
211 const bool update_meshsets = false,
212 const bool debug = false);
213
214 /**
215 * @brief Create front from the surface
216 *
217 * @param debug
218 * @return MoFEMErrorCode
219 */
220 MoFEMErrorCode makeFront(const bool debug = false);
221
222 /**
223 * @brief Calculate distance from mesh nodes to cut surface
224 *
225 * @param intersect_vol
226 * @param verb
227 * @param debug
228 * @return MoFEMErrorCode
229 */
231 const bool debug = false);
232
233 /**
234 * @brief Calculate distance from mesh nodes to surface front
235 *
236 * @param vol
237 * @param th if not null take coordinates from tag
238 * @param verb
239 * @param debug
240 * @return MoFEMErrorCode
241 */
242 MoFEMErrorCode createFrontLevelSets(Range vol, Tag th = nullptr,
243 int verb = QUIET,
244 const bool debug = false);
245
246 /**
247 * @brief Find level set volumes
248 *
249 * @param th
250 * @param vol_edges
251 * @param remove_adj_prims_edges
252 * @param verb
253 * @param debug
254 * @param edges_file_name
255 * @return MoFEMErrorCode
256 */
258 const bool remove_adj_prims_edges,
259 int verb = QUIET, 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,
353 const double tol = 1e-4,
354 const bool debug = false);
355
356 /**
357 * \brief trim edges
358 * @param bit bit level of the trimmed mesh
359 * @return error code
360 */
362 const bool debug = false);
363
364 /**
365 * \brief move trimmed edges mid nodes
366 * @return error code
367 */
369
370 /**
371 * @brief Trim surface from faces beyond front
372 *
373 * @param fixed_edge
374 * @param corner_nodes
375 * @param debug
376 * @return MoFEMErrorCode
377 */
378 MoFEMErrorCode trimSurface(Range *fixed_edge, Range *corner_nodes,
379 const double tol = 1e-4, Tag th = NULL,
380 const bool debug = false);
381
382 /**
383 * \brief Remove pathological elements on surface internal front
384 *
385 * Internal surface skin is a set of edges in iterior of the body on boundary
386 * of surface. This set of edges is called surface front. If surface face has
387 * three nodes on surface front, non of the face nodes is split and should be
388 * removed from surface if it is going to be split.
389 *
390 * @param split_bit split bit level
391 * @param bit bit level of split mesh
392 * @param ents ents on the surface which is going to be split
393 * @return error code
394 */
396 Range &ents);
397
398 /**
399 * \brief split sides
400 * @param split_bit split bit level
401 * @param bit bit level of split mesh
402 * @param ents ents on the surface which is going to be split
403 * @return error code
404 */
405 MoFEMErrorCode splitSides(const BitRefLevel split_bit, const BitRefLevel bit,
406 const Range &ents, Tag th = NULL);
407
408 /**
409 * @brief Merge edges
410 *
411 * Sort all edges, where sorting is by quality calculated as edge length times
412 * quality of tets adjacent to the edge. Edge is merged if quality if the mesh
413 * is improved.
414 *
415 * @param fraction_level Fraction of edges attemt to be merged at iteration
416 * @param tets Tets of the mesh which edges are merged
417 * @param surface Surface created by edge spliting
418 * @param fixed_edges edges which are geometrical corners of the body
419 * @param corner_nodes vertices on the corners
420 * @param merged_nodes merged nodes
421 * @param out_tets returned test after merge
422 * @param new_surf new surface without merged edges
423 * @param th tag with nodal positons
424 * @param bit_ptr set bit ref level to mesh without merged edges
425 * @param debug
426 * @return MoFEMErrorCode
427 */
428 MoFEMErrorCode mergeBadEdges(const int fraction_level, const Range &tets,
429 const Range &surface, const Range &fixed_edges,
430 const Range &corner_nodes, Range &merged_nodes,
431 Range &out_tets, Range &new_surf, Tag th,
432 const bool update_meshsets = false,
433 const BitRefLevel *bit_ptr = NULL,
434 const bool debug = false);
435
436 /**
437 * @brief Merge edges
438 *
439 * Sort all edges, where sorting is by quality calculated as edge length times
440 * quality of tets adjacent to the edge. Edge is merged if quality if the mesh
441 * is improved.
442 */
444 mergeBadEdges(const int fraction_level, const BitRefLevel cut_bit,
445 const BitRefLevel trim_bit, const BitRefLevel bit,
446 const Range &surface, const Range &fixed_edges,
447 const Range &corner_nodes, Tag th = NULL,
448 const bool update_meshsets = false, const bool debug = false);
449
450 /**
451 * \brief set coords to tag
452 * @param th tag handle
453 * @return error code
454 */
456
457 /**
458 * \brief set coords from tag
459 * @param th tag handle
460 * @return error code
461 */
463 const BitRefLevel mask = BitRefLevel().set());
464
465 inline const Range &getVolume() const { return vOlume; }
466 inline const Range &getSurface() const { return sUrface; }
467 inline const Range &getFront() const { return fRont; }
468
469 inline const Range &getCutEdges() const { return cutEdges; }
470 inline const Range &getNewCutVolumes() const { return cutNewVolumes; }
471 inline const Range &getNewCutSurfaces() const { return cutNewSurfaces; }
472 inline const Range &getNewCutVertices() const { return cutNewVertices; }
473 inline const Range &projectZeroDistanceEnts() const {
474 return zeroDistanceEnts;
475 }
476
477 inline const Range &getTrimEdges() const { return trimEdges; }
478 inline const Range &getNewTrimVolumes() const { return trimNewVolumes; }
479 inline const Range &getNewTrimSurfaces() const { return trimNewSurfaces; }
480 inline const Range &getNewTrimVertices() const { return trimNewVertices; }
481
482 inline const Range &getMergedVolumes() const { return mergedVolumes; }
483 inline const Range &getMergedSurfaces() const { return mergedSurfaces; }
484
485 inline const Range &getCutSurfaceVolumes() const { return cutSurfaceVolumes; }
486 inline const Range &getCutFrontVolumes() const { return cutFrontVolumes; }
487
488 MoFEMErrorCode saveCutEdges(const std::string prefix = "");
489
491
492 inline boost::shared_ptr<OrientedBoxTreeTool> &getTreeSurfPtr() {
493 return treeSurfPtr;
494 }
495
497
498 struct SaveData {
499 moab::Interface &moab;
500 SaveData(moab::Interface &moab) : moab(moab) {}
501 MoFEMErrorCode operator()(const std::string name, const Range &ents) {
503 EntityHandle meshset;
504 CHKERR moab.create_meshset(MESHSET_SET, meshset);
505 CHKERR moab.add_entities(meshset, ents);
506 CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1);
507 CHKERR moab.delete_entities(&meshset, 1);
509 }
510 };
511
513 double lEngth;
514 double qUality;
516 bool skip;
517 LengthMapData(const double l, double q, const EntityHandle e)
518 : lEngth(l), qUality(q), eDge(e), skip(false) {}
519 };
520
521 typedef multi_index_container<
522 LengthMapData,
523 indexed_by<
524
525 ordered_non_unique<
526 member<LengthMapData, double, &LengthMapData::lEngth>>,
527
528 hashed_unique<
529 member<LengthMapData, EntityHandle, &LengthMapData::eDge>>,
530
531 ordered_non_unique<
532 member<LengthMapData, double, &LengthMapData::qUality>>
533
534 >>
536
537private:
541
542 boost::shared_ptr<OrientedBoxTreeTool> treeSurfPtr;
544
551
556
559
560 struct TreeData {
561 double dIst;
562 double lEngth;
565 };
566
567 map<EntityHandle, TreeData> edgesToCut;
568 map<EntityHandle, TreeData> verticesOnCutEdges;
569 map<EntityHandle, TreeData> edgesToTrim;
570 map<EntityHandle, TreeData> verticesOnTrimEdges;
571
572 double aveLength; ///< Average edge length
573 double maxLength; ///< Maximal edge length
574
578};
579} // namespace MoFEM
580
581#endif //__CUTMESHINTERFACE_HPP__
582
583/**
584 * \defgroup mesh_cut CutMeshInterface
585 * \brief Interface to mesh cut mesh
586 *
587 * \ingroup mofem
588 */
MoFEM interface.
@ QUIET
Definition: definitions.h:208
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
static const bool debug
auto bit
set bit
FTensor::Index< 'l', 3 > l
double tol
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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
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
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
MoFEMErrorCode saveCutEdges(const std::string prefix="")
MoFEMErrorCode saveTrimEdges()
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
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)
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