v0.13.1
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:1086

◆ ~CutMeshInterface()

MoFEM::CutMeshInterface::~CutMeshInterface ( )

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
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
virtual moab::Interface & get_moab()=0
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}
Tensor2_Expr< transform_Tensor2< A, B, T, Dim0, Dim1, i, j >, T, Dim0, Dim1, i, j > transform(const Tensor2_Expr< A, T, Dim0, Dim1, i, j > &a, B function)
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:304
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 1310 of file CutMeshInterface.cpp.

1314 {
1315 CoreInterface &m_field = cOre;
1316 moab::Interface &moab = m_field.get_moab();
1317 MeshRefinement *refiner;
1318 auto ref_ents_ptr = m_field.get_ref_ents();
1320
1321 if (cutEdges.size() != edgesToCut.size())
1322 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1323
1324 auto refine_mesh = [&]() {
1326 CHKERR m_field.getInterface(refiner);
1327 CHKERR refiner->addVerticesInTheMiddleOfEdges(cutEdges, bit);
1328 CHKERR refiner->refineTets(vOlume, bit, QUIET, debug);
1330 };
1331
1332 CHKERR refine_mesh();
1333
1334 if (debug)
1335 if (cutEdges.size() != edgesToCut.size())
1336 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1337
1338 cut_vols.clear();
1339 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1340 bit, BitRefLevel().set(), MBTET, cut_vols);
1341 cut_surf.clear();
1342 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1343 bit, bit, MBTRI, cut_surf);
1344
1345 if (debug)
1346 CHKERR SaveData(moab)("cut_surf_from_bit.vtk", cut_surf);
1347
1348 // Find new vertices on cut edges
1349 cut_verts.clear();
1350 cut_verts.merge(zeroDistanceVerts);
1351
1352 for (auto &m : edgesToCut) {
1353 auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1354 boost::make_tuple(MBVERTEX, m.first));
1355 if (vit ==
1356 ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1357 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1358 "No vertex on cut edges, that make no sense");
1359 }
1360 const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1361 if ((ref_ent->getBitRefLevel() & bit).any()) {
1362 EntityHandle vert = ref_ent->getEnt();
1363 cut_verts.insert(vert);
1364 verticesOnCutEdges[vert] = m.second;
1365 } else {
1366 SETERRQ1(
1367 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1368 "Vertex has wrong bit ref level %s",
1369 boost::lexical_cast<std::string>(ref_ent->getBitRefLevel()).c_str());
1370 }
1371 }
1372
1373 // Add zero distance entities faces
1374 Range tets_skin;
1375 Skinner skin(&moab);
1376 CHKERR skin.find_skin(0, cut_vols, false, tets_skin);
1377 tets_skin.merge(constrainSurface);
1378
1379 // At that point cut_surf has all newly created faces, now take all
1380 // nodes on those faces and subtract nodes on cut edges. Faces adjacent to
1381 // nodes which left are not part of surface.
1382 cut_surf.merge(zeroDistanceEnts.subset_by_type(MBTRI));
1383 Range diff_verts;
1384 CHKERR moab.get_connectivity(cut_surf, diff_verts, true);
1385 diff_verts = subtract(diff_verts, cut_verts);
1386
1387 Range subtract_faces;
1388 CHKERR moab.get_adjacencies(diff_verts, 2, false, subtract_faces,
1389 moab::Interface::UNION);
1390 cut_surf = subtract(cut_surf, unite(subtract_faces, tets_skin));
1391
1392 cut_verts.clear();
1393 CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1394
1395 // Check non-mainfolds
1396 auto check_for_non_minfold = [&]() {
1398 Range surf_edges;
1399 CHKERR moab.get_adjacencies(cut_surf, 1, false, surf_edges,
1400 moab::Interface::UNION);
1401 for (auto e : surf_edges) {
1402
1403 Range faces;
1404 CHKERR moab.get_adjacencies(&e, 1, 2, false, faces);
1405 faces = intersect(faces, cut_surf);
1406 if (faces.size() > 2) {
1407
1408 bool resolved = false;
1409
1410 // Check for haning node
1411 Range nodes;
1412 CHKERR moab.get_connectivity(faces, nodes, true);
1413 for (auto n : nodes) {
1414 Range adj_faces;
1415 CHKERR moab.get_adjacencies(&n, 1, 2, false, adj_faces);
1416 adj_faces = intersect(adj_faces, cut_surf);
1417 if (adj_faces.size() == 1) {
1418 cut_surf.erase(adj_faces[0]);
1419 resolved = true;
1420 }
1421 }
1422
1423 // Check for two edges minfold
1424 Range adj_edges;
1425 CHKERR moab.get_adjacencies(faces, 1, false, adj_edges,
1426 moab::Interface::UNION);
1427 adj_edges = intersect(adj_edges, surf_edges);
1428 adj_edges.erase(e);
1429 for (auto other_e : adj_edges) {
1430 Range other_faces;
1431 CHKERR moab.get_adjacencies(&other_e, 1, 2, false, other_faces);
1432 other_faces = intersect(other_faces, cut_surf);
1433 if (other_faces.size() > 2) {
1434 other_faces = intersect(other_faces, faces);
1435 cut_surf = subtract(cut_surf, other_faces);
1436 resolved = true;
1437 }
1438 }
1439
1440 if (!resolved && !debug)
1441 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1442 "Non-mainfold surface");
1443
1444 cut_verts.clear();
1445 CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1446 }
1447 }
1449 };
1450
1451 CHKERR check_for_non_minfold();
1452
1453 if (debug)
1454 CHKERR SaveData(moab)("cut_surf.vtk", cut_surf);
1455
1457}
@ 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 get_range = [](std::vector<EntityHandle> v) {
932 Range r;
933 r.insert_list(v.begin(), v.end());
934 return r;
935 };
936
937 auto vol_edges = get_adj(vol, 1);
938
939 aveLength = 0;
940 maxLength = 0;
941 int nb_ave_length = 0;
942 for (auto e : vol_edges) {
943
944 auto conn = get_conn(e);
945 auto dist = get_tag_edge_dist(th_signed_dist, conn);
946 const auto dist_max = std::max(dist[0], dist[1]);
947 const auto dist_min = std::min(dist[0], dist[1]);
948 const auto dot = dist[0] * dist[1];
949
950 VectorDouble6 coords(6);
951 CHKERR moab.get_coords(conn, 2, &coords[0]);
952 VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
953 VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
954 VectorDouble3 ray = n1 - n0;
955 const double ray_length = norm_2(ray);
956 ray /= ray_length;
957
958 if (
959
960 (dot < 0 && dist_max > std::numeric_limits<double>::epsilon()) ||
961 (std::abs(dist_min) < std::numeric_limits<double>::epsilon() &&
962 dist_max > std::numeric_limits<double>::epsilon())
963
964 ) {
965
966 // Edges is on two sides of the surface
967 const double s = dist[0] / (dist[0] - dist[1]);
968 const double dist = s * ray_length;
969
970 auto add_edge = [&](auto dist) {
971 edgesToCut[e].dIst = dist;
972 edgesToCut[e].lEngth = ray_length;
973 edgesToCut[e].unitRayDir = ray;
974 edgesToCut[e].rayPoint = n0;
975 cutEdges.insert(e);
976
977 aveLength += norm_2(ray);
978 maxLength = fmax(maxLength, norm_2(ray));
979 ++nb_ave_length;
980 };
981
982 add_edge(dist);
983 }
984 }
985 aveLength /= nb_ave_length;
986
987 if (debug)
988 CHKERR SaveData(m_field.get_moab())("cut_edges.vtk", cutEdges);
989
990 if (debug)
991 CHKERR SaveData(m_field.get_moab())("cut_edges_vol.vtk",
992 get_adj(cutEdges, 3));
993
995}
constexpr double a
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:95
const double r
rate factor
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

