v0.14.0
Loading...
Searching...
No Matches
Classes | Public Types | Public Member Functions | Public Attributes | Private Attributes | List of all members
MoFEM::CutMeshInterface Struct Reference

Interface to cut meshes. More...

#include <src/interfaces/CutMeshInterface.hpp>

Inheritance diagram for MoFEM::CutMeshInterface:
[legend]
Collaboration diagram for MoFEM::CutMeshInterface:
[legend]

Classes

struct  LengthMapData
 
struct  SaveData
 
struct  TreeData
 

Public Types

typedef 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
 

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 
 CutMeshInterface (const MoFEM::Core &core)
 
 ~CutMeshInterface ()
 
MoFEMErrorCode getSubInterfaceOptions ()
 
MoFEMErrorCode getOptions ()
 Get options from command line. More...
 
MoFEMErrorCode setFront (const Range surface)
 
MoFEMErrorCode setSurface (const Range surface)
 set surface entities More...
 
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 More...
 
MoFEMErrorCode setVolume (const Range volume)
 set volume entities More...
 
MoFEMErrorCode setConstrainSurface (const Range surf)
 Set the constrain surface object. More...
 
MoFEMErrorCode mergeSurface (const Range surface)
 merge surface entities More...
 
MoFEMErrorCode mergeVolumes (const Range volume)
 merge volume entities More...
 
MoFEMErrorCode snapSurfaceSkinToEdges (const Range fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
 
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)
 
MoFEMErrorCode buildTree ()
 build tree More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
MoFEMErrorCode makeFront (const bool debug=false)
 Create front from the surface. More...
 
MoFEMErrorCode createSurfaceLevelSets (int verb=QUIET, const bool debug=false)
 Calculate distance from mesh nodes to cut surface. More...
 
MoFEMErrorCode createFrontLevelSets (Range vol, Tag th=nullptr, int verb=QUIET, const bool debug=false)
 Calculate distance from mesh nodes to surface front. More...
 
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. More...
 
MoFEMErrorCode findLevelSetVolumes (int verb=QUIET, const bool debug=false)
 Find level set volumes. More...
 
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. More...
 
MoFEMErrorCode findEdgesToCut (Range vol, int verb=QUIET, const bool debug=false)
 find edges to cut More...
 
MoFEMErrorCode projectZeroDistanceEnts (Range *fixed_edges, Range *corner_nodes, const double geometry_tol=0, const double close_tol=0, const int verb=QUIET, const bool debug=false)
 Find entities on cut surface which can be projected. More...
 
MoFEMErrorCode cutEdgesInMiddle (const BitRefLevel bit, Range &cut_vols, Range &cut_surf, Range &cut_verts, const bool debug=false)
 cut edges More...
 
MoFEMErrorCode moveMidNodesOnCutEdges (Tag th=NULL)
 projecting of mid edge nodes on new mesh on surface More...
 
MoFEMErrorCode findEdgesToTrim (Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, const bool debug=false)
 Find edges to trimEdges. More...
 
MoFEMErrorCode trimEdgesInTheMiddle (const BitRefLevel bit, const bool debug=false)
 trim edges More...
 
MoFEMErrorCode moveMidNodesOnTrimmedEdges (Tag th=NULL)
 move trimmed edges mid nodes More...
 
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. More...
 
MoFEMErrorCode removePathologicalFrontTris (const BitRefLevel split_bit, Range &ents)
 Remove pathological elements on surface internal front. More...
 
MoFEMErrorCode splitSides (const BitRefLevel split_bit, const BitRefLevel bit, const Range &ents, Tag th=NULL)
 split sides More...
 
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. More...
 
MoFEMErrorCode mergeBadEdges (const int fraction_level, const BitRefLevel cut_bit, const BitRefLevel trim_bit, const BitRefLevel bit, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Tag th=NULL, const bool update_meshsets=false, const bool debug=false)
 Merge edges. More...
 
MoFEMErrorCode setTagData (Tag th, const BitRefLevel bit=BitRefLevel())
 set coords to tag More...
 
MoFEMErrorCode setCoords (Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
 set coords from tag More...
 
const RangegetVolume () const
 
const RangegetSurface () const
 
const RangegetFront () const
 
const RangegetCutEdges () const
 
const RangegetNewCutVolumes () const
 
const RangegetNewCutSurfaces () const
 
const RangegetNewCutVertices () const
 
const RangeprojectZeroDistanceEnts () const
 
const RangegetTrimEdges () const
 
const RangegetNewTrimVolumes () const
 
const RangegetNewTrimSurfaces () const
 
const RangegetNewTrimVertices () const
 
const RangegetMergedVolumes () const
 
const RangegetMergedSurfaces () const
 
const RangegetCutSurfaceVolumes () const
 
const RangegetCutFrontVolumes () const
 
MoFEMErrorCode saveCutEdges (const std::string prefix="")
 
MoFEMErrorCode saveTrimEdges ()
 
boost::shared_ptr< OrientedBoxTreeTool > & getTreeSurfPtr ()
 
MoFEMErrorCode clearMap ()
 
- Public Member Functions inherited from MoFEM::UnknownInterface
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

Public Attributes

MoFEM::CorecOre
 
int lineSearchSteps
 
int nbMaxMergingCycles
 
double projectEntitiesQualityTrashold
 

Private Attributes

Range fRont
 
Range sUrface
 
Range vOlume
 
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
 
EntityHandle rootSetSurf
 
Range cutEdges
 
Range cutNewVolumes
 
Range cutNewSurfaces
 
Range zeroDistanceEnts
 
Range zeroDistanceVerts
 
Range cutNewVertices
 
Range trimNewVolumes
 
Range trimNewVertices
 
Range trimNewSurfaces
 
Range trimEdges
 
Range mergedVolumes
 
Range mergedSurfaces
 
map< EntityHandle, TreeDataedgesToCut
 
map< EntityHandle, TreeDataverticesOnCutEdges
 
map< EntityHandle, TreeDataedgesToTrim
 
map< EntityHandle, TreeDataverticesOnTrimEdges
 
double aveLength
 Average edge length. More...
 
double maxLength
 Maximal edge length. More...
 
Range cutSurfaceVolumes
 
Range cutFrontVolumes
 
Range constrainSurface
 

Additional Inherited Members

- Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode getLibVersion (Version &version)
 Get library version. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

Detailed Description

Interface to cut meshes.

Examples
mesh_cut.cpp.

Definition at line 19 of file CutMeshInterface.hpp.

Member Typedef Documentation

◆ LengthMapData_multi_index

typedef 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> > > > MoFEM::CutMeshInterface::LengthMapData_multi_index

Definition at line 535 of file CutMeshInterface.hpp.

Constructor & Destructor Documentation

◆ CutMeshInterface()

MoFEM::CutMeshInterface::CutMeshInterface ( const MoFEM::Core core)

Definition at line 23 of file CutMeshInterface.cpp.

24 : cOre(const_cast<Core &>(core)) {
25 lineSearchSteps = 10;
28}
CoreTmp< 0 > Core
Definition: Core.hpp:1094

◆ ~CutMeshInterface()

MoFEM::CutMeshInterface::~CutMeshInterface ( )
inline

Definition at line 25 of file CutMeshInterface.hpp.

25{}

Member Function Documentation

◆ buildTree()

MoFEMErrorCode MoFEM::CutMeshInterface::buildTree ( )

build tree

Returns
error code
Examples
mesh_cut.cpp.

Definition at line 215 of file CutMeshInterface.cpp.

215 {
216 CoreInterface &m_field = cOre;
217 moab::Interface &moab = m_field.get_moab();
219 treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
220 new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
223}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#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
moab::Interface & get_moab()
Definition: Core.hpp:322
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr

◆ clearMap()

MoFEMErrorCode MoFEM::CutMeshInterface::clearMap ( )

Definition at line 30 of file CutMeshInterface.cpp.

30 {
32 treeSurfPtr.reset();
34}

◆ copySurface()

MoFEMErrorCode MoFEM::CutMeshInterface::copySurface ( const Range  surface,
Tag  th = NULL,
double shift = NULL,
double origin = NULL,
double transform = NULL,
const std::string  save_mesh = "" 
)

copy surface entities

Parameters
surfaceentities which going to be added
Returns
error code
Examples
mesh_cut.cpp.

Definition at line 48 of file CutMeshInterface.cpp.

51 {
52 CoreInterface &m_field = cOre;
53 moab::Interface &moab = m_field.get_moab();
55 sUrface.clear();
56 std::map<EntityHandle, EntityHandle> verts_map;
57 for (Range::const_iterator tit = surface.begin(); tit != surface.end();
58 tit++) {
59 int num_nodes;
60 const EntityHandle *conn;
61 CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
62 MatrixDouble coords(num_nodes, 3);
63 if (th) {
64 CHKERR moab.tag_get_data(th, conn, num_nodes, &coords(0, 0));
65 } else {
66 CHKERR moab.get_coords(conn, num_nodes, &coords(0, 0));
67 }
68 EntityHandle new_verts[num_nodes];
69 for (int nn = 0; nn != num_nodes; nn++) {
70 if (verts_map.find(conn[nn]) != verts_map.end()) {
71 new_verts[nn] = verts_map[conn[nn]];
72 } else {
73 if (transform) {
74 ublas::matrix_row<MatrixDouble> mr(coords, nn);
75 if (origin) {
76 VectorAdaptor vec_origin(
77 3, ublas::shallow_array_adaptor<double>(3, origin));
78 mr = mr - vec_origin;
79 }
80 MatrixAdaptor mat_transform = MatrixAdaptor(
81 3, 3, ublas::shallow_array_adaptor<double>(9, transform));
82 mr = prod(mat_transform, mr);
83 if (origin) {
84 VectorAdaptor vec_origin(
85 3, ublas::shallow_array_adaptor<double>(3, origin));
86 mr = mr + vec_origin;
87 }
88 }
89 if (shift) {
90 ublas::matrix_row<MatrixDouble> mr(coords, nn);
91 VectorAdaptor vec_shift(
92 3, ublas::shallow_array_adaptor<double>(3, shift));
93 mr = mr + vec_shift;
94 }
95 CHKERR moab.create_vertex(&coords(nn, 0), new_verts[nn]);
96 verts_map[conn[nn]] = new_verts[nn];
97 }
98 }
99 EntityHandle ele;
100 CHKERR moab.create_element(MBTRI, new_verts, num_nodes, ele);
101 sUrface.insert(ele);
102 }
103 if (!save_mesh.empty())
104 CHKERR SaveData(m_field.get_moab())(save_mesh, sUrface);
106}
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:132

◆ createFrontLevelSets()

MoFEMErrorCode MoFEM::CutMeshInterface::createFrontLevelSets ( Range  vol,
Tag  th = nullptr,
int  verb = QUIET,
const bool  debug = false 
)

Calculate distance from mesh nodes to surface front.

Parameters
vol
thif not null take coordinates from tag
verb
debug
Returns
MoFEMErrorCode

Definition at line 609 of file CutMeshInterface.cpp.

611 {
612 CoreInterface &m_field = cOre;
613 moab::Interface &moab = m_field.get_moab();
615
616 auto create_tag = [&](const std::string name, const int dim) {
617 Tag th;
618 rval = moab.tag_get_handle(name.c_str(), th);
619 if (rval == MB_SUCCESS)
620 return th;
621 std::vector<double> def_val(dim, 0);
622 CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
623 MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
624 return th;
625 };
626
627 Range vol_vertices;
628 CHKERR moab.get_connectivity(vol, vol_vertices, true);
629
630 std::vector<double> min_distances_from_front(vol_vertices.size(), -1);
631 std::vector<double> points_on_edges(3 * vol_vertices.size(), 0);
632 std::vector<EntityHandle> closest_edges(vol_vertices.size(), 0);
633
634 std::vector<double> coords(3 * vol_vertices.size());
635 if (th)
636 CHKERR moab.tag_get_data(th, vol_vertices, &*coords.begin());
637 else
638 CHKERR moab.get_coords(vol_vertices, &*coords.begin());
639
640 CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
641 &*coords.begin(), vol_vertices.size(), fRont,
642 &*min_distances_from_front.begin(), &*points_on_edges.begin(),
643 &*closest_edges.begin());
644
645 if (!points_on_edges.empty()) {
646 for (int i = 0; i != min_distances_from_front.size(); ++i) {
647 Range faces;
648 CHKERR moab.get_adjacencies(&closest_edges[0], 1, 2, false, faces);
649 auto point_in = getVectorAdaptor(&coords[3 * i], 3);
650 auto point_out = getVectorAdaptor(&points_on_edges[3 * i], 3);
651 point_out -= point_in;
652 }
653 }
654
655 auto th_dist_front_vec = create_tag("DIST_FRONT_VECTOR", 3);
656 CHKERR moab.tag_set_data(th_dist_front_vec, vol_vertices,
657 &*points_on_edges.begin());
658
660}
#define CutMeshFunctionBegin
const int dim
FTensor::Index< 'i', SPACE_DIM > i
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ createSurfaceLevelSets()

MoFEMErrorCode MoFEM::CutMeshInterface::createSurfaceLevelSets ( int  verb = QUIET,
const bool  debug = false 
)

Calculate distance from mesh nodes to cut surface.

Parameters
intersect_vol
verb
debug
Returns
MoFEMErrorCode

Definition at line 477 of file CutMeshInterface.cpp.

478 {
479 CoreInterface &m_field = cOre;
480 moab::Interface &moab = m_field.get_moab();
482 auto tools_interface = m_field.getInterface<Tools>();
483
484 auto create_tag = [&](const std::string name, const int dim) {
485 Tag th;
486 rval = moab.tag_get_handle(name.c_str(), th);
487 if (rval == MB_SUCCESS)
488 return th;
489 std::vector<double> def_val(dim, 0);
490 CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
491 MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
492
493 return th;
494 };
495
496 auto set_vol = [&](const Range &vol_verts, std::vector<double> &coords,
497 std::vector<double> &dist_surface_vec,
498 std::vector<double> &dist_surface_normal_vec,
499 std::vector<double> &dist_surface_signed_dist_vec) {
501
502 coords.resize(3 * vol_verts.size());
503 dist_surface_vec.resize(3 * vol_verts.size());
504 dist_surface_normal_vec.resize(3 * vol_verts.size());
505 dist_surface_signed_dist_vec.resize(vol_verts.size());
506
507 CHKERR moab.get_coords(vol_verts, &*coords.begin());
508 std::srand(0);
509
510 for (auto v : vol_verts) {
511
512 const int index = vol_verts.index(v);
513 auto point_in = getVectorAdaptor(&coords[3 * index], 3);
514 VectorDouble3 point_out(3);
515 EntityHandle facets_out;
516 CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
517 &point_out[0], facets_out);
518
519 VectorDouble3 delta = point_out - point_in;
520 VectorDouble3 n_first(3);
521 CHKERR tools_interface->getTriNormal(facets_out, &*n_first.begin());
522
523 // Check of only triangle in proximity of point out is one triangle, if is
524 // more than one, use normal of lager triangle to set orientation.
525 auto check_triangle_orientation = [&](auto n) {
526 int num_nodes;
527 const EntityHandle *conn;
528 CHKERR moab.get_connectivity(facets_out, conn, num_nodes, true);
529 MatrixDouble3by3 coords(3, 3);
530 CHKERR moab.get_coords(conn, 3, &coords(0, 0));
531 VectorDouble3 center(3);
532 center.clear();
533 for (auto ii : {0, 1, 2})
534 for (auto jj : {0, 1, 2})
535 center(jj) += coords(ii, jj) / 3;
536
537 std::vector<EntityHandle> triangles_out;
538 std::vector<double> distance_out;
539 auto a_max = norm_2(n);
540 VectorDouble3 ray = -n / a_max;
541 VectorDouble3 pt = center - ray * sqrt(a_max);
542 const double ray_length = 2 * sqrt(a_max);
543
544 CHKERR treeSurfPtr->ray_intersect_triangles(
545 distance_out, triangles_out, rootSetSurf,
546 std::numeric_limits<float>::epsilon(), &pt[0], &ray[0],
547 &ray_length);
548
549 if (triangles_out.size() > 1) {
550 VectorDouble3 nt(3);
551 for (auto t : triangles_out) {
552 CHKERR tools_interface->getTriNormal(t, &*nt.begin());
553 double at = norm_2(nt);
554 if (at > a_max) {
555 n = nt;
556 a_max = at;
557 }
558 }
559 }
560
561 return n;
562 };
563
564 auto n = check_triangle_orientation(n_first);
565 n /= norm_2(n);
566 const double dot = inner_prod(delta, n);
567
568
569 if (std::abs(dot) < std::numeric_limits<double>::epsilon()) {
570 if (std::rand() % 2 == 0)
571 delta += n * std::numeric_limits<double>::epsilon();
572 else
573 delta -= n * std::numeric_limits<double>::epsilon();
574 }
575
576 auto dist_vec = getVectorAdaptor(&dist_surface_vec[3 * index], 3);
577 noalias(dist_vec) = delta;
578
579 auto dist_normal_vec =
580 getVectorAdaptor(&dist_surface_normal_vec[3 * index], 3);
581
582 dist_surface_signed_dist_vec[index] = dot;
583 noalias(dist_normal_vec) = dot * n;
584 }
585
587 };
588
589 std::vector<double> coords;
590 std::vector<double> dist_surface_vec;
591 std::vector<double> dist_surface_normal_vec;
592 std::vector<double> dist_surface_signed_dist_vec;
593 Range vol_verts;
594 CHKERR moab.get_connectivity(vOlume, vol_verts, true);
595
596 CHKERR set_vol(vol_verts, coords, dist_surface_vec, dist_surface_normal_vec,
597 dist_surface_signed_dist_vec);
598
599 CHKERR moab.tag_set_data(create_tag("DIST_SURFACE_VECTOR", 3), vol_verts,
600 &*dist_surface_vec.begin());
601 CHKERR moab.tag_set_data(create_tag("DIST_SURFACE_NORMAL_VECTOR", 3),
602 vol_verts, &*dist_surface_normal_vec.begin());
603 CHKERR moab.tag_set_data(create_tag("DIST_SURFACE_NORMAL_SIGNED", 1),
604 vol_verts, &*dist_surface_signed_dist_vec.begin());
605
607}
FTensor::Index< 'n', SPACE_DIM > n
const double v
phase velocity of light in medium (cm/ns)
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
constexpr double t
plate stiffness
Definition: plate.cpp:59
static constexpr double delta