Definition at line 1492 of file CutMeshInterface.cpp.

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

Definition at line 469 of file CutMeshInterface.hpp.

469{ return cutEdges; }

◆ getCutFrontVolumes()

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

Definition at line 486 of file CutMeshInterface.hpp.

486{ return cutFrontVolumes; }

◆ getCutSurfaceVolumes()

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

Definition at line 485 of file CutMeshInterface.hpp.

485{ return cutSurfaceVolumes; }

◆ getFront()

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

Definition at line 467 of file CutMeshInterface.hpp.

467{ return fRont; }

◆ getMergedSurfaces()

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

Definition at line 483 of file CutMeshInterface.hpp.

483{ return mergedSurfaces; }

◆ getMergedVolumes()

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

Definition at line 482 of file CutMeshInterface.hpp.

◆ getNewCutSurfaces()

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

Definition at line 471 of file CutMeshInterface.hpp.

471{ return cutNewSurfaces; }

◆ getNewCutVertices()

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

Definition at line 472 of file CutMeshInterface.hpp.

472{ return cutNewVertices; }

◆ getNewCutVolumes()

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

Definition at line 470 of file CutMeshInterface.hpp.

470{ return cutNewVolumes; }

◆ getNewTrimSurfaces()

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

Definition at line 479 of file CutMeshInterface.hpp.