◆ cutAndTrim()

MoFEMErrorCode MoFEM::CutMeshInterface::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.

Parameters
first_bit
th
tol_geometry
tol_cut_close
tol_trim_close
fixed_edges
corner_nodes
update_meshsets
debug
Returns
MoFEMErrorCode

Definition at line 309 of file CutMeshInterface.cpp.

312 {
313 CoreInterface &m_field = cOre;
315
316 std::vector<BitRefLevel> bit_levels;
317
318 auto get_back_bit_levels = [&]() {
319 bit_levels.push_back(BitRefLevel());
320 bit_levels.back().set(first_bit + bit_levels.size() - 1);
321 return bit_levels.back();
322 };
323
324 BitRefLevel cut_bit = get_back_bit_levels();
325
327 tol_geometry, tol_cut_close, fixed_edges, corner_nodes,
328 update_meshsets, debug);
329
330 auto get_min_quality = [&m_field](const BitRefLevel bit, Tag th) {
331 Range tets_level; // test at level
332 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
333 bit, BitRefLevel().set(), MBTET, tets_level);
334 double min_q = 1;
335 CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
336 return min_q;
337 };
338
339 MOFEM_LOG_C("WORLD", Sev::inform, "Min quality cut %6.4g",
340 get_min_quality(cut_bit, th));
341
342 Range starting_volume = cutNewVolumes;
343 Range starting_volume_nodes;
344 CHKERR m_field.get_moab().get_connectivity(starting_volume,
345 starting_volume_nodes, true);
346 std::vector<double> staring_volume_coords(3 * starting_volume_nodes.size());
347 CHKERR m_field.get_moab().get_coords(starting_volume_nodes,
348 &*staring_volume_coords.begin());
349
350 if (th) {
351 std::vector<double> staring_volume_th_coords(3 *
352 starting_volume_nodes.size());
353 CHKERR m_field.get_moab().tag_get_data(th, starting_volume_nodes,
354 &*staring_volume_th_coords.begin());
355 CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
356 &*staring_volume_th_coords.begin());
357 }
358
359 if (th)
361
362 BitRefLevel trim_bit = get_back_bit_levels();
363
364 CHKERR trimOnly(trim_bit, th, tol_trim_close, fixed_edges, corner_nodes,
365 update_meshsets, debug);
366
367 MOFEM_LOG_C("WORLD", Sev::inform, "Min quality trim %3.2g",
368 get_min_quality(trim_bit, th));
369
370 first_bit += bit_levels.size() - 1;
371
372 if (th)
373 CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
374 &*staring_volume_coords.begin());
375
377}
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
static const bool debug
auto bit
set bit
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
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.
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 setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag

◆ cutEdgesInMiddle()

MoFEMErrorCode MoFEM::CutMeshInterface::cutEdgesInMiddle ( const BitRefLevel  bit,
Range cut_vols,
Range cut_surf,
Range cut_verts,
const bool  debug = false 
)

cut edges

For edges to cut (calculated by findEdgesToCut), edges are split in the middle and then using MoFEM::MeshRefinement interface, tetrahedra mesh are cut.

Parameters
bitBitRefLevel of new mesh created by cutting edges
Returns
error code

Definition at line 1298 of file CutMeshInterface.cpp.

1302 {
1303 CoreInterface &m_field = cOre;
1304 moab::Interface &moab = m_field.get_moab();
1305 MeshRefinement *refiner;
1306 auto ref_ents_ptr = m_field.get_ref_ents();
1308
1309 if (cutEdges.size() != edgesToCut.size())
1310 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1311
1312 auto refine_mesh = [&]() {
1314 CHKERR m_field.getInterface(refiner);
1315 CHKERR refiner->addVerticesInTheMiddleOfEdges(cutEdges, bit);
1316 CHKERR refiner->refineTets(vOlume, bit, QUIET, debug);
1318 };
1319
1320 CHKERR refine_mesh();
1321
1322 if (debug)
1323 if (cutEdges.size() != edgesToCut.size())
1324 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1325
1326 cut_vols.clear();
1327 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1328 bit, BitRefLevel().set(), MBTET, cut_vols);
1329 cut_surf.clear();
1330 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1331 bit, bit, MBTRI, cut_surf);
1332
1333 if (debug)
1334 CHKERR SaveData(moab)("cut_surf_from_bit.vtk", cut_surf);
1335
1336 // Find new vertices on cut edges
1337 cut_verts.clear();
1338 cut_verts.merge(zeroDistanceVerts);
1339
1340 for (auto &m : edgesToCut) {
1341 auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1342 boost::make_tuple(MBVERTEX, m.first));
1343 if (vit ==
1344 ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1345 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1346 "No vertex on cut edges, that make no sense");
1347 }
1348 const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1349 if ((ref_ent->getBitRefLevel() & bit).any()) {
1350 EntityHandle vert = ref_ent->getEnt();
1351 cut_verts.insert(vert);
1352 verticesOnCutEdges[vert] = m.second;
1353 } else {
1354 SETERRQ1(
1355 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1356 "Vertex has wrong bit ref level %s",
1357 boost::lexical_cast<std::string>(ref_ent->getBitRefLevel()).c_str());
1358 }
1359 }
1360
1361 // Add zero distance entities faces
1362 Range tets_skin;
1363 Skinner skin(&moab);
1364 CHKERR skin.find_skin(0, cut_vols, false, tets_skin);
1365 tets_skin.merge(constrainSurface);
1366
1367 // At that point cut_surf has all newly created faces, now take all
1368 // nodes on those faces and subtract nodes on cut edges. Faces adjacent to
1369 // nodes which left are not part of surface.
1370 cut_surf.merge(zeroDistanceEnts.subset_by_type(MBTRI));
1371 Range diff_verts;
1372 CHKERR moab.get_connectivity(cut_surf, diff_verts, true);
1373 diff_verts = subtract(diff_verts, cut_verts);
1374
1375 Range subtract_faces;
1376 CHKERR moab.get_adjacencies(diff_verts, 2, false, subtract_faces,
1377 moab::Interface::UNION);
1378 cut_surf = subtract(cut_surf, unite(subtract_faces, tets_skin));
1379
1380 cut_verts.clear();
1381 CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1382
1383 // Check non-mainfolds
1384 auto check_for_non_minfold = [&]() {
1386 Range surf_edges;
1387 CHKERR moab.get_adjacencies(cut_surf, 1, false, surf_edges,
1388 moab::Interface::UNION);
1389 for (auto e : surf_edges) {
1390
1391 Range faces;
1392 CHKERR moab.get_adjacencies(&e, 1, 2, false, faces);
1393 faces = intersect(faces, cut_surf);
1394 if (faces.size() > 2) {
1395
1396 bool resolved = false;
1397
1398 // Check for haning node
1399 Range nodes;
1400 CHKERR moab.get_connectivity(faces, nodes, true);
1401 for (auto n : nodes) {
1402 Range adj_faces;
1403 CHKERR moab.get_adjacencies(&n, 1, 2, false, adj_faces);
1404 adj_faces = intersect(adj_faces, cut_surf);
1405 if (adj_faces.size() == 1) {
1406 cut_surf.erase(adj_faces[0]);
1407 resolved = true;
1408 }
1409 }
1410
1411 // Check for two edges minfold
1412 Range adj_edges;
1413 CHKERR moab.get_adjacencies(faces, 1, false, adj_edges,
1414 moab::Interface::UNION);
1415 adj_edges = intersect(adj_edges, surf_edges);
1416 adj_edges.erase(e);
1417 for (auto other_e : adj_edges) {
1418 Range other_faces;
1419 CHKERR moab.get_adjacencies(&other_e, 1, 2, false, other_faces);
1420 other_faces = intersect(other_faces, cut_surf);
1421 if (other_faces.size() > 2) {
1422 other_faces = intersect(other_faces, faces);
1423 cut_surf = subtract(cut_surf, other_faces);
1424 resolved = true;
1425 }
1426 }
1427
1428 if (!resolved && !debug)
1429 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1430 "Non-mainfold surface");
1431
1432 cut_verts.clear();
1433 CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1434 }
1435 }
1437 };
1438
1439 CHKERR check_for_non_minfold();
1440
1441 if (debug)
1442 CHKERR SaveData(moab)("cut_surf.vtk", cut_surf);
1443
1445}
@ QUIET
Definition: definitions.h:208
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FTensor::Index< 'm', SPACE_DIM > m
map< EntityHandle, TreeData > verticesOnCutEdges
map< EntityHandle, TreeData > edgesToCut

◆ cutOnly()

MoFEMErrorCode MoFEM::CutMeshInterface::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.

Parameters
vol
cut_bit
th
tol_geometry
tol_cut_close
fixed_edges
corner_nodes
update_meshsets
debug
Returns
MoFEMErrorCode

Definition at line 226 of file CutMeshInterface.cpp.

229 {
230 CoreInterface &m_field = cOre;
231 moab::Interface &moab = m_field.get_moab();
233
234 // cut mesh
236 CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_geometry,
237 tol_cut_close, QUIET, debug);
240
241 CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(constrainSurface,
243 if (fixed_edges)
244 CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(*fixed_edges,
245 *fixed_edges);
246 if (corner_nodes)
247 CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(*corner_nodes,
248 *corner_nodes);
249 if (update_meshsets)
250 CHKERR m_field.getInterface<MeshsetsManager>()
251 ->updateAllMeshsetsByEntitiesChildren(cut_bit);
253
254 if (debug) {
256 if (fixed_edges)
257 CHKERR SaveData(moab)("out_cut_new_fixed_edges.vtk", *fixed_edges);
258 }
259
261}
const Range & projectZeroDistanceEnts() const
MoFEMErrorCode findEdgesToCut(Range vol, int verb=QUIET, const bool debug=false)
find edges to cut
MoFEMErrorCode cutEdgesInMiddle(const BitRefLevel bit, Range &cut_vols, Range &cut_surf, Range &cut_verts, const bool debug=false)
cut edges
MoFEMErrorCode moveMidNodesOnCutEdges(Tag th=NULL)
projecting of mid edge nodes on new mesh on surface
MoFEMErrorCode saveCutEdges(const std::string prefix="")

◆ cutTrimAndMerge()

MoFEMErrorCode MoFEM::CutMeshInterface::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.

Parameters
first_bitfirst bit of bit revel, subsequent set bits are for trim and merge
fraction_levelfraction of edges merged at each merge step
thtag storring mesh node positions
tol_geometrytolerance how mesh node should be close to cut surface (mesh node is moved), should be small
tol_cut_closehow crack node should be close to mesh (cut surface node is moved), can be big
tol_trim_closehow front node should be close to mesh, can be big
fixed_edgesedges on which nodes can not be moved
corner_nodesnodes which can not be moved
update_meshsetsupdate meshsets by parents
debugswich on debugging
Returns
MoFEMErrorCode
Examples
mesh_cut.cpp.

Definition at line 379 of file CutMeshInterface.cpp.

382 {
383 CoreInterface &m_field = cOre;
385
386 std::vector<BitRefLevel> bit_levels;
387
388 auto get_back_bit_levels = [&]() {
389 bit_levels.push_back(BitRefLevel());
390 bit_levels.back().set(first_bit + bit_levels.size() - 1);
391 return bit_levels.back();
392 };
393
394 if (debug) {
395 CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
396 "ents_not_in_database.vtk", "VTK", "");
397 }
398
399 CHKERR cutAndTrim(first_bit, th, tol_geometry, tol_cut_close, tol_trim_close,
400 &fixed_edges, &corner_nodes, update_meshsets, debug);
401 if (debug)
402 CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
403 "cut_trim_ents_not_in_database.vtk", "VTK", "");
404
405 BitRefLevel bit_level1 = BitRefLevel().set(first_bit - 1);
406 BitRefLevel bit_level2 = get_back_bit_levels();
407 BitRefLevel bit_level3 = get_back_bit_levels();
408
409 CHKERR mergeBadEdges(fraction_level, bit_level1, bit_level2, bit_level3,
410 getNewTrimSurfaces(), fixed_edges, corner_nodes, th,
411 update_meshsets, debug);
412
413 auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
414 Range tets_level; // test at level
415 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
416 bit, BitRefLevel().set(), MBTET, tets_level);
417 double min_q = 1;
418 CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
419 if (min_q < 0 && debug) {
420 CHKERR m_field.getInterface<Tools>()->writeTetsWithQuality(
421 "negative_tets.vtk", "VTK", "", tets_level, th);
422 }
423 return min_q;
424 };
425
426 MOFEM_LOG_C("WORLD", Sev::inform, "Min quality node merge %6.4g",
427 get_min_quality(bit_level3, th));
428
429 CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(constrainSurface,
431 CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(fixed_edges,
432 fixed_edges);
433 CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(corner_nodes,
434 corner_nodes);
435
436 first_bit += bit_levels.size() - 1;
437
438 if (debug) {
439 CHKERR cOre.getInterface<BitRefManager>()->writeBitLevelByType(
440 bit_level3, BitRefLevel().set(), MBTET, "out_tets_merged.vtk", "VTK",
441 "");
442 CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
443 "cut_trim_merge_ents_not_in_database.vtk", "VTK", "");
444 CHKERR SaveData(m_field.get_moab())("merged_surface.vtk", mergedSurfaces);
445 }
446
448}
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 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

◆ findEdgesToCut()

MoFEMErrorCode MoFEM::CutMeshInterface::findEdgesToCut ( Range  vol,
int  verb = QUIET,
const bool  debug = false 
)

find edges to cut

Parameters
volis tetrahedrons search to cut
verbverbosity level
debugdebugging
Returns
MoFEMErrorCode

Definition at line 897 of file CutMeshInterface.cpp.

898 {
899 CoreInterface &m_field = cOre;
900 moab::Interface &moab = m_field.get_moab();
902
903 edgesToCut.clear();
904 cutEdges.clear();
905
906 Tag th_signed_dist;
907 CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_SIGNED", th_signed_dist);
908
909 auto get_tag_edge_dist = [&](auto th, auto conn) {
910 std::array<double, 2> r;
911 CHKERR moab.tag_get_data(th, conn, 2, r.data());
912 return r;
913 };
914
915 auto get_conn = [&](const auto e) {
916 int num_nodes;
917 const EntityHandle *conn;
918 CHKERR moab.get_connectivity(e, conn, num_nodes, true);
919 return conn;
920 };
921
922 auto get_adj = [&moab](const Range r, const int dim) {
923 Range a;
924 if (dim)
925 CHKERR moab.get_adjacencies(r, dim, false, a, moab::Interface::UNION);
926 else
927 CHKERR moab.get_connectivity(r, a, true);
928 return a;
929 };
930
931 auto vol_edges = get_adj(vol, 1);
932
933 aveLength = 0;
934 maxLength = 0;
935 int nb_ave_length = 0;
936 for (auto e : vol_edges) {
937
938 auto conn = get_conn(e);
939 auto dist = get_tag_edge_dist(th_signed_dist, conn);
940 const auto dist_max = std::max(dist[0], dist[1]);
941 const auto dist_min = std::min(dist[0], dist[1]);
942 const auto dot = dist[0] * dist[1];
943
944 VectorDouble6 coords(6);
945 CHKERR moab.get_coords(conn, 2, &coords[0]);
946 VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
947 VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
948 VectorDouble3 ray = n1 - n0;
949 const double ray_length = norm_2(ray);
950 ray /= ray_length;
951
952 if (
953
954 (dot < 0 && dist_max > std::numeric_limits<double>::epsilon()) ||
955 (std::abs(dist_min) < std::numeric_limits<double>::epsilon() &&
956 dist_max > std::numeric_limits<double>::epsilon())
957
958 ) {
959
960 // Edges is on two sides of the surface
961 const double s = dist[0] / (dist[0] - dist[1]);
962 const double dist = s * ray_length;
963
964 auto add_edge = [&](auto dist) {
965 edgesToCut[e].dIst = dist;
966 edgesToCut[e].lEngth = ray_length;
967 edgesToCut[e].unitRayDir = ray;
968 edgesToCut[e].rayPoint = n0;
969 cutEdges.insert(e);
970
971 aveLength += norm_2(ray);
972 maxLength = fmax(maxLength, norm_2(ray));
973 ++nb_ave_length;
974 };
975
976 add_edge(dist);
977 }
978 }
979 aveLength /= nb_ave_length;
980
981 if (debug)
982 CHKERR SaveData(m_field.get_moab())("cut_edges.vtk", cutEdges);
983
984 if (debug)
985 CHKERR SaveData(m_field.get_moab())("cut_edges_vol.vtk",
986 get_adj(cutEdges, 3));
987
989}
constexpr double a
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:95
int r
Definition: sdf.py:8
double aveLength
Average edge length.
double maxLength
Maximal edge length.

◆ findEdgesToTrim()

MoFEMErrorCode MoFEM::CutMeshInterface::findEdgesToTrim ( Range fixed_edges,
Range corner_nodes,
Tag  th = NULL,
const double  tol = 1e-4,
const bool  debug = false 
)

Find edges to trimEdges.

To make this work, you need to find edges to cut (findEdgesToCut), then cut edges in the middle (cutEdgesInMiddle) and finally project edges on the surface (moveMidNodesOnCutEdges)

Parameters
verbverbosity level
Returns
error code

Project front entities on on the cut surface plane

Definition at line 1480 of file CutMeshInterface.cpp.

1483 {
1484 CoreInterface &m_field = cOre;
1485 moab::Interface &moab = m_field.get_moab();
1487
1488 // takes body skin
1489 Skinner skin(&moab);
1490 Range tets_skin;
1491 CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
1492 tets_skin.merge(constrainSurface);
1493
1494 // vertices on the skin
1495 Range tets_skin_verts;
1496 CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1497 // edges on the skin
1498 Range tets_skin_edges;
1499 CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1500 moab::Interface::UNION);
1501
1502 // get edges on new surface
1503 Range cut_surface_edges;
1504 CHKERR moab.get_adjacencies(cutNewSurfaces, 1, false, cut_surface_edges,
1505 moab::Interface::UNION);
1506 Range cut_surface_verts;
1507 CHKERR moab.get_connectivity(cutNewSurfaces, cut_surface_verts, true);
1508
1509 Range corners;
1510 if (corner_nodes)
1511 corners = *corner_nodes;
1512
1513 Range fix_edges;
1514 if (fixed_edges)
1515 fix_edges = *fixed_edges;
1516
1517 Range fixed_edges_verts;
1518 if (fixed_edges)
1519 CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts, true);
1520
1521 Range surface_skin;
1522 if (fRont.empty())
1523 CHKERR skin.find_skin(0, sUrface, false, surface_skin);
1524 else
1525 surface_skin = fRont;
1526
1527 auto get_point_coords = [&](EntityHandle v) {
1528 VectorDouble3 point(3);
1529 if (th)
1530 CHKERR moab.tag_get_data(th, &v, 1, &point[0]);
1531 else
1532 CHKERR moab.get_coords(&v, 1, &point[0]);
1533 return point;
1534 };
1535
1536 struct VertMap {
1537 const double d;
1538 const EntityHandle v;
1539 const EntityHandle e;
1540 VertMap(const double d, const EntityHandle v, const EntityHandle e)
1541 : d(d), v(v), e(e) {}
1542 };
1543
1544 typedef multi_index_container<
1545 VertMap,
1546 indexed_by<
1547 ordered_non_unique<member<VertMap, const double, &VertMap::d>>,
1548 ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::v>>,
1549 ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::e>>
1550
1551 >>
1552 VertMapMultiIndex;
1553
1554 VertMapMultiIndex verts_map;
1555
1556 auto add_vert = [&](const EntityHandle v, const EntityHandle e,
1557 const double dist) {
1558 verts_map.insert(VertMap(dist, v, e));
1559 };
1560
1561 // clear data ranges
1562 trimEdges.clear();
1563 edgesToTrim.clear();
1564 verticesOnTrimEdges.clear();
1565 trimNewVertices.clear();
1566
1567 auto cut_this_edge = [&](const EntityHandle e, const double length, auto &ray,
1568 auto &ray_point) {
1569 trimEdges.insert(e);
1570 edgesToTrim[e].dIst = 1;
1571 edgesToTrim[e].lEngth = length;
1572 edgesToTrim[e].unitRayDir.resize(3, false);
1573 edgesToTrim[e].rayPoint.resize(3, false);
1574 for (int ii = 0; ii != 3; ++ii)
1575 edgesToTrim[e].unitRayDir[ii] = ray(ii);
1576 for (int ii = 0; ii != 3; ++ii)
1577 edgesToTrim[e].rayPoint[ii] = ray_point(ii);
1578 };
1579
1583
1584 int num_nodes;
1585
1586 auto get_tag_data = [&](auto th, auto conn) {
1588 CHKERR moab.tag_get_data(th, &conn, 1, &t(0));
1589 return t;
1590 };
1591
1592 double max_edge_length = 0;
1593
1594 /// Project front entities on on the cut surface plane
1595 if (!fRont.empty()) {
1596 // Calculate distances
1597 treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
1598 new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
1600
1601 for (auto s : surface_skin) {
1602
1603 auto project_on_surface = [&](auto &t) {
1605
1606 EntityHandle facets_out;
1607 CHKERR treeSurfPtr->closest_to_location(&t(0), rootSetSurf, &t_p(0),
1608 facets_out);
1609
1611 CHKERR m_field.getInterface<Tools>()->getTriNormal(facets_out,
1612 &t_normal(0));
1613 t_normal(i) /= sqrt(t_normal(i) * t_normal(i));
1614 const double dot = t_normal(i) * (t_p(i) - t(i));
1615 t_p(i) = t(i) + dot * t_normal(i);
1616
1617 return t_p;
1618 };
1619
1620 const EntityHandle *conn;
1621 CHKERR moab.get_connectivity(s, conn, num_nodes, true);
1622
1623 VectorDouble6 coords_front(6);
1624 CHKERR moab.get_coords(conn, num_nodes, &coords_front[0]);
1625
1626 FTensor::Tensor1<double *, 3> t_f0(&coords_front[0], &coords_front[1],
1627 &coords_front[2]);
1628 FTensor::Tensor1<double *, 3> t_f1(&coords_front[3], &coords_front[4],
1629 &coords_front[5]);
1630
1631 auto t_p_f0 = project_on_surface(t_f0);
1632 auto t_p_f1 = project_on_surface(t_f1);
1633
1634 CHKERR moab.set_coords(&conn[0], 1, &t_p_f0(0));
1635 CHKERR moab.set_coords(&conn[1], 1, &t_p_f1(0));
1636 }
1637 }
1638
1639 if (debug)
1640 CHKERR SaveData(moab)("surface_skin_to_trim.vtk", surface_skin);
1641
1643 Tag th_dist_front_vec;
1644 CHKERR moab.tag_get_handle("DIST_FRONT_VECTOR", th_dist_front_vec);
1645
1646 // iterate over edges on cut surface
1647 for (auto e : cut_surface_edges) {
1648
1649 auto get_conn_edge = [&moab](auto e) {
1650 // Get edge connectivity and coords
1651 const EntityHandle *conn_edge;
1652 int num_nodes;
1653 CHKERR moab.get_connectivity(e, conn_edge, num_nodes, true);
1654 return conn_edge;
1655 };
1656
1657 auto get_coords_edge = [&moab](auto conn_edge) {
1658 std::array<double, 6> coords_edge;
1659 CHKERR moab.get_coords(conn_edge, 2, coords_edge.data());
1660 return coords_edge;
1661 };
1662
1663 auto get_ftensor_coords = [](const double *ptr) {
1664 return FTensor::Tensor1<double, 3>{ptr[0], ptr[1], ptr[2]};
1665 };
1666
1667 auto conn_edge = get_conn_edge(e);
1668 auto coords_edge = get_coords_edge(conn_edge);
1669
1670 auto t_dist0 = get_tag_data(th_dist_front_vec, conn_edge[0]);
1671 auto t_dist1 = get_tag_data(th_dist_front_vec, conn_edge[1]);
1672 auto t_e0 = get_ftensor_coords(&coords_edge[0]);
1673 auto t_e1 = get_ftensor_coords(&coords_edge[3]);
1674
1675 FTensor::Tensor1<double, 3> t_edge_delta;
1676 t_edge_delta(i) = t_e1(i) - t_e0(i);
1677 const double edge_length2 = t_edge_delta(i) * t_edge_delta(i);
1678 const double edge_length = sqrt(edge_length2);
1679 if (edge_length == 0)
1680 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Zero edge length");
1681
1682 auto add_edge = [&](auto t_cut) {
1684 if (t_cut < 0.5) {
1685 t_ray(i) = t_cut * t_edge_delta(i);
1686 add_vert(conn_edge[0], e, fabs(t_cut));
1687 add_vert(conn_edge[1], e, fabs(t_cut - 1));
1688 cut_this_edge(e, edge_length, t_ray, t_e0);
1689 } else {
1690 FTensor::Tensor1<double, 3> t_edge_point;
1691 t_edge_point(i) = t_e0(i) + t_cut * t_edge_delta(i);
1692 t_ray(i) = t_edge_point(i) - t_e1(i);
1693 add_vert(conn_edge[0], e, fabs(t_cut));
1694 add_vert(conn_edge[1], e, fabs(t_cut - 1));
1695 cut_this_edge(e, edge_length, t_ray, t_e1);
1696 }
1697 max_edge_length = std::max(max_edge_length, edge_length);
1698 };
1699
1700 auto trim_cross_edges = [&]() {
1701 if (std::min(t_dist0(i) * t_dist0(i), t_dist1(i) * t_dist1(i)) <
1702 edge_length2) {
1703
1704 auto make_pair = [&](const double d, const double te) {
1705 return std::make_pair(d, te);
1706 };
1707
1708 auto min_pair = make_pair(-1, -1);
1709
1710 for (auto f : surface_skin) {
1711
1712 auto conn_front = get_conn_edge(f);
1713 auto coords_front = get_coords_edge(conn_front);
1714 auto t_f0 = get_ftensor_coords(&coords_front[0]);
1715 auto t_f1 = get_ftensor_coords(&coords_front[3]);
1716
1717 double te;
1718 double tf;
1719
1721 &t_e0(0), &t_e1(0), &t_f0(0), &t_f1(0), &te, &tf);
1722
1723 if (res != Tools::NO_SOLUTION) {
1724
1725 auto check = [](auto v) {
1726 return (v > -std::numeric_limits<double>::epsilon() &&
1727 (v - 1) < std::numeric_limits<double>::epsilon());
1728 };
1729
1730 if (check(te) && check(tf)) {
1731
1732 FTensor::Tensor1<double, 3> t_delta, t_cross;
1733 t_delta(i) = (t_e0(i) + te * t_edge_delta(i)) -
1734 (t_f0(0) + tf * (t_f1(i) - t_f0(i)));
1735
1736 t_delta(i) /= edge_length;
1737
1738 const double dot = t_delta(i) * t_delta(i);
1739
1740 if (min_pair.first < 0 || min_pair.first > dot) {
1741
1742 if (dot < tol_trim_close * tol_trim_close)
1743 min_pair = make_pair(dot, te);
1744 }
1745 }
1746 }
1747 }
1748
1749 if (min_pair.first > -std::numeric_limits<double>::epsilon()) {
1750 add_edge(min_pair.second);
1751 return true;
1752 }
1753 }
1754
1755 return false;
1756 };
1757
1758 if (!trim_cross_edges()) {
1759
1760 const double dot = t_dist0(i) * t_dist1(i);
1761 const double dot_direction0 = t_dist0(i) * t_edge_delta(i);
1762 const double dot_direction1 = t_dist1(i) * t_edge_delta(i);
1763
1764 if (dot < std::numeric_limits<double>::epsilon() &&
1765 dot_direction0 > -std::numeric_limits<double>::epsilon() &&
1766 dot_direction1 < std::numeric_limits<double>::epsilon()) {
1767
1768 const double s0 = t_dist0(i) * t_edge_delta(i) / edge_length;
1769 const double s1 = t_dist1(i) * t_edge_delta(i) / edge_length;
1770 const double t_cut = s0 / (s0 - s1);
1771
1772 add_edge(t_cut);
1773 }
1774 }
1775 }
1776
1777 if (debug)
1778 CHKERR SaveData(moab)("edges_potentially_to_trim.vtk", cut_surface_edges);
1779
1780 if (debug)
1781 CHKERR SaveData(moab)("edges_selected_to_trim.vtk", trimEdges);
1782
1783 auto get_quality_change = [&](const Range &adj_tets, const EntityHandle &v,
1784 const VectorDouble3 &new_pos) {
1785 double q0 = 1;
1786 CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0, th);
1787
1788 double min_q = 1;
1789 for (auto t : adj_tets) {
1790 int num_nodes;
1791 const EntityHandle *conn;
1792 CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1793 VectorDouble12 coords(12);
1794 if (th)
1795 CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1796 else
1797 CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1798
1799 for (int n = 0; n != 4; ++n) {
1800 auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1801 if (conn[n] == v) {
1802 noalias(n_coords) = new_pos;
1803 } else {
1804 auto m = verticesOnTrimEdges.find(conn[n]);
1805 if (m != verticesOnTrimEdges.end())
1806 noalias(n_coords) =
1807 m->second.rayPoint + m->second.dIst * m->second.unitRayDir;
1808 }
1809 }
1810
1811 const double q = Tools::volumeLengthQuality(&coords[0]);
1812 if (!std::isnormal(q)) {
1813 min_q = -2;
1814 break;
1815 }
1816 min_q = std::min(min_q, q);
1817 }
1818 return min_q / q0;
1819 };
1820
1821 Range checked_verts;
1822 auto add_trim_vert = [&](const EntityHandle v, const EntityHandle e) {
1823 if (!(trimNewVertices.find(v) != trimNewVertices.end())) {
1824 auto &r = edgesToTrim.at(e);
1825 VectorDouble3 ray_point = get_point_coords(v);
1826 VectorDouble3 new_pos = r.rayPoint + r.dIst * r.unitRayDir;
1827 VectorDouble3 unit_ray_dir = ray_point - new_pos;
1828 Range adj_tets;
1829 CHKERR moab.get_adjacencies(&v, 1, 3, false, adj_tets);
1830 adj_tets = intersect(adj_tets, cutNewVolumes);
1831 double q = get_quality_change(adj_tets, v, new_pos);
1833 VectorDouble3 ray_dir = new_pos - ray_point;
1834 double dist = norm_2(ray_dir);
1835 verticesOnTrimEdges[v].dIst = 1;
1836 verticesOnTrimEdges[v].lEngth = dist;
1837 verticesOnTrimEdges[v].unitRayDir = ray_dir;
1838 verticesOnTrimEdges[v].rayPoint = ray_point;
1839 trimNewVertices.insert(v);
1840 }
1841 checked_verts.insert(v);
1842 }
1843 };
1844
1845 auto add_no_move_trim = [&](const EntityHandle v, const EntityHandle e) {
1846 if (!(trimNewVertices.find(v) != trimNewVertices.end())) {
1847 auto &m = edgesToTrim.at(e);
1849 verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
1850 verticesOnTrimEdges[v].lEngth = 0;
1851 verticesOnTrimEdges[v].dIst = 0;
1852 trimNewVertices.insert(v);
1853 checked_verts.insert(v);
1854 }
1855 };
1856
1857 // Iterate over all vertives close to surface front and check if those can
1858 // be moved
1859
1860 // filer nodes which distance is in given tolerance
1861 auto hi = verts_map.get<0>().upper_bound(tol_trim_close);
1862 verts_map.get<0>().erase(hi, verts_map.end());
1863
1864 auto remove_verts = [&](Range nodes) {
1865 for (auto n : nodes) {
1866 auto r = verts_map.get<1>().equal_range(n);
1867 verts_map.get<1>().erase(r.first, r.second);
1868 }
1869 };
1870
1871 auto trim_verts = [&](const Range verts, const bool move) {
1872 VertMapMultiIndex verts_map_tmp;
1873 for (auto p = corners.pair_begin(); p != corners.pair_end(); ++p) {
1874 auto lo = verts_map.get<1>().lower_bound(p->first);
1875 auto hi = verts_map.get<1>().upper_bound(p->second);
1876 verts_map_tmp.insert(lo, hi);
1877 }
1878 if (move) {
1879 for (auto &m : verts_map_tmp.get<0>())
1880 add_trim_vert(m.v, m.e);
1881 } else {
1882 for (auto &m : verts_map_tmp.get<0>())
1883 add_no_move_trim(m.v, m.e);
1884 }
1885 };
1886
1887 auto trim_edges = [&](const Range ents, const bool move) {
1888 VertMapMultiIndex verts_map_tmp;
1889 for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
1890 auto lo = verts_map.get<2>().lower_bound(p->first);
1891 auto hi = verts_map.get<2>().upper_bound(p->second);
1892 verts_map_tmp.insert(lo, hi);
1893 verts_map.get<2>().erase(lo, hi);
1894 }
1895 if (move) {
1896 for (auto &m : verts_map_tmp.get<0>())
1897 add_trim_vert(m.v, m.e);
1898 } else {
1899 for (auto &m : verts_map_tmp.get<0>())
1900 add_no_move_trim(m.v, m.e);
1901 }
1902 };
1903
1904 auto intersect_self = [&](Range &a, const Range b) { a = intersect(a, b); };
1905
1906 map<std::string, Range> range_maps;
1907
1908 // Edges on surface skin
1909 CHKERR skin.find_skin(0, cutNewSurfaces, false, range_maps["surface_skin"]);
1910 // Edges as trimmed
1911 intersect_self(range_maps["surface_skin"], trimEdges);
1912 // Edges are on fixed edges
1913 range_maps["fixed_edges_on_surface_skin"] =
1914 intersect(range_maps["surface_skin"], fix_edges);
1915
1916 // Edges adjacent to fixed edges
1917 CHKERR moab.get_adjacencies(fixed_edges_verts, 1, false,
1918 range_maps["fixed_edges_verts_edges"],
1919 moab::Interface::UNION);
1920 intersect_self(range_maps["fixed_edges_verts_edges"], trimEdges);
1921
1922 CHKERR moab.get_connectivity(
1923 range_maps["fixed_edges_verts_edges"],
1924 range_maps["fixed_edges_verts_edges_verts_on_the_skin"], false);
1925 intersect_self(range_maps["fixed_edges_verts_edges_verts_on_the_skin"],
1926 tets_skin_verts);
1927
1928 // Do not move nodes at the corners
1929 trim_verts(corners, false);
1930 remove_verts(corners);
1931
1932 // Trim edges on the body skin move
1933 trim_edges(range_maps["fixed_edges_on_surface_skin"], true);
1934 remove_verts(range_maps["fixed_edges_on_surface_skin_verts"]);
1935
1936 // Trim edges nodes on the skin but edge not on the skin do not
1937 trim_verts(range_maps["fixed_edges_verts_edges_verts_on_the_skin"], false);
1938 remove_verts(range_maps["fixed_edges_verts_edges_verts_on_the_skin"]);
1939
1940 // Trim edges on the cut skin move
1941 trim_edges(range_maps["surface_skin"], true);
1942 trim_verts(tets_skin_verts, false);
1943 remove_verts(tets_skin_verts);
1944
1945 for (auto &m : verts_map.get<0>())
1946 add_trim_vert(m.v, m.e);
1947
1948 for (auto m : verticesOnTrimEdges) {
1949 EntityHandle v = m.first;
1950 Range adj;
1951 CHKERR moab.get_adjacencies(&v, 1, 1, false, adj);
1952 for (auto e : adj) {
1953 edgesToTrim.erase(e);
1954 trimEdges.erase(e);
1955 }
1956 }
1957
1958 if (debug)
1959 if (!trimNewVertices.empty())
1960 CHKERR SaveData(moab)("trim_close_vertices.vtk", trimNewVertices);
1961
1962 if (debug)
1963 if (!trimEdges.empty())
1964 CHKERR SaveData(moab)("trim_edges.vtk", trimEdges);
1965
1967}
static Index< 'p', 3 > p
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Types.hpp:97
map< EntityHandle, TreeData > edgesToTrim
map< EntityHandle, TreeData > verticesOnTrimEdges
MoFEMErrorCode createFrontLevelSets(Range vol, Tag th=nullptr, int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to surface front.
static SEGMENT_MIN_DISTANCE minDistanceFromSegments(const double *w_ptr, const double *v_ptr, const double *k_ptr, const double *l_ptr, double *const tvw_ptr=nullptr, double *const tlk_ptr=nullptr)
Find points on two segments in closest distance.
Definition: Tools.cpp:436
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:17

◆ findLevelSetVolumes() [1/2]

MoFEMErrorCode MoFEM::CutMeshInterface::findLevelSetVolumes ( int  verb = QUIET,
const bool  debug = false 
)

Find level set volumes.

Parameters
update_front
verb
debug
Returns
MoFEMErrorCode

Definition at line 778 of file CutMeshInterface.cpp.

779 {
780 CoreInterface &m_field = cOre;
781 moab::Interface &moab = m_field.get_moab();
783
784 CHKERR createFrontLevelSets(vOlume, nullptr, verb, debug);
785 Tag th_dist_front_vec;
786 CHKERR moab.tag_get_handle("DIST_FRONT_VECTOR", th_dist_front_vec);
787 CHKERR findLevelSetVolumes(th_dist_front_vec, cutFrontVolumes, true, verb,
788 debug, "cutFrontEdges.vtk");
789
791
792 Tag th_dist_surface_vec;
793 CHKERR moab.tag_get_handle("DIST_SURFACE_VECTOR", th_dist_surface_vec);
794 cutSurfaceVolumes.clear();
795 CHKERR findLevelSetVolumes(th_dist_surface_vec, cutSurfaceVolumes, true, verb,
796 debug, "cutSurfaceEdges.vtk");
797
798 if (debug)
799 CHKERR SaveData(m_field.get_moab())("level_sets.vtk", vOlume);
800 if (debug)
801 CHKERR SaveData(m_field.get_moab())("cutSurfaceVolumes.vtk",
803 if (debug)
804 CHKERR SaveData(m_field.get_moab())("cutFrontVolumes.vtk", cutFrontVolumes);
805
807}
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 createSurfaceLevelSets(int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to cut surface.

◆ findLevelSetVolumes() [2/2]

MoFEMErrorCode MoFEM::CutMeshInterface::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.

Parameters
th
vol_edges
remove_adj_prims_edges
verb
debug
edges_file_name
Returns
MoFEMErrorCode

Definition at line 662 of file CutMeshInterface.cpp.

664 {
665 CoreInterface &m_field = cOre;
666 moab::Interface &moab = m_field.get_moab();
668
669 auto get_tag_dim = [&](auto th) {
670 int dim;
671 CHKERR moab.tag_get_length(th, dim);
672 return dim;
673 };
674 auto dim = get_tag_dim(th);
675
676 auto get_tag_data = [&](auto th, auto conn) {
677 const void *ptr;
678 CHKERR moab.tag_get_by_ptr(th, &conn, 1, &ptr);
679 return getVectorAdaptor(
680 const_cast<double *>(static_cast<const double *>(ptr)), 3);
681 };
682
683 auto get_edge_ray = [&](auto conn) {
684 VectorDouble6 coords(6);
685 CHKERR moab.get_coords(conn, 2, &coords[0]);
686 VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
687 VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
688 VectorDouble3 ray = n1 - n0;
689 return ray;
690 };
691
692 Range edges;
693 CHKERR moab.get_adjacencies(vOlume, 1, true, edges, moab::Interface::UNION);
694
695 auto remove_prisms_edges = [&](Range &edges) {
697 Range prisms;
698 CHKERR moab.get_adjacencies(edges, 3, false, prisms,
699 moab::Interface::UNION);
700 prisms = prisms.subset_by_type(MBPRISM);
701 Range prisms_verts;
702 CHKERR moab.get_connectivity(prisms, prisms_verts, true);
703 Range prism_edges;
704 CHKERR moab.get_adjacencies(prisms_verts, 1, false, prism_edges,
705 moab::Interface::UNION);
706 edges = subtract(edges, prism_edges);
708 };
709 if (remove_adj_prims_edges)
710 CHKERR remove_prisms_edges(edges);
711
712 std::vector<EntityHandle> cut_edges;
713 cut_edges.reserve(edges.size());
714
715 auto get_cut_edges_vec = [&](auto th, auto &cut_edges, auto e, auto &&conn) {
717
718 auto ray = get_edge_ray(conn);
719 const double length = norm_2(ray);
720 ray /= length;
721 const auto dist0 = get_tag_data(th, conn[0]);
722 const auto dist1 = get_tag_data(th, conn[1]);
723 const double max_dist = std::max(norm_2(dist0), norm_2(dist1));
724 if (max_dist < length) {
725 cut_edges.push_back(e);
726 }
727
729 };
730
731 auto get_cut_edges_signed_dist = [&](auto th, auto &cut_edges, auto e,
732 auto &&conn) {
734 auto get_tag_signed_dist = [&](auto conn) {
735 double dist;
736 CHKERR moab.tag_get_data(th, &conn, 1, &dist);
737 return dist;
738 };
739 const auto dist0 = get_tag_signed_dist(conn[0]);
740 const auto dist1 = get_tag_signed_dist(conn[1]);
741 const double opposite = dist0 * dist1;
742 if (opposite <= 0)
743 cut_edges.push_back(e);
745 };
746
747 auto get_conn = [&](auto e) {
748 int num_nodes;
749 const EntityHandle *conn;
750 CHKERR moab.get_connectivity(e, conn, num_nodes, true);
751 return conn;
752 };
753
754 if (dim == 3)
755 for (auto e : edges)
756 CHKERR get_cut_edges_vec(th, cut_edges, e, get_conn(e));
757 else
758 for (auto e : edges)
759 CHKERR get_cut_edges_signed_dist(th, cut_edges, e, get_conn(e));
760
761 CHKERR moab.get_adjacencies(&*cut_edges.begin(), cut_edges.size(), 3, false,
762 vol_edges, moab::Interface::UNION);
763 vol_edges = intersect(vol_edges, vOlume);
764
765 auto convert_to_ranger = [](auto &v) {
766 Range r;
767 r.insert_list(v.begin(), v.end());
768 return r;
769 };
770
771 if (debug && !edges_file_name.empty())
772 CHKERR SaveData(m_field.get_moab())(edges_file_name,
773 convert_to_ranger(cut_edges));
774
776}

◆ getCutEdges()

const Range & MoFEM::CutMeshInterface::getCutEdges ( ) const
inline

Definition at line 469 of file CutMeshInterface.hpp.

469{ return cutEdges; }

◆ getCutFrontVolumes()

const Range & MoFEM::CutMeshInterface::getCutFrontVolumes ( ) const
inline

Definition at line 486 of file CutMeshInterface.hpp.

486{ return cutFrontVolumes; }

◆ getCutSurfaceVolumes()

const Range & MoFEM::CutMeshInterface::getCutSurfaceVolumes ( ) const
inline

Definition at line 485 of file CutMeshInterface.hpp.

485{ return cutSurfaceVolumes; }

◆ getFront()

const Range & MoFEM::CutMeshInterface::getFront ( ) const
inline

Definition at line 467 of file CutMeshInterface.hpp.

467{ return fRont; }

◆ getMergedSurfaces()

const Range & MoFEM::CutMeshInterface::getMergedSurfaces ( ) const
inline
Examples
mesh_cut.cpp.

Definition at line 483 of file CutMeshInterface.hpp.

483{ return mergedSurfaces; }

◆ getMergedVolumes()

const Range & MoFEM::CutMeshInterface::getMergedVolumes ( ) const
inline

Definition at line 482 of file CutMeshInterface.hpp.

◆ getNewCutSurfaces()

const Range & MoFEM::CutMeshInterface::getNewCutSurfaces ( ) const
inline

Definition at line 471 of file CutMeshInterface.hpp.

471{ return cutNewSurfaces; }

◆ getNewCutVertices()

const Range & MoFEM::CutMeshInterface::getNewCutVertices ( ) const
inline

Definition at line 472 of file CutMeshInterface.hpp.

472{ return cutNewVertices; }

◆ getNewCutVolumes()

const Range & MoFEM::CutMeshInterface::getNewCutVolumes ( ) const
inline

Definition at line 470 of file CutMeshInterface.hpp.

470{ return cutNewVolumes; }

◆ getNewTrimSurfaces()

const Range & MoFEM::CutMeshInterface::getNewTrimSurfaces ( ) const
inline

Definition at line 479 of file CutMeshInterface.hpp.

◆ getNewTrimVertices()

const Range & MoFEM::CutMeshInterface::getNewTrimVertices ( ) const
inline

Definition at line 480 of file CutMeshInterface.hpp.

480{ return trimNewVertices; }

◆ getNewTrimVolumes()

const Range & MoFEM::CutMeshInterface::getNewTrimVolumes ( ) const
inline

Definition at line 478 of file CutMeshInterface.hpp.

◆ getOptions()

MoFEMErrorCode MoFEM::CutMeshInterface::getOptions ( )
inline

Get options from command line.

Returns
error code

Definition at line 37 of file CutMeshInterface.hpp.

37 {
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 }
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76

◆ getSubInterfaceOptions()

MoFEMErrorCode MoFEM::CutMeshInterface::getSubInterfaceOptions ( )
inline

Definition at line 31 of file CutMeshInterface.hpp.

31{ return getOptions(); }
MoFEMErrorCode getOptions()
Get options from command line.

◆ getSurface()

const Range & MoFEM::CutMeshInterface::getSurface ( ) const
inline
Examples
mesh_cut.cpp.

Definition at line 466 of file CutMeshInterface.hpp.

466{ return sUrface; }

◆ getTreeSurfPtr()

boost::shared_ptr< OrientedBoxTreeTool > & MoFEM::CutMeshInterface::getTreeSurfPtr ( )
inline

Definition at line 492 of file CutMeshInterface.hpp.

492 {
493 return treeSurfPtr;
494 }

◆ getTrimEdges()

const Range & MoFEM::CutMeshInterface::getTrimEdges ( ) const
inline

Definition at line 477 of file CutMeshInterface.hpp.

477{ return trimEdges; }

◆ getVolume()

const Range & MoFEM::CutMeshInterface::getVolume ( ) const
inline

Definition at line 465 of file CutMeshInterface.hpp.

465{ return vOlume; }

◆ makeFront()

MoFEMErrorCode MoFEM::CutMeshInterface::makeFront ( const bool  debug = false)

Create front from the surface.

Parameters
debug
Returns
MoFEMErrorCode
Examples
mesh_cut.cpp.

Definition at line 450 of file CutMeshInterface.cpp.

450 {
451 CoreInterface &m_field = cOre;
452 moab::Interface &moab = m_field.get_moab();
454 Skinner skin(&moab);
455 Range tets_skin;
456 CHKERR skin.find_skin(0, vOlume, false, tets_skin);
457 Range tets_skin_edges;
458 ErrorCode tmp_result;
459 tmp_result = moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
460 moab::Interface::UNION);
461
462 if (MB_SUCCESS != tmp_result)
463 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
464 "Duplicated edges: most likely the source of error is comming from "
465 "adding the vertices of the cracking "
466 "volume to a BLOCKSET rather than NODESET (corresponding to the "
467 "input parameter-vertex_block_set)");
468
469 Range surface_skin;
470 CHKERR skin.find_skin(0, sUrface, false, surface_skin);
471 fRont = subtract(surface_skin, tets_skin_edges);
472 if (debug)
473 CHKERR SaveData(moab)("front.vtk", fRont);
475}

◆ mergeBadEdges() [1/2]

MoFEMErrorCode MoFEM::CutMeshInterface::mergeBadEdges ( const int  fraction_level,
const BitRefLevel  cut_bit,
const BitRefLevel  trim_bit,
const BitRefLevel  bit,
const Range surface,
const Range fixed_edges,
const Range corner_nodes,
Tag  th = NULL,
const bool  update_meshsets = false,
const bool  debug = false 
)

Merge edges.

Sort all edges, where sorting is by quality calculated as edge length times quality of tets adjacent to the edge. Edge is merged if quality if the mesh is improved.

Definition at line 3062 of file CutMeshInterface.cpp.

3066 {
3067 CoreInterface &m_field = cOre;
3069 Range tets_level;
3070 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3071 trim_bit, BitRefLevel().set(), MBTET, tets_level);
3072
3073 Range edges_to_merge;
3074 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
3075 trim_bit, cut_bit | trim_bit, edges_to_merge);
3076
3077 // get all entities not in database
3078 Range all_ents_not_in_database_before;
3079 CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
3080 all_ents_not_in_database_before);
3081
3082 edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
3083 if (debug)
3084 CHKERR SaveData(m_field.get_moab())("edges_to_merge.vtk", edges_to_merge);
3085
3086 Range out_new_tets, out_new_surf;
3087 CHKERR mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
3088 corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
3089 th, update_meshsets, &bit, debug);
3090
3091 // get all entities not in database after merge
3092 Range all_ents_not_in_database_after;
3093 CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
3094 all_ents_not_in_database_after);
3095
3096 // delete hanging entities
3097 all_ents_not_in_database_after =
3098 subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
3099 Range meshsets;
3100 CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
3101 false);
3102 for (auto m : meshsets)
3103 CHKERR m_field.get_moab().remove_entities(m,
3104 all_ents_not_in_database_after);
3105
3106 m_field.get_moab().delete_entities(all_ents_not_in_database_after);
3107
3108 mergedVolumes.swap(out_new_tets);
3109 mergedSurfaces.swap(out_new_surf);
3110
3112}