◆ getNewTrimVertices()

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

Definition at line 480 of file CutMeshInterface.hpp.

480{ return trimNewVertices; }

◆ getNewTrimVolumes()

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

Definition at line 478 of file CutMeshInterface.hpp.

◆ getOptions()

MoFEMErrorCode MoFEM::CutMeshInterface::getOptions ( )

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 ( )

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
Examples
mesh_cut.cpp.

Definition at line 466 of file CutMeshInterface.hpp.

466{ return sUrface; }

◆ getTreeSurfPtr()

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

Definition at line 492 of file CutMeshInterface.hpp.

492 {
493 return treeSurfPtr;
494 }

◆ getTrimEdges()

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

Definition at line 477 of file CutMeshInterface.hpp.

477{ return trimEdges; }

◆ getVolume()

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

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 3074 of file CutMeshInterface.cpp.

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

◆ 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

Definition at line 2314 of file CutMeshInterface.cpp.

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

◆ 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 1459 of file CutMeshInterface.cpp.

1459 {
1460 CoreInterface &m_field = cOre;
1461 moab::Interface &moab = m_field.get_moab();
1463
1464 // Range out_side_vertices;
1465 for (auto m : verticesOnCutEdges) {
1466 double dist = m.second.dIst;
1467 VectorDouble3 new_coors = m.second.rayPoint + dist * m.second.unitRayDir;
1468 if (th)
1469 CHKERR moab.tag_set_data(th, &m.first, 1, &new_coors[0]);
1470 else
1471 CHKERR moab.set_coords(&m.first, 1, &new_coors[0]);
1472 }
1473
1475}

◆ moveMidNodesOnTrimmedEdges()

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

move trimmed edges mid nodes

Returns
error code

Definition at line 1477 of file CutMeshInterface.cpp.

1477 {
1478 CoreInterface &m_field = cOre;
1479 moab::Interface &moab = m_field.get_moab();
1481 for (auto &v : verticesOnTrimEdges) {
1482 double dist = v.second.dIst;
1483 VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1484 if (th)
1485 CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1486 else
1487 CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1488 }
1490}

◆ projectZeroDistanceEnts() [1/2]

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

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 997 of file CutMeshInterface.cpp.

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

2248 {
2249 CoreInterface &m_field = cOre;
2250 moab::Interface &moab = m_field.get_moab();
2251 PrismInterface *interface;
2253 CHKERR m_field.getInterface(interface);
2254 // Remove tris on surface front
2255 {
2256 Range front_tris;
2257 EntityHandle meshset;
2258 CHKERR moab.create_meshset(MESHSET_SET, meshset);
2259 CHKERR moab.add_entities(meshset, ents);
2260 CHKERR interface->findFacesWithThreeNodesOnInternalSurfaceSkin(
2261 meshset, split_bit, true, front_tris);
2262 CHKERR moab.delete_entities(&meshset, 1);
2263 ents = subtract(ents, front_tris);
2264 }
2265 // Remove entities on skin
2266 Range tets;
2267 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2268 split_bit, BitRefLevel().set(), MBTET, tets);
2269 // Remove entities on skin
2270 Skinner skin(&moab);
2271 Range tets_skin;
2272 CHKERR skin.find_skin(0, tets, false, tets_skin);
2273 ents = subtract(ents, tets_skin);
2274
2276}