◆ mergeBadEdges() [2/2]

MoFEMErrorCode MoFEM::CutMeshInterface::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.

Sort all edges, where sorting is by quality calculated as edge length times quality of tets adjacent to the edge. Edge is merged if quality if the mesh is improved.

Parameters
fraction_levelFraction of edges attemt to be merged at iteration
tetsTets of the mesh which edges are merged
surfaceSurface created by edge spliting
fixed_edgesedges which are geometrical corners of the body
corner_nodesvertices on the corners
merged_nodesmerged nodes
out_tetsreturned test after merge
new_surfnew surface without merged edges
thtag with nodal positons
bit_ptrset bit ref level to mesh without merged edges
debug
Returns
MoFEMErrorCode

Merge nodes

Calculate edge length

< Controls importance of length when ranking edges for merge

< Controls importance of quality when ranking edges for merge

< Do not merge quality if quality above accepted thrashold

Topological relations

Definition at line 2302 of file CutMeshInterface.cpp.

2306 {
2307 CoreInterface &m_field = cOre;
2308 moab::Interface &moab = m_field.get_moab();
2310
2311 /**
2312 * \brief Merge nodes
2313 */
2314 struct MergeNodes {
2315 CoreInterface &mField;
2316 const bool onlyIfImproveQuality;
2317 Tag tH;
2318 bool updateMehsets;
2319
2320 MergeNodes(CoreInterface &m_field, const bool only_if_improve_quality,
2321 Tag th, bool update_mehsets)
2322 : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
2323 tH(th), updateMehsets(update_mehsets) {
2324 mField.getInterface(nodeMergerPtr);
2325 }
2326 NodeMergerInterface *nodeMergerPtr;
2327 MoFEMErrorCode mergeNodes(int line_search, EntityHandle father,
2328 EntityHandle mother, Range &proc_tets,
2329 bool add_child = true) {
2330 moab::Interface &moab = mField.get_moab();
2332 const EntityHandle conn[] = {father, mother};
2333 Range vert_tets;
2334 CHKERR moab.get_adjacencies(conn, 2, 3, false, vert_tets,
2335 moab::Interface::UNION);
2336 vert_tets = intersect(vert_tets, proc_tets);
2337 Range out_tets;
2338 CHKERR nodeMergerPtr->mergeNodes(father, mother, out_tets, &vert_tets,
2339 onlyIfImproveQuality, 0, line_search,
2340 tH);
2341
2342 if (add_child && nodeMergerPtr->getSuccessMerge()) {
2343
2344 Range::iterator lo, hi = proc_tets.begin();
2345 for (auto pt = vert_tets.pair_begin(); pt != vert_tets.pair_end();
2346 ++pt) {
2347 lo = proc_tets.lower_bound(hi, proc_tets.end(), pt->first);
2348 if (lo != proc_tets.end()) {
2349 hi = proc_tets.upper_bound(lo, proc_tets.end(), pt->second);
2350 proc_tets.erase(lo, hi);
2351 } else
2352 break;
2353 }
2354 proc_tets.merge(out_tets);
2355
2356 auto &parent_child_map = nodeMergerPtr->getParentChildMap();
2357
2358 struct ChangeChild {
2359 EntityHandle child;
2360 ChangeChild(const EntityHandle child) : child(child) {}
2361 void operator()(NodeMergerInterface::ParentChild &p) {
2362 p.cHild = child;
2363 }
2364 };
2365
2366 std::vector<decltype(parentsChildMap.get<0>().begin())> it_vec;
2367 it_vec.reserve(parentsChildMap.size());
2368
2369 for (auto &p : parent_child_map) {
2370
2371 it_vec.clear();
2372 for (auto it = parentsChildMap.get<0>().equal_range(p.pArent);
2373 it.first != it.second; ++it.first)
2374 it_vec.emplace_back(it.first);
2375
2376 for (auto it = parentsChildMap.get<1>().equal_range(p.pArent);
2377 it.first != it.second; ++it.first)
2378 it_vec.emplace_back(parentsChildMap.project<0>(it.first));
2379
2380 for (auto &it : it_vec)
2381 parentsChildMap.modify(it, ChangeChild(p.cHild));
2382 }
2383
2384 parentsChildMap.insert(parent_child_map.begin(),
2385 parent_child_map.end());
2386 }
2388 }
2389
2390 MoFEMErrorCode updateRangeByChilds(Range &new_surf, Range &edges_to_merge,
2391 Range &not_merged_edges,
2392 bool add_child) {
2393 moab::Interface &moab = mField.get_moab();
2395 if (add_child) {
2396
2397 std::vector<EntityHandle> parents_ents_vec(parentsChildMap.size());
2398 for (auto &it : parentsChildMap)
2399 parents_ents_vec.emplace_back(it.pArent);
2400 Range parent_ents;
2401 parent_ents.insert_list(parents_ents_vec.begin(),
2402 parents_ents_vec.end());
2403
2404 Range surf_parent_ents = intersect(new_surf, parent_ents);
2405 new_surf = subtract(new_surf, surf_parent_ents);
2406 Range child_surf_ents;
2407 CHKERR updateRangeByChilds(parentsChildMap, surf_parent_ents,
2408 child_surf_ents);
2409 new_surf.merge(child_surf_ents);
2410
2411 Range edges_to_merge_parent_ents =
2412 intersect(edges_to_merge, parent_ents);
2413 edges_to_merge = subtract(edges_to_merge, edges_to_merge_parent_ents);
2414 Range merged_child_edge_ents;
2415 CHKERR updateRangeByChilds(parentsChildMap, edges_to_merge_parent_ents,
2416 merged_child_edge_ents);
2417
2418 Range not_merged_edges_child_ents =
2419 intersect(not_merged_edges, parent_ents);
2420 not_merged_edges =
2421 subtract(not_merged_edges, not_merged_edges_child_ents);
2422 Range not_merged_child_edge_ents;
2423 CHKERR updateRangeByChilds(parentsChildMap, not_merged_edges_child_ents,
2424 not_merged_child_edge_ents);
2425
2426 merged_child_edge_ents =
2427 subtract(merged_child_edge_ents, not_merged_child_edge_ents);
2428 edges_to_merge.merge(merged_child_edge_ents);
2429 not_merged_edges.merge(not_merged_child_edge_ents);
2430
2431 if (updateMehsets) {
2433 (*mField.getInterface<MeshsetsManager>()), cubit_it)) {
2434 EntityHandle cubit_meshset = cubit_it->meshset;
2435 Range meshset_ents;
2436 CHKERR moab.get_entities_by_handle(cubit_meshset, meshset_ents,
2437 true);
2438 Range child_ents;
2439 CHKERR updateRangeByChilds(parentsChildMap, meshset_ents,
2440 child_ents);
2441 CHKERR moab.add_entities(cubit_meshset, child_ents);
2442 }
2443 }
2444 }
2445
2447 };
2448
2449 private:
2450 NodeMergerInterface::ParentChildMap parentsChildMap;
2451 std::vector<EntityHandle> childsVec;
2452
2453 inline MoFEMErrorCode updateRangeByChilds(
2454 const NodeMergerInterface::ParentChildMap &parent_child_map,
2455 const Range &parents, Range &childs) {
2457 childsVec.clear();
2458 childsVec.reserve(parents.size());
2459 for (auto pit = parents.pair_begin(); pit != parents.pair_end(); pit++) {
2460 auto it = parent_child_map.lower_bound(pit->first);
2461 if (it != parent_child_map.end()) {
2462 for (auto hi_it = parent_child_map.upper_bound(pit->second);
2463 it != hi_it; ++it)
2464 childsVec.emplace_back(it->cHild);
2465 }
2466 }
2467 childs.insert_list(childsVec.begin(), childsVec.end());
2469 }
2470 };
2471
2472 /**
2473 * \brief Calculate edge length
2474 */
2475 struct LengthMap {
2476 Tag tH;
2477 CoreInterface &mField;
2478 moab::Interface &moab;
2479 const double maxLength;
2480 LengthMap(CoreInterface &m_field, Tag th, double max_length)
2481 : tH(th), mField(m_field), moab(m_field.get_moab()),
2482 maxLength(max_length) {
2483 gammaL = 1.;
2484 gammaQ = 3.;
2485 acceptedThrasholdMergeQuality = 0.5;
2486 }
2487
2488 double
2489 gammaL; ///< Controls importance of length when ranking edges for merge
2490 double
2491 gammaQ; ///< Controls importance of quality when ranking edges for merge
2492 double acceptedThrasholdMergeQuality; ///< Do not merge quality if quality
2493 ///< above accepted thrashold
2494
2495 MoFEMErrorCode operator()(const Range &tets, const Range &edges,
2496 LengthMapData_multi_index &length_map,
2497 double &ave) const {
2498 int num_nodes;
2499 const EntityHandle *conn;
2500 std::array<double, 12> coords;
2502 VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
2503 VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
2505
2506 struct NodeQuality {
2507 EntityHandle node;
2508 double quality;
2509 NodeQuality(const EntityHandle node) : node(node), quality(1) {}
2510 };
2511
2512 typedef multi_index_container<
2513 NodeQuality, indexed_by<
2514
2515 sequenced<>,
2516
2517 hashed_non_unique<tag<Ent_mi_tag>,
2518 member<NodeQuality, EntityHandle,
2519 &NodeQuality::node>>
2520
2521 >>
2522 NodeQuality_sequence;
2523
2524 NodeQuality_sequence node_quality_sequence;
2525
2526 Range edges_nodes;
2527 CHKERR moab.get_connectivity(tets, edges_nodes, false);
2528 Range edges_tets;
2529 CHKERR moab.get_adjacencies(edges, 3, false, edges_tets,
2530 moab::Interface::UNION);
2531 edges_tets = intersect(edges_tets, tets);
2532
2533 for (auto node : edges_nodes)
2534 node_quality_sequence.emplace_back(node);
2535
2536 for (auto tet : edges_tets) {
2537
2538 CHKERR moab.get_connectivity(tet, conn, num_nodes, true);
2539 if (tH)
2540 CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2541 else
2542 CHKERR moab.get_coords(conn, num_nodes, coords.data());
2543
2544 const double q = Tools::volumeLengthQuality(coords.data());
2545
2546 for (auto n : {0, 1, 2, 3}) {
2547 auto it = node_quality_sequence.get<1>().find(conn[n]);
2548 if (it != node_quality_sequence.get<1>().end())
2549 const_cast<double &>(it->quality) = std::min(q, it->quality);
2550 }
2551 }
2552
2553 for (auto edge : edges) {
2554
2555 CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
2556
2557 if (tH)
2558 CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2559 else
2560 CHKERR moab.get_coords(conn, num_nodes, coords.data());
2561
2562 double q = 1;
2563 for (auto n : {0, 1}) {
2564 auto it = node_quality_sequence.get<1>().find(conn[n]);
2565 if (it != node_quality_sequence.get<1>().end())
2566 q = std::min(q, it->quality);
2567 }
2568
2569 if (q < acceptedThrasholdMergeQuality) {
2570 noalias(delta) = (s0 - s1) / maxLength;
2571 double dot = inner_prod(delta, delta);
2572 double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
2573 length_map.insert(LengthMapData(val, q, edge));
2574 }
2575 }
2576
2577 ave = 0;
2578 for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
2579 length_map.get<0>().begin();
2580 mit != length_map.get<0>().end(); mit++) {
2581 ave += mit->qUality;
2582 }
2583 ave /= length_map.size();
2585 }
2586 };
2587
2588 /**
2589 * \brief Topological relations
2590 */
2591 struct Toplogy {
2592
2593 CoreInterface &mField;
2594 Tag tH;
2595 Toplogy(CoreInterface &m_field, Tag th) : mField(m_field), tH(th) {}
2596
2597 enum TYPE {
2598 FREE = 0,
2599 SKIN = 1 << 0,
2600 SURFACE = 1 << 1,
2601 SURFACE_SKIN = 1 << 2,
2602 FRONT_ENDS = 1 << 3,
2603 FIX_EDGES = 1 << 4,
2604 FIX_CORNERS = 1 << 5
2605 };
2606
2607 typedef map<int, Range> SetsMap;
2608
2609 MoFEMErrorCode classifyVerts(const Range &surface, const Range &tets,
2610 const Range &fixed_edges,
2611 const Range &corner_nodes,
2612 const Range &constrain_surface,
2613 SetsMap &sets_map) const {
2614 moab::Interface &moab(mField.get_moab());
2615 Skinner skin(&moab);
2617
2618 sets_map[FIX_CORNERS].merge(corner_nodes);
2619 Range fixed_verts;
2620 CHKERR moab.get_connectivity(fixed_edges, fixed_verts, true);
2621 sets_map[FIX_EDGES].swap(fixed_verts);
2622
2623 Range tets_skin;
2624 CHKERR skin.find_skin(0, tets, false, tets_skin);
2625 Range tets_skin_edges;
2626 CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2627 moab::Interface::UNION);
2628
2629 // constrain surface
2630 Range constrain_surface_verts;
2631 CHKERR moab.get_connectivity(constrain_surface, constrain_surface_verts,
2632 true);
2633 Range constrain_surface_edges;
2634 CHKERR moab.get_adjacencies(constrain_surface, 1, false,
2635 constrain_surface_edges,
2636 moab::Interface::UNION);
2637
2638 // surface skin
2639 Range surface_skin;
2640 CHKERR skin.find_skin(0, surface, false, surface_skin);
2641 Range front_in_the_body;
2642 front_in_the_body = subtract(surface_skin, tets_skin_edges);
2643 Range front_in_the_body_verts;
2644 CHKERR moab.get_connectivity(front_in_the_body, front_in_the_body_verts,
2645 true);
2646 Range front_ends;
2647 CHKERR skin.find_skin(0, front_in_the_body, false, front_ends);
2648 front_ends.merge(
2649 intersect(front_in_the_body_verts, constrain_surface_verts));
2650 sets_map[FRONT_ENDS].swap(front_ends);
2651
2652 Range surface_skin_verts;
2653 CHKERR moab.get_connectivity(surface_skin, surface_skin_verts, true);
2654 sets_map[SURFACE_SKIN].swap(surface_skin_verts);
2655
2656 // surface
2657 Range surface_verts;
2658 CHKERR moab.get_connectivity(surface, surface_verts, true);
2659 sets_map[SURFACE_SKIN].merge(
2660 intersect(constrain_surface_verts, surface_verts));
2661 sets_map[SURFACE].swap(surface_verts);
2662
2663 // skin
2664 Range tets_skin_verts;
2665 CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
2666 sets_map[SKIN].swap(tets_skin_verts);
2667 sets_map[SKIN].merge(constrain_surface_verts);
2668
2669 Range tets_verts;
2670 CHKERR moab.get_connectivity(tets, tets_verts, true);
2671 sets_map[FREE].swap(tets_verts);
2672
2674 }
2675
2676 MoFEMErrorCode getProcTets(const Range &tets, const Range &edges_to_merge,
2677 Range &proc_tets) const {
2678 moab::Interface &moab(mField.get_moab());
2680 Range edges_to_merge_verts;
2681 CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts, true);
2682 Range edges_to_merge_verts_tets;
2683 CHKERR moab.get_adjacencies(edges_to_merge_verts, 3, false,
2684 edges_to_merge_verts_tets,
2685 moab::Interface::UNION);
2686 edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
2687 proc_tets.swap(edges_to_merge_verts_tets);
2689 }
2690
2691 MoFEMErrorCode removeBadEdges(const Range &surface, const Range &tets,
2692 const Range &fixed_edges,
2693 const Range &corner_nodes,
2694 const Range &constrain_surface,
2695 Range &edges_to_merge,
2696 Range &not_merged_edges) {
2697 moab::Interface &moab(mField.get_moab());
2699
2700 // find skin
2701 Skinner skin(&moab);
2702 Range tets_skin;
2703 CHKERR skin.find_skin(0, tets, false, tets_skin);
2704 Range surface_skin;
2705 CHKERR skin.find_skin(0, surface, false, surface_skin);
2706
2707 // constrain surface
2708 Range constrain_surface_verts;
2709 CHKERR moab.get_connectivity(constrain_surface, constrain_surface_verts,
2710 true);
2711 Range constrain_surface_edges;
2712 CHKERR moab.get_adjacencies(constrain_surface, 1, false,
2713 constrain_surface_edges,
2714 moab::Interface::UNION);
2715
2716 // end nodes
2717 Range tets_skin_edges;
2718 CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2719 moab::Interface::UNION);
2720
2721 Range surface_front;
2722 surface_front = subtract(surface_skin, tets_skin_edges);
2723 Range surface_front_nodes;
2724 CHKERR moab.get_connectivity(surface_front, surface_front_nodes, true);
2725
2726 Range ends_nodes;
2727 CHKERR skin.find_skin(0, surface_front, false, ends_nodes);
2728 ends_nodes.merge(intersect(surface_front_nodes, constrain_surface_verts));
2729
2730 // remove bad merges
2731 surface_skin.merge(constrain_surface);
2732 tets_skin_edges.merge(constrain_surface_edges);
2733
2734 // get surface and body skin verts
2735 Range surface_edges;
2736 CHKERR moab.get_adjacencies(surface, 1, false, surface_edges,
2737 moab::Interface::UNION);
2738 // get nodes on the surface
2739 Range surface_edges_verts;
2740 CHKERR moab.get_connectivity(surface_edges, surface_edges_verts, true);
2741 // get vertices on the body skin
2742 Range tets_skin_edges_verts;
2743 CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts,
2744 true);
2745
2746 Range edges_to_remove;
2747
2748 // remove edges self connected to body skin
2749 {
2750 Range ents_nodes_and_edges;
2751 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2752 ents_nodes_and_edges.merge(tets_skin_edges);
2753 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2754 false);
2755 }
2756 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2757 not_merged_edges.merge(edges_to_remove);
2758
2759 // remove edges self connected to surface
2760 {
2761 Range ents_nodes_and_edges;
2762 ents_nodes_and_edges.merge(surface_edges_verts);
2763 ents_nodes_and_edges.merge(surface_edges);
2764 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2765 ents_nodes_and_edges.merge(tets_skin_edges);
2766 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2767 false);
2768 }
2769 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2770 not_merged_edges.merge(edges_to_remove);
2771
2772 // remove edges adjacent corner_nodes execpt those on fixed edges
2773 Range fixed_edges_nodes;
2774 CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes, true);
2775 {
2776 Range ents_nodes_and_edges;
2777 ents_nodes_and_edges.merge(fixed_edges_nodes);
2778 ents_nodes_and_edges.merge(ends_nodes);
2779 ents_nodes_and_edges.merge(corner_nodes);
2780 ents_nodes_and_edges.merge(fixed_edges);
2781 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2782 false);
2783 }
2784 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2785 not_merged_edges.merge(edges_to_remove);
2786
2787 // remove edges self connected to surface
2788 CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove, false);
2789 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2790 not_merged_edges.merge(edges_to_remove);
2791
2792 // remove edges self contented on surface skin
2793 {
2794 Range ents_nodes_and_edges;
2795 ents_nodes_and_edges.merge(surface_skin);
2796 ents_nodes_and_edges.merge(fixed_edges_nodes);
2797 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2798 false);
2799 }
2800 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2801 not_merged_edges.merge(edges_to_remove);
2802
2803 // remove edges connecting crack front and fixed edges, except those short
2804 {
2805 Range ents_nodes_and_edges;
2806 ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
2807 ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2808 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2809 false);
2810 }
2811 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2812 not_merged_edges.merge(edges_to_remove);
2813
2814 // remove crack front nodes connected to the surface, except those short
2815 {
2816 Range ents_nodes_and_edges;
2817 ents_nodes_and_edges.merge(surface_front_nodes);
2818 ents_nodes_and_edges.merge(surface_front);
2819 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2820 ents_nodes_and_edges.merge(tets_skin_edges);
2821 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2822 false);
2823 }
2824 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2825 not_merged_edges.merge(edges_to_remove);
2826
2828 }
2829
2830 private:
2831 MoFEMErrorCode removeSelfConectingEdges(const Range &ents,
2832 Range &edges_to_remove,
2833 bool debug) const {
2834 moab::Interface &moab(mField.get_moab());
2836 // get nodes
2837 Range ents_nodes = ents.subset_by_type(MBVERTEX);
2838 if (ents_nodes.empty()) {
2839 CHKERR moab.get_connectivity(ents, ents_nodes, true);
2840 }
2841 // edges adj. to nodes
2842 Range ents_nodes_edges;
2843 CHKERR moab.get_adjacencies(ents_nodes, 1, false, ents_nodes_edges,
2844 moab::Interface::UNION);
2845 // nodes of adj. edges
2846 Range ents_nodes_edges_nodes;
2847 CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2848 true);
2849 // hanging nodes
2850 ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2851 Range ents_nodes_edges_nodes_edges;
2852 CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1, false,
2853 ents_nodes_edges_nodes_edges,
2854 moab::Interface::UNION);
2855 // remove edges adj. to hanging edges
2856 ents_nodes_edges =
2857 subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2858 ents_nodes_edges =
2859 subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2860
2861 edges_to_remove.swap(ents_nodes_edges);
2862 if (debug) {
2863 CHKERR SaveData(moab)("edges_to_remove.vtk", edges_to_remove);
2864 }
2866 }
2867 };
2868
2869 Range not_merged_edges;
2870 CHKERR Toplogy(m_field, th)
2871 .removeBadEdges(surface, tets, fixed_edges, corner_nodes,
2872 constrainSurface, edges_to_merge, not_merged_edges);
2873 Toplogy::SetsMap sets_map;
2874 CHKERR Toplogy(m_field, th)
2875 .classifyVerts(surface, tets, fixed_edges, corner_nodes, constrainSurface,
2876 sets_map);
2877 if (debug) {
2878 for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2879 sit != sets_map.rend(); sit++) {
2880 std::string name = "classification_verts_" +
2881 boost::lexical_cast<std::string>(sit->first) + ".vtk";
2882 if (!sit->second.empty())
2883 CHKERR SaveData(moab)(name, sit->second);
2884 }
2885 }
2886 Range proc_tets;
2887 CHKERR Toplogy(m_field, th).getProcTets(tets, edges_to_merge, proc_tets);
2888 out_tets = subtract(tets, proc_tets);
2889
2890 if (bit_ptr) {
2891 Range all_out_ents = out_tets;
2892 for (int dd = 2; dd >= 0; dd--) {
2893 CHKERR moab.get_adjacencies(out_tets, dd, false, all_out_ents,
2894 moab::Interface::UNION);
2895 }
2896 CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(all_out_ents,
2897 *bit_ptr);
2898 }
2899
2900 int nb_nodes_merged = 0;
2901 LengthMapData_multi_index length_map;
2902 new_surf = surface;
2903
2904 auto save_merge_step = [&](const int pp, const Range collapsed_edges) {
2906 Range adj_faces;
2907 CHKERR moab.get_adjacencies(proc_tets, 2, false, adj_faces,
2908 moab::Interface::UNION);
2909 std::string name;
2910 name = "node_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2911 CHKERR SaveData(moab)(
2912 name, unite(intersect(new_surf, adj_faces), collapsed_edges));
2913 name =
2914 "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2915 CHKERR SaveData(moab)(
2916 name, unite(intersect(new_surf, adj_faces), edges_to_merge));
2918 };
2919
2920 if (debug)
2921 CHKERR save_merge_step(0, Range());
2922
2923 double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2924 for (int pp = 0; pp != nbMaxMergingCycles; ++pp) {
2925
2926 int nb_nodes_merged_p = nb_nodes_merged;
2927 length_map.clear();
2928 min_pp = min_p;
2929 min_p = min;
2930 CHKERR LengthMap(m_field, th, aveLength)(proc_tets, edges_to_merge,
2931 length_map, ave);
2932
2933 if(!length_map.empty())
2934 min = length_map.get<2>().begin()->qUality;
2935 if (pp == 0) {
2936 ave0 = ave;
2937 }
2938
2939 int nn = 0;
2940 Range collapsed_edges;
2941 MergeNodes merge_nodes(m_field, true, th, update_meshsets);
2942
2943 for (auto mit = length_map.get<0>().begin();
2944 mit != length_map.get<0>().end(); mit++, nn++) {
2945
2946 if (!mit->skip) {
2947
2948 auto get_father_and_mother =
2949 [&](EntityHandle &father, EntityHandle &mother, int &line_search) {
2951 int num_nodes;
2952 const EntityHandle *conn;
2953 CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes, true);
2954 std::array<int, 2> conn_type = {0, 0};
2955 for (int nn = 0; nn != 2; nn++) {
2956 for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2957 sit != sets_map.rend(); sit++) {
2958 if (sit->second.find(conn[nn]) != sit->second.end()) {
2959 conn_type[nn] |= sit->first;
2960 }
2961 }
2962 }
2963 int type_father, type_mother;
2964 if (conn_type[0] > conn_type[1]) {
2965 father = conn[0];
2966 mother = conn[1];
2967 type_father = conn_type[0];
2968 type_mother = conn_type[1];
2969 } else {
2970 father = conn[1];
2971 mother = conn[0];
2972 type_father = conn_type[1];
2973 type_mother = conn_type[0];
2974 }
2975 if (type_father == type_mother) {
2976 line_search = lineSearchSteps;
2977 }
2979 };
2980
2981 int line_search = 0;
2982 EntityHandle father, mother;
2983 CHKERR get_father_and_mother(father, mother, line_search);
2984 CHKERR merge_nodes.mergeNodes(line_search, father, mother, proc_tets);
2985 if (m_field.getInterface<NodeMergerInterface>()->getSuccessMerge()) {
2986 const EntityHandle father_and_mother[] = {father, mother};
2987 Range adj_tets;
2988 CHKERR moab.get_adjacencies(father_and_mother, 1, 3, false, adj_tets);
2989 Range adj_tets_nodes;
2990 CHKERR moab.get_connectivity(adj_tets, adj_tets_nodes, true);
2991 Range adj_edges;
2992 CHKERR moab.get_adjacencies(adj_tets_nodes, 1, false, adj_edges,
2993 moab::Interface::UNION);
2994 for (auto ait : adj_edges) {
2995 auto miit = length_map.get<1>().find(ait);
2996 if (miit != length_map.get<1>().end())
2997 (const_cast<LengthMapData &>(*miit)).skip = true;
2998 }
2999 nb_nodes_merged++;
3000 collapsed_edges.insert(mit->eDge);
3001 }
3002
3003 if (nn > static_cast<int>(length_map.size() / fraction_level))
3004 break;
3005 if (mit->qUality > ave)
3006 break;
3007 }
3008 }
3009
3010 CHKERR merge_nodes.updateRangeByChilds(new_surf, edges_to_merge,
3011 not_merged_edges, true);
3012
3013 MOFEM_LOG_C("WORLD", Sev::verbose,
3014 "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e", pp,
3015 nb_nodes_merged, ave, min);
3016
3017 if (debug)
3018 CHKERR save_merge_step(pp + 1, collapsed_edges);
3019
3020 if (nb_nodes_merged == nb_nodes_merged_p)
3021 break;
3022 if (min > 1e-2 && min == min_pp)
3023 break;
3024 if (min > ave0)
3025 break;
3026
3027 Range adj_edges;
3028 CHKERR moab.get_adjacencies(proc_tets, 1, false, adj_edges,
3029 moab::Interface::UNION);
3030 edges_to_merge = intersect(edges_to_merge, adj_edges);
3031 CHKERR Toplogy(m_field, th)
3032 .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
3033 constrainSurface, edges_to_merge, not_merged_edges);
3034 }
3035
3036 auto reconstruct_refined_ents = [&]() {
3038 CHKERR reconstructMultiIndex(*m_field.get_ref_ents());
3040 };
3041
3042 // Add function which reconstruct core multi-index. Node merging is messy
3043 // process and entity parent could be changed without notification to
3044 // multi-index. TODO Issue has to be tracked down better, however in principle
3045 // is better not to modify multi-index each time parent is changed, that makes
3046 // code slow. Is better to do it in the bulk as below.
3047 CHKERR reconstruct_refined_ents();
3048
3049 if (bit_ptr)
3050 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(proc_tets,
3051 *bit_ptr);
3052
3053 out_tets.merge(proc_tets);
3054 Range adj_faces;
3055 CHKERR moab.get_adjacencies(out_tets, 2, false, adj_faces,
3056 moab::Interface::UNION);
3057 new_surf = intersect(new_surf, adj_faces);
3058
3060}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:1716
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

◆ mergeSurface()

MoFEMErrorCode MoFEM::CutMeshInterface::mergeSurface ( const Range  surface)

merge surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 120 of file CutMeshInterface.cpp.

120 {
122 sUrface.merge(surface);
124}

◆ mergeVolumes()

MoFEMErrorCode MoFEM::CutMeshInterface::mergeVolumes ( const Range  volume)

merge volume entities

Parameters
volumeentities which going to be added
Returns
error code

Definition at line 126 of file CutMeshInterface.cpp.

126 {
128 vOlume.merge(volume);
130}

◆ moveMidNodesOnCutEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::moveMidNodesOnCutEdges ( Tag  th = NULL)

projecting of mid edge nodes on new mesh on surface

Returns
error code

Definition at line 1447 of file CutMeshInterface.cpp.

1447 {
1448 CoreInterface &m_field = cOre;
1449 moab::Interface &moab = m_field.get_moab();
1451
1452 // Range out_side_vertices;
1453 for (auto m : verticesOnCutEdges) {
1454 double dist = m.second.dIst;
1455 VectorDouble3 new_coors = m.second.rayPoint + dist * m.second.unitRayDir;
1456 if (th)
1457 CHKERR moab.tag_set_data(th, &m.first, 1, &new_coors[0]);
1458 else
1459 CHKERR moab.set_coords(&m.first, 1, &new_coors[0]);
1460 }
1461
1463}

◆ moveMidNodesOnTrimmedEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::moveMidNodesOnTrimmedEdges ( Tag  th = NULL)

move trimmed edges mid nodes

Returns
error code

Definition at line 1465 of file CutMeshInterface.cpp.

1465 {
1466 CoreInterface &m_field = cOre;
1467 moab::Interface &moab = m_field.get_moab();
1469 for (auto &v : verticesOnTrimEdges) {
1470 double dist = v.second.dIst;
1471 VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1472 if (th)
1473 CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1474 else
1475 CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1476 }
1478}

◆ projectZeroDistanceEnts() [1/2]

const Range & MoFEM::CutMeshInterface::projectZeroDistanceEnts ( ) const
inline

Definition at line 473 of file CutMeshInterface.hpp.

473 {
474 return zeroDistanceEnts;
475 }

◆ projectZeroDistanceEnts() [2/2]

MoFEMErrorCode MoFEM::CutMeshInterface::projectZeroDistanceEnts ( Range fixed_edges,
Range corner_nodes,
const double  geometry_tol = 0,
const double  close_tol = 0,
const int  verb = QUIET,
const bool  debug = false 
)

Find entities on cut surface which can be projected.

Parameters
fixed_edgespointer to fix edges
corner_nodespointer to corner nodes
closetolerance is tolerance how close entities has to be
verbverbosity level
debugtrue for debuging purposes

Definition at line 991 of file CutMeshInterface.cpp.

993 {
994 CoreInterface &m_field = cOre;
995 moab::Interface &moab = m_field.get_moab();
996 Skinner skin(&moab);
998
999 auto get_ent_adj = [&moab](const EntityHandle v, const int dim) {
1000 Range a;
1001 if (dim)
1002 CHKERR moab.get_adjacencies(&v, 1, dim, true, a);
1003 return a;
1004 };
1005
1006 auto get_adj = [&moab](const Range r, const int dim) {
1007 Range a;
1008 if (dim)
1009 CHKERR moab.get_adjacencies(r, dim, false, a, moab::Interface::UNION);
1010 else
1011 CHKERR moab.get_connectivity(r, a, true);
1012 return a;
1013 };
1014
1015 auto get_skin = [&skin](const auto r) {
1016 Range s;
1017 CHKERR skin.find_skin(0, r, false, s);
1018 return s;
1019 };
1020
1021
1022 auto get_range = [](std::vector<EntityHandle> v) {
1023 Range r;
1024 r.insert_list(v.begin(), v.end());
1025 return r;
1026 };
1027
1028 auto get_coords = [&](const auto v) {
1029 VectorDouble3 coords(3);
1030 CHKERR moab.get_coords(&v, 1, &coords[0]);
1031 return coords;
1032 };
1033
1034 Tag th_normal_dist;
1035 CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_VECTOR", th_normal_dist);
1036
1037 auto get_normal_dits = [&](const EntityHandle v) {
1038 VectorDouble3 dist_vec(3);
1039 CHKERR moab.tag_get_data(th_normal_dist, &v, 1, &*dist_vec.begin());
1040 return dist_vec;
1041 };
1042
1043 auto get_prj_point = [&](const EntityHandle v, const Range edges,
1044 const double geometry_tol) {
1045 auto get_tuple = [](const EntityHandle e, const double dist,
1046 const double l) { return std::make_tuple(e, dist, l); };
1047
1048 std::tuple<EntityHandle, double, double> min_tuple{0, 0, 0};
1049
1050 for (auto e : edges) {
1051 auto eit = edgesToCut.find(e);
1052 if (eit != edgesToCut.end()) {
1053
1054 auto get_dist = [&](auto eit) {
1055 int num_nodes;
1056 const EntityHandle *conn;
1057 CHKERR moab.get_connectivity(eit->first, conn, num_nodes, true);
1058 if (conn[0] == v)
1059 return eit->second.dIst;
1060 else if (conn[1] == v)
1061 return (eit->second.lEngth - eit->second.dIst);
1062 else
1063 THROW_MESSAGE("Data inconsistency");
1064 };
1065
1066 const auto d = get_dist(eit);
1067 if (std::get<0>(min_tuple) == 0) {
1068 min_tuple = get_tuple(e, d, eit->second.lEngth);
1069 } else if (std::get<1>(min_tuple) > d) {
1070 min_tuple = get_tuple(e, d, eit->second.lEngth);
1071 }
1072 }
1073 }
1074
1075 const auto geom_dist_vec = get_normal_dits(v);
1076 const auto geom_tol = norm_2(geom_dist_vec);
1077 const auto l = std::get<2>(min_tuple);
1078
1079 if (geom_tol < l * geometry_tol) {
1080 return std::make_pair(get_coords(v), l);
1081
1082 } else {
1083 const auto &d = edgesToCut.at(std::get<0>(min_tuple));
1084 return std::make_pair(VectorDouble3(d.rayPoint + d.dIst * d.unitRayDir),
1085 l);
1086 }
1087 };
1088
1089 auto get_in_range = [](auto v, auto &r) { return (r.find(v) != r.end()); };
1090
1091 auto project_nodes = [&](auto nodes_to_check) {
1092 if (debug)
1093 CHKERR SaveData(moab)("nodes_to_check.vtk", nodes_to_check);
1094
1095 auto get_fix_e = [](auto fixed_edges) {
1096 if (fixed_edges)
1097 return *fixed_edges;
1098 else
1099 return Range();
1100 };
1101
1102 const auto fix_e = get_fix_e(fixed_edges);
1103 const auto skin_e = get_adj(unite(get_skin(vOlume), constrainSurface), 1);
1104 const auto cut_fix = intersect(cutEdges, fix_e);
1105 const auto cut_skin = intersect(cutEdges, skin_e);
1106
1107 Range corner_n;
1108 if (corner_nodes)
1109 corner_n = intersect(*corner_nodes, nodes_to_check);
1110 auto fixe_n = intersect(get_adj(fix_e, 0), nodes_to_check);
1111 auto skin_n = intersect(get_adj(skin_e, 0), nodes_to_check);
1112
1113 std::vector<std::pair<EntityHandle, TreeData>> vertices_on_cut_edges;
1114 vertices_on_cut_edges.reserve(nodes_to_check.size());
1115
1116 auto add_zero_vertex = [&](const EntityHandle v, auto &&new_pos) {
1117 auto coords = get_coords(v);
1118 auto ray = new_pos - coords;
1119 auto l = norm_2(ray);
1120 TreeData data;
1121 data.dIst = l;
1122 if (l > std::numeric_limits<double>::epsilon()) {
1123 data.unitRayDir = ray / l;
1124 } else {
1125 data.dIst = 0;
1126 data.unitRayDir = ray;
1127 }
1128 data.rayPoint = coords;
1129 return std::make_pair(v, data);
1130 };
1131
1132 auto intersect_v = [&](const auto v, const auto r) {
1133 return intersect(r, get_ent_adj(v, 1));
1134 };
1135
1136 for (auto v : nodes_to_check) {
1137
1138 const auto e = intersect_v(v, cutEdges);
1139 if (!e.empty()) {
1140
1141 if (get_in_range(v, corner_n)) {
1142 auto p = get_prj_point(v, e, 0);
1143 if (norm_2(get_coords(v) - p.first) < close_tol * p.second) {
1144 vertices_on_cut_edges.push_back(add_zero_vertex(v, get_coords(v)));
1145 continue;
1146 }
1147
1148 } else if (get_in_range(v, fixe_n)) {
1149 const auto b = intersect_v(v, cut_fix);
1150 if (!b.empty()) {
1151 auto p = get_prj_point(v, b, geometry_tol);
1152 if (norm_2(get_coords(v) - p.first) < close_tol * p.second) {
1153 vertices_on_cut_edges.push_back(add_zero_vertex(v, p.first));
1154 continue;
1155 }
1156
1157 } else if (get_in_range(v, skin_n)) {
1158 const auto b = intersect_v(v, cut_skin);
1159 if (!b.empty()) {
1160 auto p = get_prj_point(v, b, geometry_tol);
1161 if (norm_2(get_coords(v) - p.first) < close_tol * p.second) {
1162 vertices_on_cut_edges.push_back(add_zero_vertex(v, p.first));
1163 continue;
1164 }
1165 }
1166
1167 }
1168 }
1169
1170 auto p = get_prj_point(v, e, geometry_tol);
1171 if (norm_2(get_coords(v) - p.first) < close_tol * p.second) {
1172 if (get_in_range(v, fixe_n) || get_in_range(v, skin_n))
1173 vertices_on_cut_edges.push_back(add_zero_vertex(v, get_coords(v)));
1174 else
1175 vertices_on_cut_edges.push_back(add_zero_vertex(v, p.first));
1176 }
1177 }
1178
1179 }
1180
1181 auto get_distances = [&](auto &data) {
1182 std::map<EntityHandle, double> dist_map;
1183 if (!data.empty()) {
1184 Tag th;
1185 CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_SIGNED", th);
1186
1187 std::vector<EntityHandle> verts;
1188 verts.reserve(data.size());
1189 for (auto d : data)
1190 verts.emplace_back(d.first);
1191 std::vector<double> distances(verts.size());
1192 CHKERR moab.tag_get_data(th, &*verts.begin(), verts.size(),
1193 &*distances.begin());
1194 for (size_t i = 0; i != verts.size(); ++i)
1195 dist_map[verts[i]] = distances[i];
1196 }
1197 return dist_map;
1198 };
1199
1200 auto dist_map = get_distances(vertices_on_cut_edges);
1201 if (!dist_map.empty()) {
1202 auto cmp = [&dist_map](const auto &a, const auto &b) {
1203 return dist_map[a.first] < dist_map[b.first];
1204 };
1205 std::sort(vertices_on_cut_edges.begin(), vertices_on_cut_edges.end(),
1206 cmp);
1207 }
1208
1209 return vertices_on_cut_edges;
1210 };
1211
1212 auto get_min_quality =
1213 [&](const Range &adj_tets,
1214 const map<EntityHandle, TreeData> &vertices_on_cut_edges,
1215 const std::pair<EntityHandle, TreeData> *test_data_ptr) {
1216 double min_q = 1;
1217 for (auto t : adj_tets) {
1218 int num_nodes;
1219 const EntityHandle *conn;
1220 CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1221 VectorDouble12 coords(12);
1222 CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1223
1224 auto set_new_coord = [](auto d) {
1225 return d.rayPoint + d.dIst * d.unitRayDir;
1226 };
1227
1228 for (auto n : {0, 1, 2, 3}) {
1229 auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1230 if (test_data_ptr) {
1231 if (test_data_ptr->first == conn[n]) {
1232 noalias(n_coords) = set_new_coord(test_data_ptr->second);
1233 }
1234 }
1235 auto mit = vertices_on_cut_edges.find(conn[n]);
1236 if (mit != vertices_on_cut_edges.end()) {
1237 noalias(n_coords) = set_new_coord(mit->second);
1238 }
1239 }
1240 min_q = std::min(min_q, Tools::volumeLengthQuality(&coords[0]));
1241 }
1242 return min_q;
1243 };
1244
1245 auto get_zero_distance_verts = [&](const auto &&vertices_on_cut_edges) {
1246 std::vector<EntityHandle> zero_dist_vec;
1247 zero_dist_vec.reserve(vertices_on_cut_edges.size());
1248 for (auto t : vertices_on_cut_edges) {
1249
1250 auto adj_tet = intersect(vOlume, get_ent_adj(t.first, 3));
1251 const double q0 = get_min_quality(adj_tet, verticesOnCutEdges, nullptr);
1252 const double q = get_min_quality(adj_tet, verticesOnCutEdges, &t);
1253 if (std::isnormal(q) && (q / q0) > projectEntitiesQualityTrashold) {
1254 verticesOnCutEdges[t.first] = t.second;
1255 zero_dist_vec.push_back(t.first);
1256 }
1257 }
1258 return zero_dist_vec;
1259 };
1260
1261 auto vol_cut_ents = intersect(vOlume, get_adj(cutEdges, 3));
1262
1263 auto get_zero_distant_ents = [&](auto bridge_ents, const int dim,
1264 const int bridge_dim) {
1265 auto ents =
1266 intersect(get_adj(bridge_ents, dim), get_adj(vol_cut_ents, dim));
1267 auto ents_to_remove =
1268 get_adj(subtract(get_adj(ents, bridge_dim), bridge_ents), dim);
1269 return subtract(ents, ents_to_remove);
1270 };
1271
1272 verticesOnCutEdges.clear();
1273 zeroDistanceVerts.clear();
1274 zeroDistanceVerts.merge(
1275 get_range(get_zero_distance_verts(project_nodes(get_adj(cutEdges, 0)))));
1276 zeroDistanceEnts = subtract(get_zero_distant_ents(zeroDistanceVerts, 2, 0),
1277 get_skin(vOlume));
1278
1279 if (debug)
1280 CHKERR SaveData(moab)("zero_distance_verts.vtk", zeroDistanceVerts);
1281 if (debug)
1282 CHKERR SaveData(moab)("zero_distance_ents.vtk", zeroDistanceEnts);
1283
1284 for (auto f : zeroDistanceVerts) {
1285 auto adj_edges = get_ent_adj(f, 1);
1286 for (auto e : adj_edges) {
1287 cutEdges.erase(e);
1288 edgesToCut.erase(e);
1289 }
1290 }
1291
1292 if (debug)
1293 CHKERR SaveData(moab)("cut_edges_to_cut.vtk", cutEdges);
1294
1296}
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
FTensor::Index< 'l', 3 > l