◆ saveCutEdges()

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

Definition at line 3159 of file CutMeshInterface.cpp.

3159 {
3160 CoreInterface &m_field = cOre;
3161 moab::Interface &moab = m_field.get_moab();
3163 CHKERR SaveData(moab)(prefix + "out_vol.vtk", vOlume);
3164 CHKERR SaveData(moab)(prefix + "out_surface.vtk", sUrface);
3165 CHKERR SaveData(moab)(prefix + "out_cut_edges.vtk", cutEdges);
3166 CHKERR SaveData(moab)(prefix + "out_cut_new_volumes.vtk", cutNewVolumes);
3167 CHKERR SaveData(moab)(prefix + "out_cut_new_surfaces.vtk", cutNewSurfaces);
3168 CHKERR SaveData(moab)(prefix + "out_cut_zero_distance_ents.vtk",
3171}

◆ saveTrimEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveTrimEdges ( )

Definition at line 3173 of file CutMeshInterface.cpp.

3173 {
3174 moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
3176 CHKERR SaveData(moab)("out_trim_surface.vtk", sUrface);
3177 CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
3178 CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
3179 CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
3181}

◆ 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 3142 of file CutMeshInterface.cpp.

3143 {
3144 CoreInterface &m_field = cOre;
3145 moab::Interface &moab = m_field.get_moab();
3147 Range nodes;
3148 if (bit.none())
3149 CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3150 else
3151 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3152 bit, mask, MBVERTEX, nodes);
3153 std::vector<double> coords(3 * nodes.size());
3154 CHKERR moab.tag_get_data(th, nodes, &coords[0]);
3155 CHKERR moab.set_coords(nodes, &coords[0]);
3157}

◆ 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.

42 {
44 sUrface = surface;
46}

◆ 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 3126 of file CutMeshInterface.cpp.

3126 {
3127 CoreInterface &m_field = cOre;
3128 moab::Interface &moab = m_field.get_moab();
3130 Range nodes;
3131 if (bit.none())
3132 CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3133 else
3134 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3135 bit, BitRefLevel().set(), MBVERTEX, nodes);
3136 std::vector<double> coords(3 * nodes.size());
3137 CHKERR moab.get_coords(nodes, &coords[0]);
3138 CHKERR moab.tag_set_data(th, nodes, &coords[0]);
3140}

◆ 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:301

◆ 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 2278 of file CutMeshInterface.cpp.

2280 {
2281 CoreInterface &m_field = cOre;
2282 moab::Interface &moab = m_field.get_moab();
2283 PrismInterface *interface;
2285 CHKERR m_field.getInterface(interface);
2286 EntityHandle meshset_volume;
2287 CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
2288 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2289 split_bit, BitRefLevel().set(), MBTET, meshset_volume);
2290 EntityHandle meshset_surface;
2291 CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
2292 CHKERR moab.add_entities(meshset_surface, ents);
2293 CHKERR interface->getSides(meshset_surface, split_bit, true);
2294 CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
2295 true);
2296 CHKERR moab.delete_entities(&meshset_volume, 1);
2297 CHKERR moab.delete_entities(&meshset_surface, 1);
2298 if (th) {
2299 Range prisms;
2300 ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2301 bit, BitRefLevel().set(), MBPRISM, prisms);
2302 for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
2303 int num_nodes;
2304 const EntityHandle *conn;
2305 CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
2306 MatrixDouble data(3, 3);
2307 CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
2308 CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
2309 }
2310 }
2312}

◆ 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 1981 of file CutMeshInterface.cpp.

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

◆ 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 2040 of file CutMeshInterface.cpp.

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