◆ query_interface()

MoFEMErrorCode MoFEM::CutMeshInterface::query_interface ( boost::typeindex::type_index  type_index,
UnknownInterface **  iface 
) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 16 of file CutMeshInterface.cpp.

17 {
19 *iface = const_cast<CutMeshInterface *>(this);
21}
CutMeshInterface(const MoFEM::Core &core)

◆ refineMesh()

MoFEMErrorCode MoFEM::CutMeshInterface::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.

Note
Should be run befor cutting
Parameters
refine_frontrefine nodes at front
update_frontupdate level set at front
init_bit_levelinital bit ref level to store refined meshes
surf_levelsnumber of mesh surface refinement
front_levels
fixed_edges
verb
debug
Returns
MoFEMErrorCode
Examples
mesh_cut.cpp.

Definition at line 809 of file CutMeshInterface.cpp.

813 {
814 CoreInterface &m_field = cOre;
815 moab::Interface &moab = m_field.get_moab();
816 MeshRefinement *refiner;
817 BitRefManager *bit_ref_manager;
819 CHKERR m_field.getInterface(refiner);
820 CHKERR m_field.getInterface(bit_ref_manager);
821
822 auto add_bit = [&](const int bit) {
824 CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
825 verb);
826 Range adj_ents;
827 for (auto d : {2, 1, 0}) {
828 CHKERR moab.get_adjacencies(vOlume, d, true, adj_ents,
829 moab::Interface::UNION);
830 CHKERR bit_ref_manager->addBitRefLevel(adj_ents, BitRefLevel().set(bit),
831 verb);
832 }
834 };
835 CHKERR add_bit(init_bit_level);
836
837 auto update_range = [&](Range *r_ptr) {
839 if (r_ptr) {
840 Range childs;
841 CHKERR bit_ref_manager->updateRangeByChildren(*r_ptr, childs);
842 r_ptr->merge(childs);
843 }
845 };
846
847 auto refine = [&](const BitRefLevel bit, const Range tets) {
849 Range verts;
850 CHKERR moab.get_connectivity(tets, verts, true);
851 Range ref_edges;
852 CHKERR moab.get_adjacencies(verts, 1, true, ref_edges,
853 moab::Interface::UNION);
854
855 CHKERR refiner->addVerticesInTheMiddleOfEdges(ref_edges, bit);
856 CHKERR refiner->refineTets(vOlume, bit, verb);
857
858 CHKERR update_range(fixed_edges);
859 CHKERR update_range(&vOlume);
860 CHKERR m_field.getInterface<MeshsetsManager>()
861 ->updateAllMeshsetsByEntitiesChildren(bit);
862
863 Range bit_tets;
864 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
865 bit, BitRefLevel().set(), MBTET, bit_tets);
866 vOlume = intersect(vOlume, bit_tets);
867
868 Range bit_edges;
869 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
870 bit, BitRefLevel().set(), MBEDGE, bit_edges);
871 if (fixed_edges)
872 *fixed_edges = intersect(*fixed_edges, bit_edges);
873
875 };
876
877 for (int ll = init_bit_level; ll != init_bit_level + surf_levels; ++ll) {
879 CHKERR refine(BitRefLevel().set(ll + 1),
881 }
882
883 for (int ll = init_bit_level + surf_levels;
884 ll != init_bit_level + surf_levels + front_levels; ++ll) {
886 CHKERR refine(BitRefLevel().set(ll + 1), cutFrontVolumes);
887 }
888
890
891 if (debug)
892 CHKERR SaveData(m_field.get_moab())("refinedVolume.vtk", vOlume);
893
895}

◆ removePathologicalFrontTris()

MoFEMErrorCode MoFEM::CutMeshInterface::removePathologicalFrontTris ( const BitRefLevel  split_bit,
Range ents 
)

Remove pathological elements on surface internal front.

Internal surface skin is a set of edges in iterior of the body on boundary of surface. This set of edges is called surface front. If surface face has three nodes on surface front, non of the face nodes is split and should be removed from surface if it is going to be split.

Parameters
split_bitsplit bit level
bitbit level of split mesh
entsents on the surface which is going to be split
Returns
error code

Definition at line 2235 of file CutMeshInterface.cpp.

2236 {
2237 CoreInterface &m_field = cOre;
2238 moab::Interface &moab = m_field.get_moab();
2239 PrismInterface *interface;
2241 CHKERR m_field.getInterface(interface);
2242 // Remove tris on surface front
2243 {
2244 Range front_tris;
2245 EntityHandle meshset;
2246 CHKERR moab.create_meshset(MESHSET_SET, meshset);
2247 CHKERR moab.add_entities(meshset, ents);
2248 CHKERR interface->findFacesWithThreeNodesOnInternalSurfaceSkin(
2249 meshset, split_bit, true, front_tris);
2250 CHKERR moab.delete_entities(&meshset, 1);
2251 ents = subtract(ents, front_tris);
2252 }
2253 // Remove entities on skin
2254 Range tets;
2255 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2256 split_bit, BitRefLevel().set(), MBTET, tets);
2257 // Remove entities on skin
2258 Skinner skin(&moab);
2259 Range tets_skin;
2260 CHKERR skin.find_skin(0, tets, false, tets_skin);
2261 ents = subtract(ents, tets_skin);
2262
2264}

◆ saveCutEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveCutEdges ( const std::string  prefix = "")

Definition at line 3147 of file CutMeshInterface.cpp.

3147 {
3148 CoreInterface &m_field = cOre;
3149 moab::Interface &moab = m_field.get_moab();
3151 CHKERR SaveData(moab)(prefix + "out_vol.vtk", vOlume);
3152 CHKERR SaveData(moab)(prefix + "out_surface.vtk", sUrface);
3153 CHKERR SaveData(moab)(prefix + "out_cut_edges.vtk", cutEdges);
3154 CHKERR SaveData(moab)(prefix + "out_cut_new_volumes.vtk", cutNewVolumes);
3155 CHKERR SaveData(moab)(prefix + "out_cut_new_surfaces.vtk", cutNewSurfaces);
3156 CHKERR SaveData(moab)(prefix + "out_cut_zero_distance_ents.vtk",
3159}

◆ saveTrimEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveTrimEdges ( )

Definition at line 3161 of file CutMeshInterface.cpp.

3161 {
3162 moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
3164 CHKERR SaveData(moab)("out_trim_surface.vtk", sUrface);
3165 CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
3166 CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
3167 CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
3169}

◆ setConstrainSurface()

MoFEMErrorCode MoFEM::CutMeshInterface::setConstrainSurface ( const Range  surf)

Set the constrain surface object.

Add surfaces which are restricted by mesh cut. Example of surface which needs to be respected is an interface, the boundary between two materials, or crack surface.

Parameters
surf
Returns
MoFEMErrorCode

Definition at line 114 of file CutMeshInterface.cpp.

114 {
116 constrainSurface = surf;
118}

◆ setCoords()

MoFEMErrorCode MoFEM::CutMeshInterface::setCoords ( Tag  th,
const BitRefLevel  bit = BitRefLevel(),
const BitRefLevel  mask = BitRefLevel().set() 
)

set coords from tag

Parameters
thtag handle
Returns
error code
Examples
mesh_cut.cpp.

Definition at line 3130 of file CutMeshInterface.cpp.

3131 {
3132 CoreInterface &m_field = cOre;
3133 moab::Interface &moab = m_field.get_moab();
3135 Range nodes;
3136 if (bit.none())
3137 CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3138 else
3139 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3140 bit, mask, MBVERTEX, nodes);
3141 std::vector<double> coords(3 * nodes.size());
3142 CHKERR moab.tag_get_data(th, nodes, &coords[0]);
3143 CHKERR moab.set_coords(nodes, &coords[0]);
3145}

◆ setFront()

MoFEMErrorCode MoFEM::CutMeshInterface::setFront ( const Range  surface)

Definition at line 36 of file CutMeshInterface.cpp.

36 {
38 fRont = front;
40}

◆ setSurface()

MoFEMErrorCode MoFEM::CutMeshInterface::setSurface ( const Range  surface)

set surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 42 of file CutMeshInterface.cpp.

◆ setTagData()

MoFEMErrorCode MoFEM::CutMeshInterface::setTagData ( Tag  th,
const BitRefLevel  bit = BitRefLevel() 
)

set coords to tag

Parameters
thtag handle
Returns
error code
Examples
mesh_cut.cpp.

Definition at line 3114 of file CutMeshInterface.cpp.

3114 {
3115 CoreInterface &m_field = cOre;
3116 moab::Interface &moab = m_field.get_moab();
3118 Range nodes;
3119 if (bit.none())
3120 CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3121 else
3122 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3123 bit, BitRefLevel().set(), MBVERTEX, nodes);
3124 std::vector<double> coords(3 * nodes.size());
3125 CHKERR moab.get_coords(nodes, &coords[0]);
3126 CHKERR moab.tag_set_data(th, nodes, &coords[0]);
3128}

◆ setVolume()

MoFEMErrorCode MoFEM::CutMeshInterface::setVolume ( const Range  volume)

set volume entities

Parameters
volumeentities which going to be added
Returns
error code
Examples
mesh_cut.cpp.

Definition at line 108 of file CutMeshInterface.cpp.

108 {
110 vOlume = volume;
112}

◆ snapSurfaceSkinToEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::snapSurfaceSkinToEdges ( const Range  fixed_edges,
const double  rel_tol,
const double  abs_tol,
Tag  th = nullptr,
const bool  debug = false 
)

Definition at line 132 of file CutMeshInterface.cpp.

136 {
137 CoreInterface &m_field = cOre;
138 moab::Interface &moab = m_field.get_moab();
140
141 // Get cutting surface skin
142 Skinner skin(&moab);
143 Range surface_skin;
144 CHKERR skin.find_skin(0, sUrface, false, surface_skin);
145
146 CHKERR snapSurfaceToEdges(surface_skin, fixed_edges, rel_tol, abs_tol, th,
147 debug);
148
150}
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)

◆ snapSurfaceToEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::snapSurfaceToEdges ( const Range  surface_edges,
const Range  fixed_edges,
const double  rel_tol,
const double  abs_tol,
Tag  th = nullptr,
const bool  debug = false 
)

Definition at line 152 of file CutMeshInterface.cpp.

156 {
157 CoreInterface &m_field = cOre;
158 moab::Interface &moab = m_field.get_moab();
161
162 map<EntityHandle, double> map_verts_length;
163
164 for (auto f : surface_edges) {
165 int num_nodes;
166 const EntityHandle *conn_skin;
167 CHKERR moab.get_connectivity(f, conn_skin, num_nodes, true);
168 VectorDouble6 coords_skin(6);
169 if (th)
170 CHKERR moab.tag_get_data(th, conn_skin, num_nodes, &coords_skin[0]);
171 else
172 CHKERR moab.get_coords(conn_skin, num_nodes, &coords_skin[0]);
174 &coords_skin[0], &coords_skin[1], &coords_skin[2]);
176 &coords_skin[3], &coords_skin[4], &coords_skin[5]);
177 FTensor::Tensor1<double, 3> t_skin_delta;
178 t_skin_delta(i) = t_n1(i) - t_n0(i);
179 const double skin_edge_length = sqrt(t_skin_delta(i) * t_skin_delta(i));
180 for (int nn = 0; nn != 2; ++nn) {
181 auto it = map_verts_length.find(conn_skin[nn]);
182 if (it != map_verts_length.end())
183 it->second = std::min(it->second, skin_edge_length);
184 else
185 map_verts_length[conn_skin[nn]] = skin_edge_length;
186 }
187 }
188
189 for (auto m : map_verts_length) {
190
192 if (th)
193 CHKERR moab.tag_get_data(th, &m.first, 1, &t_n(0));
194 else
195 CHKERR moab.get_coords(&m.first, 1, &t_n(0));
196
197 double min_dist = rel_tol * m.second;
198 FTensor::Tensor1<double, 3> t_min_coords;
199 CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
200 &t_n(0), 1, fixed_edges, &min_dist, &t_min_coords(0));
201
202 if (min_dist < rel_tol * m.second || min_dist < abs_tol) {
203 if (debug)
204 MOFEM_LOG("WORLD", Sev::noisy) << "Snap " << min_dist;
205 if (th)
206 CHKERR moab.tag_set_data(th, &m.first, 1, &t_min_coords(0));
207 else
208 CHKERR moab.set_coords(&m.first, 1, &t_min_coords(0));
209 }
210 }
211
213}
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308

◆ splitSides()

MoFEMErrorCode MoFEM::CutMeshInterface::splitSides ( const BitRefLevel  split_bit,
const BitRefLevel  bit,
const Range ents,
Tag  th = NULL 
)

split sides

Parameters
split_bitsplit bit level
bitbit level of split mesh
entsents on the surface which is going to be split
Returns
error code

Definition at line 2266 of file CutMeshInterface.cpp.

2268 {
2269 CoreInterface &m_field = cOre;
2270 moab::Interface &moab = m_field.get_moab();
2271 PrismInterface *interface;
2273 CHKERR m_field.getInterface(interface);
2274 EntityHandle meshset_volume;
2275 CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
2276 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2277 split_bit, BitRefLevel().set(), MBTET, meshset_volume);
2278 EntityHandle meshset_surface;
2279 CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
2280 CHKERR moab.add_entities(meshset_surface, ents);
2281 CHKERR interface->getSides(meshset_surface, split_bit, true);
2282 CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
2283 true);
2284 CHKERR moab.delete_entities(&meshset_volume, 1);
2285 CHKERR moab.delete_entities(&meshset_surface, 1);
2286 if (th) {
2287 Range prisms;
2288 ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2289 bit, BitRefLevel().set(), MBPRISM, prisms);
2290 for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
2291 int num_nodes;
2292 const EntityHandle *conn;
2293 CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
2294 MatrixDouble data(3, 3);
2295 CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
2296 CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
2297 }
2298 }
2300}

◆ trimEdgesInTheMiddle()

MoFEMErrorCode MoFEM::CutMeshInterface::trimEdgesInTheMiddle ( const BitRefLevel  bit,
const bool  debug = false 
)

trim edges

Parameters
bitbit level of the trimmed mesh
Returns
error code

Definition at line 1969 of file CutMeshInterface.cpp.

1970 {
1971 CoreInterface &m_field = cOre;
1972 moab::Interface &moab = m_field.get_moab();
1973 MeshRefinement *refiner;
1974 auto ref_ents_ptr = m_field.get_ref_ents();
1976
1977 CHKERR m_field.getInterface(refiner);
1978 CHKERR refiner->addVerticesInTheMiddleOfEdges(trimEdges, bit);
1979 CHKERR refiner->refineTets(cutNewVolumes, bit, QUIET, debug);
1980
1981 trimNewVolumes.clear();
1982 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1983 bit, BitRefLevel().set(), MBTET, trimNewVolumes);
1984
1985 for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
1986 mit != edgesToTrim.end(); mit++) {
1987 auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1988 boost::make_tuple(MBVERTEX, mit->first));
1989 if (vit ==
1990 ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end())
1991 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1992 "No vertex on trim edges, that make no sense");
1993
1994 const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1995 if ((ref_ent->getBitRefLevel() & bit).any()) {
1996 EntityHandle vert = ref_ent->getEnt();
1997 trimNewVertices.insert(vert);
1998 verticesOnTrimEdges[vert] = mit->second;
1999 }
2000 }
2001
2002 // Get faces which are trimmed
2003 trimNewSurfaces.clear();
2004 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2005 bit, bit, MBTRI, trimNewSurfaces);
2006
2007 Range trim_new_surfaces_nodes;
2008 CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
2009 trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
2010 trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
2011 Range faces_not_on_surface;
2012 CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
2013 faces_not_on_surface, moab::Interface::UNION);
2014 trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
2015
2016 // Get surfaces which are not trimmed and add them to surface
2017 Range all_surfaces_on_bit_level;
2018 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2019 bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
2020 all_surfaces_on_bit_level =
2021 intersect(all_surfaces_on_bit_level, cutNewSurfaces);
2022
2023 trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
2024
2026}

◆ trimOnly()

MoFEMErrorCode MoFEM::CutMeshInterface::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.

Parameters
trim_bit
th
tol_cut_close
fixed_edges
corner_nodes
update_meshsets
debug
Returns
MoFEMErrorCode

Definition at line 263 of file CutMeshInterface.cpp.

268 {
269 CoreInterface &m_field = cOre;
270 moab::Interface &moab = m_field.get_moab();
272
273 // trim mesh
274 CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim_close, debug);
276
277 CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(constrainSurface,
279 if (fixed_edges)
280 CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(*fixed_edges,
281 *fixed_edges);
282
283 if (corner_nodes)
284 CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(*corner_nodes,
285 *corner_nodes);
286
287 if (update_meshsets)
288 CHKERR m_field.getInterface<MeshsetsManager>()
289 ->updateAllMeshsetsByEntitiesChildren(trim_bit);
290
291 // move nodes
293
294 // remove faces
295 CHKERR trimSurface(fixed_edges, corner_nodes, tol_trim_close, th, debug);
296
297 if (debug) {
299 Range bit2_edges;
300 CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
301 trim_bit, BitRefLevel().set(), MBEDGE, bit2_edges);
302 CHKERR SaveData(moab)("trim_fixed_edges.vtk",
303 intersect(*fixed_edges, bit2_edges));
304 }
305
307}
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 moveMidNodesOnTrimmedEdges(Tag th=NULL)
move trimmed edges mid nodes
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.
MoFEMErrorCode trimEdgesInTheMiddle(const BitRefLevel bit, const bool debug=false)
trim edges
MoFEMErrorCode saveTrimEdges()

◆ trimSurface()

MoFEMErrorCode MoFEM::CutMeshInterface::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.

Parameters
fixed_edge
corner_nodes
debug
Returns
MoFEMErrorCode

Definition at line 2028 of file CutMeshInterface.cpp.

2031 {
2032 CoreInterface &m_field = cOre;
2033 moab::Interface &moab = m_field.get_moab();
2034 Skinner skin(&moab);
2035
2036 auto get_adj = [&moab](const Range r, const int dim) {
2037 Range a;
2038 if (dim)
2039 CHKERR moab.get_adjacencies(r, dim, false, a, moab::Interface::UNION);
2040 else
2041 CHKERR moab.get_connectivity(r, a, true);
2042 return a;
2043 };
2044
2045 auto get_skin = [&skin](const auto r) {
2046 Range s;
2047 CHKERR skin.find_skin(0, r, false, s);
2048 return s;
2049 };
2050
2052
2053 auto trim_tets_skin = get_skin(trimNewVolumes);
2054 auto trim_tets_skin_edges = get_adj(trim_tets_skin, 1);
2055 auto trim_surface_edges = get_adj(trimNewSurfaces, 1);
2056
2057 auto contarain_edges =
2058 intersect(trim_surface_edges, get_adj(constrainSurface, 1));
2059 if (fixed_edges)
2060 contarain_edges.merge(
2061 intersect(fixed_edges->subset_by_type(MBEDGE), trim_surface_edges));
2062
2063 Range barrier_vertices(trimNewVertices);
2064 if (corner_nodes) {
2065 // Add nodes which re adjacent to corner nodes
2066 barrier_vertices.merge(get_adj(
2067 intersect(get_adj(intersect(*corner_nodes, barrier_vertices), 2),
2069 0));
2070 }
2071
2072 auto get_nodes_with_one_node_on_fixed_edge_other_not = [&]() {
2073 const auto fixed_edges_on_surface =
2074 intersect(*fixed_edges, trim_surface_edges);
2075 const auto skin_fixed_edges_on_surface = get_skin(fixed_edges_on_surface);
2076 const auto barrier_nodes = subtract(skin_fixed_edges_on_surface,
2077 get_adj(get_skin(trimNewSurfaces), 0));
2078 return barrier_nodes;
2079 };
2080
2081 if (fixed_edges)
2082 barrier_vertices.merge(get_nodes_with_one_node_on_fixed_edge_other_not());
2083
2084 auto add_close_surface_barrier = [&]() {
2086
2087 // Ten algorytm dedykuje destylarni Bowmore. Jeżeli coś pójdzie nie tak, coś
2088 // pęknie inaczej niż trzeba, to pełną odpowiedzialność ponosi Beam Suntory
2089 // UK Ltd.
2090
2091 CHKERR buildTree();
2092
2093 auto test_edges =
2094 subtract(trim_surface_edges, get_adj(constrainSurface, 1));
2095 if (fixed_edges)
2096 test_edges = subtract(test_edges, *fixed_edges);
2097 auto trim_surface_nodes = get_adj(test_edges, 0);
2098 trim_surface_nodes = subtract(trim_surface_nodes, barrier_vertices);
2099 if (fixed_edges)
2100 trim_surface_nodes =
2101 subtract(trim_surface_nodes, get_adj(*fixed_edges, 0));
2102
2103 Range trim_skin;
2104 trim_skin.merge(contarain_edges);
2105 if (fRont.empty())
2106 trim_skin.merge(get_skin(sUrface));
2107 else
2108 trim_skin.merge(fRont);
2109 if (fixed_edges)
2110 trim_skin.merge(intersect(*fixed_edges, trim_surface_edges));
2111
2112 if (debug && !trim_skin.empty())
2113 CHKERR SaveData(m_field.get_moab())("trim_skin.vtk", trim_skin);
2114
2115 VectorDouble coords(3 * trim_surface_nodes.size());
2116 if (th)
2117 CHKERR moab.tag_get_data(th, trim_surface_nodes, &*coords.begin());
2118 else
2119 CHKERR moab.get_coords(trim_surface_nodes, &*coords.begin());
2120
2121 if (!trim_skin.empty()) {
2122
2123 std::vector<double> min_distances(trim_surface_nodes.size(), -1);
2124 CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
2125 &*coords.begin(), trim_surface_nodes.size(), trim_skin,
2126 &*min_distances.begin());
2127
2128 auto coords_ptr = &*coords.begin();
2129 auto min_dist = &*min_distances.begin();
2130
2131 std::vector<EntityHandle> add_nodes;
2132 add_nodes.reserve(trim_surface_nodes.size());
2133
2134 for (auto n : trim_surface_nodes) {
2135
2136 if ((*min_dist) > std::numeric_limits<double>::epsilon()) {
2137 VectorDouble3 point_out(3);
2138
2139 EntityHandle facets_out;
2140 CHKERR treeSurfPtr->closest_to_location(coords_ptr, rootSetSurf,
2141 &point_out[0], facets_out);
2142
2143 VectorDouble3 delta = point_out - getVectorAdaptor(coords_ptr, 3);
2144 VectorDouble3 normal(3);
2145 CHKERR m_field.getInterface<Tools>()->getTriNormal(facets_out,
2146 &normal[0]);
2147
2148 normal /= norm_2(normal);
2149 delta -= inner_prod(normal, delta) * normal;
2150
2151 double dist = norm_2(delta);
2152 if (dist < tol * (*min_dist))
2153 add_nodes.emplace_back(n);
2154 }
2155
2156 coords_ptr += 3;
2157 min_dist += 1;
2158 }
2159
2160 barrier_vertices.insert_list(add_nodes.begin(), add_nodes.end());
2161 }
2162
2164 };
2165
2166 auto remove_faces_on_skin = [&]() {
2168 Range skin_faces = intersect(trimNewSurfaces, trim_tets_skin);
2169 if (!skin_faces.empty()) {
2170 barrier_vertices.merge(get_adj(skin_faces, 0));
2171 for (auto f : skin_faces)
2172 trimNewSurfaces.erase(f);
2173 }
2175 };
2176
2177 auto get_trim_free_edges = [&]() {
2178 // get current surface skin
2179 Range trim_surf_skin;
2180 CHKERR skin.find_skin(0, trimNewSurfaces, false, trim_surf_skin);
2181 trim_surf_skin = subtract(trim_surf_skin, trim_tets_skin_edges);
2182
2183 Range trim_surf_skin_verts;
2184 CHKERR moab.get_connectivity(trim_surf_skin, trim_surf_skin_verts, true);
2185 for (auto e : barrier_vertices)
2186 trim_surf_skin_verts.erase(e);
2187
2188 Range free_edges;
2189 CHKERR moab.get_adjacencies(trim_surf_skin_verts, 1, false, free_edges,
2190 moab::Interface::UNION);
2191 free_edges =
2192 subtract(intersect(free_edges, trim_surf_skin), contarain_edges);
2193
2194 return free_edges;
2195 };
2196
2197 CHKERR remove_faces_on_skin();
2198 CHKERR add_close_surface_barrier();
2199
2200 if (debug && !barrier_vertices.empty())
2201 CHKERR SaveData(m_field.get_moab())("barrier_vertices.vtk",
2202 barrier_vertices);
2203
2204 int nn = 0;
2205 Range out_edges;
2206 while (!(out_edges = get_trim_free_edges()).empty()) {
2207
2208 if (debug && !trimNewSurfaces.empty())
2209 CHKERR SaveData(m_field.get_moab())(
2210 "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
2212
2213 if (debug && !out_edges.empty())
2214 CHKERR SaveData(m_field.get_moab())(
2215 "trimNewSurfacesOutsideVerts_" +
2216 boost::lexical_cast<std::string>(nn) + ".vtk",
2217 out_edges);
2218
2219 Range outside_faces;
2220 CHKERR moab.get_adjacencies(out_edges, 2, false, outside_faces,
2221 moab::Interface::UNION);
2222 trimNewSurfaces = subtract(trimNewSurfaces, outside_faces);
2223 ++nn;
2224 }
2225
2226 if (debug && !trimNewSurfaces.empty())
2227 CHKERR SaveData(m_field.get_moab())(
2228 "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
2230
2232}
double tol
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEMErrorCode buildTree()
build tree

Member Data Documentation

◆ aveLength

double MoFEM::CutMeshInterface::aveLength
private

Average edge length.

Definition at line 572 of file CutMeshInterface.hpp.

◆ constrainSurface

Range MoFEM::CutMeshInterface::constrainSurface
private

Definition at line 577 of file CutMeshInterface.hpp.

◆ cOre

MoFEM::Core& MoFEM::CutMeshInterface::cOre

Definition at line 23 of file CutMeshInterface.hpp.

◆ cutEdges

Range MoFEM::CutMeshInterface::cutEdges
private

Definition at line 545 of file CutMeshInterface.hpp.

◆ cutFrontVolumes

Range MoFEM::CutMeshInterface::cutFrontVolumes
private

Definition at line 576 of file CutMeshInterface.hpp.

◆ cutNewSurfaces

Range MoFEM::CutMeshInterface::cutNewSurfaces
private

Definition at line 547 of file CutMeshInterface.hpp.

◆ cutNewVertices

Range MoFEM::CutMeshInterface::cutNewVertices
private

Definition at line 550 of file CutMeshInterface.hpp.

◆ cutNewVolumes

Range MoFEM::CutMeshInterface::cutNewVolumes
private

Definition at line 546 of file CutMeshInterface.hpp.

◆ cutSurfaceVolumes

Range MoFEM::CutMeshInterface::cutSurfaceVolumes
private

Definition at line 575 of file CutMeshInterface.hpp.

◆ edgesToCut

map<EntityHandle, TreeData> MoFEM::CutMeshInterface::edgesToCut
private

Definition at line 567 of file CutMeshInterface.hpp.

◆ edgesToTrim

map<EntityHandle, TreeData> MoFEM::CutMeshInterface::edgesToTrim
private

Definition at line 569 of file CutMeshInterface.hpp.

◆ fRont

Range MoFEM::CutMeshInterface::fRont
private

Definition at line 538 of file CutMeshInterface.hpp.

◆ lineSearchSteps

int MoFEM::CutMeshInterface::lineSearchSteps

Definition at line 27 of file CutMeshInterface.hpp.

◆ maxLength

double MoFEM::CutMeshInterface::maxLength
private

Maximal edge length.

Definition at line 573 of file CutMeshInterface.hpp.

◆ mergedSurfaces

Range MoFEM::CutMeshInterface::mergedSurfaces
private

Definition at line 558 of file CutMeshInterface.hpp.

◆ mergedVolumes

Range MoFEM::CutMeshInterface::mergedVolumes
private

Definition at line 557 of file CutMeshInterface.hpp.

◆ nbMaxMergingCycles

int MoFEM::CutMeshInterface::nbMaxMergingCycles

Definition at line 28 of file CutMeshInterface.hpp.

◆ projectEntitiesQualityTrashold

double MoFEM::CutMeshInterface::projectEntitiesQualityTrashold

Definition at line 29 of file CutMeshInterface.hpp.

◆ rootSetSurf

EntityHandle MoFEM::CutMeshInterface::rootSetSurf
private

Definition at line 543 of file CutMeshInterface.hpp.

◆ sUrface

Range MoFEM::CutMeshInterface::sUrface
private

Definition at line 539 of file CutMeshInterface.hpp.

◆ treeSurfPtr

boost::shared_ptr<OrientedBoxTreeTool> MoFEM::CutMeshInterface::treeSurfPtr
private

Definition at line 542 of file CutMeshInterface.hpp.

◆ trimEdges

Range MoFEM::CutMeshInterface::trimEdges
private

Definition at line 555 of file CutMeshInterface.hpp.

◆ trimNewSurfaces

Range MoFEM::CutMeshInterface::trimNewSurfaces
private

Definition at line 554 of file CutMeshInterface.hpp.

◆ trimNewVertices

Range MoFEM::CutMeshInterface::trimNewVertices
private

Definition at line 553 of file CutMeshInterface.hpp.

◆ trimNewVolumes

Range MoFEM::CutMeshInterface::trimNewVolumes
private

Definition at line 552 of file CutMeshInterface.hpp.

◆ verticesOnCutEdges

map<EntityHandle, TreeData> MoFEM::CutMeshInterface::verticesOnCutEdges
private

Definition at line 568 of file CutMeshInterface.hpp.

◆ verticesOnTrimEdges

map<EntityHandle, TreeData> MoFEM::CutMeshInterface::verticesOnTrimEdges
private

Definition at line 570 of file CutMeshInterface.hpp.

◆ vOlume

Range MoFEM::CutMeshInterface::vOlume
private

Definition at line 540 of file CutMeshInterface.hpp.

◆ zeroDistanceEnts

Range MoFEM::CutMeshInterface::zeroDistanceEnts
private

Definition at line 548 of file CutMeshInterface.hpp.

◆ zeroDistanceVerts

Range MoFEM::CutMeshInterface::zeroDistanceVerts
private

Definition at line 549 of file CutMeshInterface.hpp.


The documentation for this struct was generated from the following files: