v0.8.23
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 (const MOFEMuuid &uuid, 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 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_cut, const double tol_cut_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
 Cut mesh onlu. More...
 
MoFEMErrorCode trimOnly (const BitRefLevel trim_bit, Tag th, const double tol_cut_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_cut, 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_cut, 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 createLevelSets (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())
 Create a level sets, i.e. distances from surface and surface front. More...
 
MoFEMErrorCode createLevelSets (int verb=QUIET, const bool debug=false)
 Create a level sets, i.e. distances from surface and surface front. 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, Range *fixed_edges, Range *corner_nodes, const double geometry_tol, int verb=QUIET, const bool debug=false)
 find edges to cut More...
 
MoFEMErrorCode projectZeroDistanceEnts (Range *fixed_edges, Range *corner_nodes, 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 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 Range & getVolume () const
 
const Range & getSurface () const
 
const Range & getFront () const
 
const Range & getCutEdges () const
 
const Range & getNewCutVolumes () const
 
const Range & getNewCutSurfaces () const
 
const Range & getNewCutVertices () const
 
const Range & projectZeroDistanceEnts () const
 
const Range & getTrimEdges () const
 
const Range & getNewTrimVolumes () const
 
const Range & getNewTrimSurfaces () const
 
const Range & getNewTrimVertices () const
 
const Range & getMergedVolumes () const
 
const Range & getMergedSurfaces () const
 
const Range & getCutSurfaceVolumes () const
 
const Range & getCutFrontVolumes () const
 
void setTrimFixedEdges (const bool b)
 
MoFEMErrorCode saveCutEdges (const std::string prefix="")
 
MoFEMErrorCode saveTrimEdges ()
 
boost::shared_ptr< OrientedBoxTreeTool > & getTreeSurfPtr ()
 
MoFEMErrorCode clearMap ()
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (const MOFEMuuid &uuid, bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE , bool VERIFY = false>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, IFACE *&iface) const
 Get interface by uuid and return reference to pointer of 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 ()
 
virtual MoFEMErrorCode getLibVersion (Version &version) const
 Get library version. More...
 
virtual const MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version) const
 Get database major version. More...
 
virtual MoFEMErrorCode getInterfaceVersion (Version &version) const
 Get database major version. More...
 
template<>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, UnknownInterface *&iface) const
 

Public Attributes

MoFEM::CorecOre
 
int lineSearchSteps
 
int nbMaxMergingCycles
 
double projectEntitiesQualityTrashold
 

Private Attributes

Range fRont
 
Range sUrface
 
Range vOlume
 
bool trimFixedEdges
 
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
 

Additional Inherited Members

- Protected Member Functions inherited from MoFEM::UnknownInterface
boost::typeindex::type_index getClassIdx (const MOFEMuuid &uid) const
 Get type name for interface Id. More...
 
MOFEMuuid getUId (const boost::typeindex::type_index &class_idx) const
 Get interface Id for class name. More...
 

Detailed Description

Interface to cut meshes.

Examples
mesh_cut.cpp.

Definition at line 32 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 538 of file CutMeshInterface.hpp.

Constructor & Destructor Documentation

◆ CutMeshInterface()

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

Definition at line 33 of file CutMeshInterface.cpp.

34  : cOre(const_cast<Core &>(core)) {
35  lineSearchSteps = 10;
36  nbMaxMergingCycles = 200;
38  trimFixedEdges = true;
39 }

◆ ~CutMeshInterface()

MoFEM::CutMeshInterface::~CutMeshInterface ( )

Definition at line 39 of file CutMeshInterface.hpp.

39 {}

Member Function Documentation

◆ buildTree()

MoFEMErrorCode MoFEM::CutMeshInterface::buildTree ( )

build tree

Returns
error code
Examples
mesh_cut.cpp.

Definition at line 218 of file CutMeshInterface.cpp.

218  {
219  CoreInterface &m_field = cOre;
220  moab::Interface &moab = m_field.get_moab();
222  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
223  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
226 }
moab::Interface & get_moab()
Definition: Core.hpp:266
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ clearMap()

MoFEMErrorCode MoFEM::CutMeshInterface::clearMap ( )

Definition at line 41 of file CutMeshInterface.cpp.

41  {
43  treeSurfPtr.reset();
45 }
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

62  {
63  CoreInterface &m_field = cOre;
64  moab::Interface &moab = m_field.get_moab();
66  sUrface.clear();
67  std::map<EntityHandle, EntityHandle> verts_map;
68  for (Range::const_iterator tit = surface.begin(); tit != surface.end();
69  tit++) {
70  int num_nodes;
71  const EntityHandle *conn;
72  CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
73  MatrixDouble coords(num_nodes, 3);
74  if (th) {
75  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords(0, 0));
76  } else {
77  CHKERR moab.get_coords(conn, num_nodes, &coords(0, 0));
78  }
79  EntityHandle new_verts[num_nodes];
80  for (int nn = 0; nn != num_nodes; nn++) {
81  if (verts_map.find(conn[nn]) != verts_map.end()) {
82  new_verts[nn] = verts_map[conn[nn]];
83  } else {
84  if (transform) {
85  ublas::matrix_row<MatrixDouble> mr(coords, nn);
86  if (origin) {
87  VectorAdaptor vec_origin(
88  3, ublas::shallow_array_adaptor<double>(3, origin));
89  mr = mr - vec_origin;
90  }
91  MatrixAdaptor mat_transform = MatrixAdaptor(
92  3, 3, ublas::shallow_array_adaptor<double>(9, transform));
93  mr = prod(mat_transform, mr);
94  if (origin) {
95  VectorAdaptor vec_origin(
96  3, ublas::shallow_array_adaptor<double>(3, origin));
97  mr = mr + vec_origin;
98  }
99  }
100  if (shift) {
101  ublas::matrix_row<MatrixDouble> mr(coords, nn);
102  VectorAdaptor vec_shift(
103  3, ublas::shallow_array_adaptor<double>(3, shift));
104  mr = mr + vec_shift;
105  }
106  CHKERR moab.create_vertex(&coords(nn, 0), new_verts[nn]);
107  verts_map[conn[nn]] = new_verts[nn];
108  }
109  }
110  EntityHandle ele;
111  CHKERR moab.create_element(MBTRI, new_verts, num_nodes, ele);
112  sUrface.insert(ele);
113  }
114  if (!save_mesh.empty())
115  CHKERR SaveData(m_field.get_moab())(save_mesh, sUrface);
117 }
moab::Interface & get_moab()
Definition: Core.hpp:266
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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)
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:118
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:101
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

545  {
546  CoreInterface &m_field = cOre;
547  moab::Interface &moab = m_field.get_moab();
549 
550  auto create_tag = [&](const std::string name, const int dim) {
551  Tag th;
552  rval = moab.tag_get_handle(name.c_str(), th);
553  if (rval == MB_SUCCESS)
554  return th;
555  std::vector<double> def_val(dim, 0);
556  CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
557  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
558  return th;
559  };
560 
561  Range vol_vertices;
562  CHKERR moab.get_connectivity(vol, vol_vertices, true);
563 
564  std::vector<double> min_distances_from_front(vol_vertices.size(), -1);
565  std::vector<double> points_on_edges(3 * vol_vertices.size(), 0);
566  std::vector<EntityHandle> closest_edges(vol_vertices.size(), 0);
567 
568  std::vector<double> coords(3 * vol_vertices.size());
569  if (th)
570  CHKERR moab.tag_get_data(th, vol_vertices, &*coords.begin());
571  else
572  CHKERR moab.get_coords(vol_vertices, &*coords.begin());
573 
574  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
575  &*coords.begin(), vol_vertices.size(), fRont,
576  &*min_distances_from_front.begin(), &*points_on_edges.begin(),
577  &*closest_edges.begin());
578 
579  if (!points_on_edges.empty()) {
580  for (int i = 0; i != min_distances_from_front.size(); ++i) {
581  Range faces;
582  CHKERR moab.get_adjacencies(&closest_edges[0], 1, 2, false, faces);
583  auto point_in = getVectorAdaptor(&coords[3 * i], 3);
584  auto point_out = getVectorAdaptor(&points_on_edges[3 * i], 3);
585  point_out -= point_in;
586  }
587  }
588 
589  auto th_dist_front_vec = create_tag("DIST_FRONT_VECTOR", 3);
590  CHKERR moab.tag_set_data(th_dist_front_vec, vol_vertices,
591  &*points_on_edges.begin());
592 
594 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:43
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ createLevelSets() [1/2]

MoFEMErrorCode MoFEM::CutMeshInterface::createLevelSets ( 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() 
)

Create a level sets, i.e. distances from surface and surface front.

Parameters
th
vol_edges
remove_adj_prims_edges
verb
debug
edges_file_name
Returns
MoFEMErrorCode

Definition at line 596 of file CutMeshInterface.cpp.

598  {
599  CoreInterface &m_field = cOre;
600  moab::Interface &moab = m_field.get_moab();
602 
603  auto get_tag_data = [&](auto th, auto conn) {
604  const void *ptr;
605  CHKERR moab.tag_get_by_ptr(th, &conn, 1, &ptr);
606  return getVectorAdaptor(
607  const_cast<double *>(static_cast<const double *>(ptr)), 3);
608  };
609 
610  auto get_edge_ray = [&](auto conn) {
611  VectorDouble6 coords(6);
612  CHKERR moab.get_coords(conn, 2, &coords[0]);
613  VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
614  VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
615  VectorDouble3 ray = n1 - n0;
616  return ray;
617  };
618 
619  Range cut_edges;
620 
621  Range edges;
622  CHKERR moab.get_adjacencies(vOlume, 1, true, edges, moab::Interface::UNION);
623 
624  auto remove_prisms_edges = [&](Range &edges) {
626  Range prisms;
627  CHKERR moab.get_adjacencies(edges, 3, false, prisms,
628  moab::Interface::UNION);
629  prisms = prisms.subset_by_type(MBPRISM);
630  Range prisms_verts;
631  CHKERR moab.get_connectivity(prisms, prisms_verts, true);
632  Range prism_edges;
633  CHKERR moab.get_adjacencies(prisms_verts, 1, false, prism_edges,
634  moab::Interface::UNION);
635  edges = subtract(edges, prism_edges);
637  };
638  if (remove_adj_prims_edges)
639  CHKERR remove_prisms_edges(edges);
640 
641  for (auto e : edges) {
642 
643  int num_nodes;
644  const EntityHandle *conn;
645  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
646  auto ray = get_edge_ray(conn);
647  const double length = norm_2(ray);
648  ray /= length;
649 
650  auto signed_norm = [&](const auto &v) { return inner_prod(ray, v); };
651 
652  auto get_cut_edges = [&](auto th, Range &cut_edges) {
654  const auto dist0 = get_tag_data(th, conn[0]);
655  const auto dist1 = get_tag_data(th, conn[1]);
656  const double min_dist = std::min(norm_2(dist0), norm_2(dist1));
657  if (min_dist < length) {
658  auto opposite = inner_prod(dist0, dist1);
659  if (opposite <= 0) {
660  const double sign_dist0 = signed_norm(dist0);
661  const double sign_dist1 = signed_norm(dist1);
662  if (sign_dist0 > -std::numeric_limits<double>::epsilon() &&
663  sign_dist1 < std::numeric_limits<double>::epsilon())
664  cut_edges.insert(e);
665  }
666  }
668  };
669 
670  CHKERR get_cut_edges(th, cut_edges);
671  }
672 
673  CHKERR moab.get_adjacencies(cut_edges, 3, false, vol_edges,
674  moab::Interface::UNION);
675 
676  vol_edges = intersect(vol_edges, vOlume);
677 
678  if (debug && !edges_file_name.empty())
679  CHKERR SaveData(m_field.get_moab())(edges_file_name, cut_edges);
680 
682 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:88
static const bool debug
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:43
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:101
#define CHKERR
Inline error check.
Definition: definitions.h:596
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ createLevelSets() [2/2]

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

Create a level sets, i.e. distances from surface and surface front.

Parameters
update_front
verb
debug
Returns
MoFEMErrorCode

Definition at line 684 of file CutMeshInterface.cpp.

684  {
685  CoreInterface &m_field = cOre;
686  moab::Interface &moab = m_field.get_moab();
688 
689  CHKERR createFrontLevelSets(vOlume, nullptr, verb, debug);
690  Tag th_dist_front_vec;
691  CHKERR moab.tag_get_handle("DIST_FRONT_VECTOR", th_dist_front_vec);
692  CHKERR createLevelSets(th_dist_front_vec, cutFrontVolumes, true, verb, debug,
693  "cutFrontEdges.vtk");
694 
696 
697  Tag th_dist_surface_vec;
698  CHKERR moab.tag_get_handle("DIST_SURFACE_VECTOR", th_dist_surface_vec);
699  cutSurfaceVolumes.clear();
700  CHKERR createLevelSets(th_dist_surface_vec, cutSurfaceVolumes, true, verb,
701  debug, "cutSurfaceEdges.vtk");
702 
703  if (debug)
704  CHKERR SaveData(m_field.get_moab())("level_sets.vtk", vOlume);
705  if (debug)
706  CHKERR SaveData(m_field.get_moab())("cutSurfaceVolumes.vtk",
708  if (debug)
709  CHKERR SaveData(m_field.get_moab())("cutFrontVolumes.vtk", cutFrontVolumes);
710 
712 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode createSurfaceLevelSets(int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to cut surface.
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
MoFEMErrorCode createLevelSets(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())
Create a level sets, i.e. distances from surface and surface front.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
MoFEMErrorCode createFrontLevelSets(Range vol, Tag th=nullptr, int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to surface front.

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

464  {
465  CoreInterface &m_field = cOre;
466  moab::Interface &moab = m_field.get_moab();
468  auto tools_interface = m_field.getInterface<Tools>();
469 
470  auto create_tag = [&](const std::string name, const int dim) {
471  Tag th;
472  rval = moab.tag_get_handle(name.c_str(), th);
473  if (rval == MB_SUCCESS)
474  return th;
475  std::vector<double> def_val(dim, 0);
476  CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
477  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
478 
479  return th;
480  };
481 
482  auto set_vol = [&](const Range &vol_verts, std::vector<double> &coords,
483  std::vector<double> &dist_surface_vec,
484  std::vector<double> &dist_surface_normal_vec) {
486 
487  coords.resize(3 * vol_verts.size());
488  dist_surface_vec.resize(3 * vol_verts.size());
489  dist_surface_normal_vec.resize(3 * vol_verts.size());
490  CHKERR moab.get_coords(vol_verts, &*coords.begin());
491  std::srand(0);
492 
493  for (auto v : vol_verts) {
494 
495  const int index = vol_verts.index(v);
496  auto point_in = getVectorAdaptor(&coords[3 * index], 3);
497  VectorDouble3 point_out(3);
498  EntityHandle facets_out;
499  CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
500  &point_out[0], facets_out);
501 
502  VectorDouble3 n(3);
503  CHKERR tools_interface->getTriNormal(facets_out, &*n.begin());
504  n /= norm_2(n);
505 
506  VectorDouble3 delta = point_out - point_in;
507  if (norm_2(delta) < std::numeric_limits<double>::epsilon()) {
508  if (std::rand() % 2 == 0)
509  delta += n * std::numeric_limits<double>::epsilon();
510  else
511  delta -= n * std::numeric_limits<double>::epsilon();
512  }
513 
514  auto dist_vec = getVectorAdaptor(&dist_surface_vec[3 * index], 3);
515  noalias(dist_vec) = delta;
516 
517  auto dist_normal_vec =
518  getVectorAdaptor(&dist_surface_normal_vec[3 * index], 3);
519  noalias(dist_normal_vec) = inner_prod(delta, n) * n;
520  }
521 
523  };
524 
525  std::vector<double> coords;
526  std::vector<double> dist_surface_vec;
527  std::vector<double> dist_surface_normal_vec;
528  Range vol_verts;
529  CHKERR moab.get_connectivity(vOlume, vol_verts, true);
530 
531  CHKERR set_vol(vol_verts, coords, dist_surface_vec, dist_surface_normal_vec);
532 
533  auto th_dist_surface_vec = create_tag("DIST_SURFACE_VECTOR", 3);
534  auto th_dist_surface_normal_vec = create_tag("DIST_SURFACE_NORMAL_VECTOR", 3);
535  CHKERR moab.tag_set_data(th_dist_surface_vec, vol_verts,
536  &*dist_surface_vec.begin());
537  CHKERR moab.tag_set_data(th_dist_surface_normal_vec, vol_verts,
538  &*dist_surface_normal_vec.begin());
539 
541 }
moab::Interface & get_moab()
Definition: Core.hpp:266
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:43
#define CHKERR
Inline error check.
Definition: definitions.h:596
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ cutAndTrim()

MoFEMErrorCode MoFEM::CutMeshInterface::cutAndTrim ( int first_bit,
Tag  th,
const double  tol_cut,
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_cut
tol_cut_close
tol_trim_close
fixed_edges
corner_nodes
update_meshsets
debug
Returns
MoFEMErrorCode

Definition at line 306 of file CutMeshInterface.cpp.

309  {
310  CoreInterface &m_field = cOre;
312 
313  std::vector<BitRefLevel> bit_levels;
314 
315  auto get_back_bit_levels = [&]() {
316  bit_levels.push_back(BitRefLevel());
317  bit_levels.back().set(first_bit + bit_levels.size() - 1);
318  return bit_levels.back();
319  };
320 
321  BitRefLevel cut_bit = get_back_bit_levels();
322 
323  CHKERR cutOnly(unite(cutSurfaceVolumes, cutFrontVolumes), cut_bit, th,
324  tol_cut, tol_cut_close, fixed_edges, corner_nodes,
325  update_meshsets, debug);
326 
327  auto get_min_quality = [&m_field](const BitRefLevel bit, Tag th) {
328  Range tets_level; // test at level
329  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
330  bit, BitRefLevel().set(), MBTET, tets_level);
331  double min_q = 1;
332  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
333  return min_q;
334  };
335 
336  PetscPrintf(PETSC_COMM_WORLD, "Min quality cut %6.4g\n",
337  get_min_quality(cut_bit, th));
338 
339  Range starting_volume = cutNewVolumes;
340  Range starting_volume_nodes;
341  CHKERR m_field.get_moab().get_connectivity(starting_volume,
342  starting_volume_nodes, true);
343  std::vector<double> staring_volume_coords(3 * starting_volume_nodes.size());
344  CHKERR m_field.get_moab().get_coords(starting_volume_nodes,
345  &*staring_volume_coords.begin());
346 
347  if (th) {
348  std::vector<double> staring_volume_th_coords(3 *
349  starting_volume_nodes.size());
350  CHKERR m_field.get_moab().tag_get_data(th, starting_volume_nodes,
351  &*staring_volume_th_coords.begin());
352  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
353  &*staring_volume_th_coords.begin());
354  }
355 
356  if (th)
357  CHKERR setTagData(th);
358 
359  BitRefLevel trim_bit = get_back_bit_levels();
360 
361  CHKERR trimOnly(trim_bit, th, tol_trim_close, fixed_edges, corner_nodes,
362  update_meshsets, debug);
363 
364  PetscPrintf(PETSC_COMM_WORLD, "Min quality trim %3.2g\n",
365  get_min_quality(trim_bit, th));
366 
367  first_bit += bit_levels.size() - 1;
368 
369  if (th)
370  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
371  &*staring_volume_coords.begin());
372 
374 }
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode trimOnly(const BitRefLevel trim_bit, Tag th, const double tol_cut_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Trim mesh only.
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
MoFEMErrorCode cutOnly(Range vol, const BitRefLevel cut_bit, Tag th, const double tol_cut, const double tol_cut_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Cut mesh onlu.

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

1366  {
1367  CoreInterface &m_field = cOre;
1368  moab::Interface &moab = m_field.get_moab();
1369  MeshRefinement *refiner;
1370  const RefEntity_multiIndex *ref_ents_ptr;
1372 
1373  if (cutEdges.size() != edgesToCut.size())
1374  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1375 
1376  CHKERR m_field.getInterface(refiner);
1377  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1378  CHKERR refiner->add_vertices_in_the_middle_of_edges(cutEdges, bit);
1379  CHKERR refiner->refine_TET(vOlume, bit, false, QUIET,
1380  debug ? &cutEdges : NULL);
1381 
1382  if (debug)
1383  if (cutEdges.size() != edgesToCut.size())
1384  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1385 
1386  cut_vols.clear();
1387  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1388  bit, BitRefLevel().set(), MBTET, cut_vols);
1389  cut_surf.clear();
1390  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1391  bit, BitRefLevel().set(), MBTRI, cut_surf);
1392 
1393  // Find new vertices on cut edges
1394  cut_verts.clear();
1395  CHKERR moab.get_connectivity(zeroDistanceEnts, cut_verts, true);
1396  cut_verts.merge(zeroDistanceVerts);
1397  for (auto &m : edgesToCut) {
1398  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1399  boost::make_tuple(MBVERTEX, m.first));
1400  if (vit ==
1401  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1402  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1403  "No vertex on cut edges, that make no sense");
1404  }
1405  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1406  if ((ref_ent->getBitRefLevel() & bit).any()) {
1407  EntityHandle vert = ref_ent->getRefEnt();
1408  cut_verts.insert(vert);
1409  verticesOnCutEdges[vert] = m.second;
1410  } else {
1411  SETERRQ1(
1412  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1413  "Vertex has wrong bit ref level %s",
1414  boost::lexical_cast<std::string>(ref_ent->getBitRefLevel()).c_str());
1415  }
1416  }
1417 
1418  // Add zero distance entities faces
1419  Range tets_skin;
1420  Skinner skin(&moab);
1421  CHKERR skin.find_skin(0, cut_vols, false, tets_skin);
1422  cut_surf.merge(zeroDistanceEnts.subset_by_type(MBTRI));
1423 
1424  // At that point cut_surf has all newly created faces, now take all
1425  // nodes on those faces and subtract nodes on cut edges. Faces adjacent to
1426  // nodes which left are not part of surface.
1427  Range diff_verts;
1428  CHKERR moab.get_connectivity(unite(cut_surf, zeroDistanceEnts), diff_verts,
1429  true);
1430  diff_verts = subtract(diff_verts, cut_verts);
1431  Range subtract_faces;
1432  CHKERR moab.get_adjacencies(diff_verts, 2, false, subtract_faces,
1433  moab::Interface::UNION);
1434  cut_surf = subtract(cut_surf, unite(subtract_faces, tets_skin));
1435  cut_verts.clear();
1436  CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1437 
1438  // Check non-mainfolds
1439  auto check_for_non_minfold = [&]() {
1441  Range surf_edges;
1442  CHKERR moab.get_adjacencies(cut_surf, 1, false, surf_edges,
1443  moab::Interface::UNION);
1444  for (auto e : surf_edges) {
1445 
1446  Range faces;
1447  CHKERR moab.get_adjacencies(&e, 1, 2, false, faces);
1448  faces = intersect(faces, cut_surf);
1449  if (faces.size() > 2) {
1450 
1451  bool resolved = false;
1452 
1453  // Check for haning node
1454  Range nodes;
1455  CHKERR moab.get_connectivity(faces, nodes, true);
1456  for (auto n : nodes) {
1457  Range adj_faces;
1458  CHKERR moab.get_adjacencies(&n, 1, 2, false, adj_faces);
1459  adj_faces = intersect(adj_faces, cut_surf);
1460  if (adj_faces.size() == 1) {
1461  cut_surf.erase(adj_faces[0]);
1462  resolved = true;
1463  }
1464  }
1465 
1466  // Check for two edges minfold
1467  Range adj_edges;
1468  CHKERR moab.get_adjacencies(faces, 1, false, adj_edges,
1469  moab::Interface::UNION);
1470  adj_edges = intersect(adj_edges, surf_edges);
1471  adj_edges.erase(e);
1472  for (auto other_e : adj_edges) {
1473  Range other_faces;
1474  CHKERR moab.get_adjacencies(&other_e, 1, 2, false, other_faces);
1475  other_faces = intersect(other_faces, cut_surf);
1476  if (other_faces.size() > 2) {
1477  other_faces = intersect(other_faces, faces);
1478  cut_surf = subtract(cut_surf, other_faces);
1479  resolved = true;
1480  }
1481  }
1482 
1483  if (!resolved && !debug)
1484  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1485  "Non-mainfold surface");
1486 
1487  cut_verts.clear();
1488  CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1489  }
1490  }
1492  };
1493 
1494  CHKERR check_for_non_minfold();
1495 
1497 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
map< EntityHandle, TreeData > edgesToCut

◆ cutOnly()

MoFEMErrorCode MoFEM::CutMeshInterface::cutOnly ( Range  vol,
const BitRefLevel  cut_bit,
Tag  th,
const double  tol_cut,
const double  tol_cut_close,
Range *  fixed_edges = NULL,
Range *  corner_nodes = NULL,
const bool  update_meshsets = false,
const bool  debug = false 
)

Cut mesh onlu.

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

Definition at line 229 of file CutMeshInterface.cpp.

232  {
233  CoreInterface &m_field = cOre;
234  moab::Interface &moab = m_field.get_moab();
236 
237  // cut mesh
238  CHKERR findEdgesToCut(vol, fixed_edges, corner_nodes, tol_cut, QUIET, debug);
239  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_cut_close,
240  QUIET, debug);
243  if (fixed_edges)
244  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
245  *fixed_edges);
246  if (corner_nodes)
247  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*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 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode cutEdgesInMiddle(const BitRefLevel bit, Range &cut_vols, Range &cut_surf, Range &cut_verts, const bool debug=false)
cut edges
const Range & projectZeroDistanceEnts() const
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode saveCutEdges(const std::string prefix="")
static const bool debug
MoFEMErrorCode moveMidNodesOnCutEdges(Tag th=NULL)
projecting of mid edge nodes on new mesh on surface
#define CHKERR
Inline error check.
Definition: definitions.h:596
MoFEMErrorCode findEdgesToCut(Range vol, Range *fixed_edges, Range *corner_nodes, const double geometry_tol, int verb=QUIET, const bool debug=false)
find edges to cut
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ cutTrimAndMerge()

MoFEMErrorCode MoFEM::CutMeshInterface::cutTrimAndMerge ( int first_bit,
const int  fraction_level,
Tag  th,
const double  tol_cut,
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_cuttolerance 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 376 of file CutMeshInterface.cpp.

379  {
380  CoreInterface &m_field = cOre;
382 
383  std::vector<BitRefLevel> bit_levels;
384 
385  auto get_back_bit_levels = [&]() {
386  bit_levels.push_back(BitRefLevel());
387  bit_levels.back().set(first_bit + bit_levels.size() - 1);
388  return bit_levels.back();
389  };
390 
391  if (debug) {
392  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
393  "ents_not_in_database.vtk", "VTK", "");
394  }
395 
396  CHKERR cutAndTrim(first_bit, th, tol_cut, tol_cut_close, tol_trim_close,
397  &fixed_edges, &corner_nodes, update_meshsets, debug);
398  if (debug)
399  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
400  "cut_trim_ents_not_in_database.vtk", "VTK", "");
401 
402  BitRefLevel bit_level1 = BitRefLevel().set(first_bit - 1);
403  BitRefLevel bit_level2 = get_back_bit_levels();
404  BitRefLevel bit_level3 = get_back_bit_levels();
405 
406  CHKERR mergeBadEdges(fraction_level, bit_level1, bit_level2, bit_level3,
407  getNewTrimSurfaces(), fixed_edges, corner_nodes, th,
408  update_meshsets, debug);
409 
410  auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
411  Range tets_level; // test at level
412  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
413  bit, BitRefLevel().set(), MBTET, tets_level);
414  double min_q = 1;
415  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
416  if (min_q < 0 && debug) {
417  CHKERR m_field.getInterface<Tools>()->writeTetsWithQuality(
418  "negative_tets.vtk", "VTK", "", tets_level, th);
419  }
420  return min_q;
421  };
422 
423  PetscPrintf(PETSC_COMM_WORLD, "Min quality node merge %6.4g\n",
424  get_min_quality(bit_level3, th));
425 
426  CHKERR cOre.getInterface<BitRefManager>()->updateRange(fixed_edges,
427  fixed_edges);
428  CHKERR cOre.getInterface<BitRefManager>()->updateRange(corner_nodes,
429  corner_nodes);
430 
431  first_bit += bit_levels.size() - 1;
432 
433  if (debug) {
434  CHKERR cOre.getInterface<BitRefManager>()->writeBitLevelByType(
435  bit_level3, BitRefLevel().set(), MBTET, "out_tets_merged.vtk", "VTK",
436  "");
437  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
438  "cut_trim_merge_ents_not_in_database.vtk", "VTK", "");
439  CHKERR SaveData(m_field.get_moab())("merged_surface.vtk", mergedSurfaces);
440  }
441 
443 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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 getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
const Range & getNewTrimSurfaces() const
MoFEMErrorCode cutAndTrim(int &first_bit, Tag th, const double tol_cut, 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.
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ findEdgesToCut()

MoFEMErrorCode MoFEM::CutMeshInterface::findEdgesToCut ( Range  vol,
Range *  fixed_edges,
Range *  corner_nodes,
const double  geometry_tol,
int  verb = QUIET,
const bool  debug = false 
)

find edges to cut

Parameters
volis tetrahedrons search to cut
fixed_edgespointer to fixed edges
corner_nodespointer to corner edges
geometry_toltolerance for close geometry fetchers
verbverbosity level
debugdebugging
Returns
MoFEMErrorCode

Definition at line 801 of file CutMeshInterface.cpp.

804  {
805  CoreInterface &m_field = cOre;
806  moab::Interface &moab = m_field.get_moab();
808 
809  edgesToCut.clear();
810  cutEdges.clear();
811 
812  zeroDistanceVerts.clear();
813  zeroDistanceEnts.clear();
814  verticesOnCutEdges.clear();
815 
816  Tag th_test_dist;
817  if (debug) {
818  double def_val[] = {0, 0, 0};
819  CHKERR moab.tag_get_handle("TETS_DIST", 1, MB_TYPE_DOUBLE, th_test_dist,
820  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
821  }
822 
823  Tag th_dist_normal;
824  CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_VECTOR", th_dist_normal);
825 
826  auto get_tag_data = [&](auto th, auto conn) {
827  const void *ptr;
828  CHKERR moab.tag_get_by_ptr(th, &conn, 1, &ptr);
829  return getVectorAdaptor(
830  const_cast<double *>(static_cast<const double *>(ptr)), 3);
831  };
832 
833  Range vol_edges;
834  CHKERR moab.get_adjacencies(vol, 1, true, vol_edges, moab::Interface::UNION);
835 
836  aveLength = 0;
837  maxLength = 0;
838  int nb_ave_length = 0;
839  for (auto e : vol_edges) {
840 
841  int num_nodes;
842  const EntityHandle *conn;
843  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
844 
845  VectorDouble6 coords(6);
846  CHKERR moab.get_coords(conn, 2, &coords[0]);
847  VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
848  VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
849  VectorDouble3 ray = n1 - n0;
850  const double ray_length = norm_2(ray);
851  ray /= ray_length;
852 
853  auto dist_vec0 = get_tag_data(th_dist_normal, conn[0]);
854  auto dist_vec1 = get_tag_data(th_dist_normal, conn[1]);
855 
856  const double s0 = norm_2(dist_vec0);
857  const double s1 = norm_2(dist_vec1);
858 
859  if (inner_prod(dist_vec0, dist_vec1) < 0) {
860 
861  // Edges is on two sides of the surface
862  const double s = s0 / (s0 + s1);
863  const double dist = s * ray_length;
864 
865  VectorDouble3 p = n0 + dist * ray;
866  VectorDouble3 w = n0 + dist_vec0;
867  VectorDouble3 v = n1 + dist_vec1;
868  double t;
869  auto res = Tools::minDistancePointFromOnSegment(&w[0], &v[0], &p[0], &t);
870  t = std::max(0., std::min(t, 1.));
871  double d = 0;
872  if (res == Tools::SOLUTION_EXIST) {
873  VectorDouble3 o = w + t * (v - w);
874  d = norm_2(o - p) / ray_length;
875  }
876 
877  if (debug)
878  CHKERR moab.tag_set_data(th_test_dist, &e, 1, &d);
879 
880  if (d < 0.25) {
881 
882  edgesToCut[e].dIst = dist;
883  edgesToCut[e].lEngth = ray_length;
884  edgesToCut[e].unitRayDir = ray;
885  edgesToCut[e].rayPoint = n0;
886  cutEdges.insert(e);
887 
888  aveLength += norm_2(ray);
889  maxLength = fmax(maxLength, norm_2(ray));
890  ++nb_ave_length;
891  }
892  }
893  }
894  aveLength /= nb_ave_length;
895 
896  auto not_project_node = [this, &moab](const EntityHandle v) {
898  VectorDouble3 s0(3);
899  CHKERR moab.get_coords(&v, 1, &s0[0]);
900  verticesOnCutEdges[v].dIst = 0;
901  verticesOnCutEdges[v].lEngth = 0;
902  verticesOnCutEdges[v].unitRayDir.resize(3, false);
903  verticesOnCutEdges[v].unitRayDir.clear();
904  verticesOnCutEdges[v].rayPoint = s0;
906  };
907 
908  auto project_node = [this, &moab](const EntityHandle v,
909  VectorDouble3 dist_normal) {
911  VectorDouble3 s0(3);
912  CHKERR moab.get_coords(&v, 1, &s0[0]);
913  verticesOnCutEdges[v].dIst = 1;
914  verticesOnCutEdges[v].lEngth = norm_2(dist_normal);
915  verticesOnCutEdges[v].unitRayDir = dist_normal;
916  verticesOnCutEdges[v].rayPoint = s0;
918  };
919 
920  auto get_ave_edge_length = [&](const EntityHandle ent,
921  const Range &vol_edges) {
922 
923  Range adj_edges;
924  if (moab.type_from_handle(ent) == MBVERTEX)
925  CHKERR moab.get_adjacencies(&ent, 1, 1, false, adj_edges);
926  else
927  adj_edges.insert(ent);
928  adj_edges = intersect(adj_edges, vol_edges);
929 
930  double ave_l = 0;
931  for (auto e : adj_edges) {
932  int num_nodes;
933  const EntityHandle *conn;
934  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
935  VectorDouble6 coords(6);
936  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
938  &coords[0], &coords[1], &coords[2]);
940  &coords[3], &coords[4], &coords[5]);
942  t_n0(i) -= t_n1(i);
943  const double l = sqrt(t_n0(i) * t_n0(i));
944  ave_l += l;
945  }
946  return ave_l / adj_edges.size();
947  };
948 
949  auto project_vertices_close_to_geometry_features = [&]() {
951 
952  Range vol_vertices;
953  CHKERR moab.get_connectivity(vol, vol_vertices, true);
954 
955  Range fixed_edges_verts;
956  if (fixed_edges)
957  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts, true);
958  if (corner_nodes)
959  fixed_edges_verts.merge(*corner_nodes);
960 
961  Range fix_vertices = intersect(fixed_edges_verts, vol_vertices);
962 
963  for (auto v : fix_vertices) {
964 
965  VectorDouble3 dist_normal(3);
966  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, &*dist_normal.begin());
967  const double dist = norm_2(dist_normal);
968 
969  const double tol = get_ave_edge_length(v, vol_edges) * geometry_tol;
970  if (dist < tol) {
971  CHKERR not_project_node(v);
972  zeroDistanceVerts.insert(v);
973  }
974  }
975 
976  Skinner skin(&moab);
977  Range tets_skin;
978  CHKERR skin.find_skin(0, vOlume, false, tets_skin);
979  Range tets_skin_verts;
980  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
981 
982  for (auto v : subtract(tets_skin_verts, fix_vertices)) {
983 
984  VectorDouble3 dist_normal(3);
985  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, &*dist_normal.begin());
986  const double dist = norm_2(dist_normal);
987 
988  const double tol =
989  get_ave_edge_length(v, vol_edges) * pow(geometry_tol, 2);
990  if (dist < tol) {
991  CHKERR not_project_node(v);
992  zeroDistanceVerts.insert(v);
993  }
994  }
995 
996  for (auto v : subtract(vol_vertices, tets_skin_verts)) {
997 
998  VectorDouble3 dist_normal(3);
999  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, &*dist_normal.begin());
1000  const double dist = norm_2(dist_normal);
1001 
1002  const double tol =
1003  get_ave_edge_length(v, vol_edges) * pow(geometry_tol, 3);
1004  if (dist < tol) {
1005  CHKERR project_node(v, dist_normal);
1006  zeroDistanceVerts.insert(v);
1007  }
1008  }
1009 
1011  };
1012 
1013  CHKERR project_vertices_close_to_geometry_features();
1014 
1015  for (auto v : zeroDistanceVerts) {
1016  Range adj_edges;
1017  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges);
1018  for (auto e : adj_edges) {
1019  cutEdges.erase(e);
1020  edgesToCut.erase(e);
1021  }
1022  }
1023 
1024  if (debug)
1025  CHKERR SaveData(m_field.get_moab())("cut_edges.vtk", cutEdges);
1026 
1027  if (debug)
1028  CHKERR SaveData(m_field.get_moab())("cut_edges_zero_distance_verts.vtk",
1030 
1032 }
double maxLength
Maximal edge length.
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
double aveLength
Average edge length.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:88
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
double tol
static const bool debug
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:43
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:101
#define CHKERR
Inline error check.
Definition: definitions.h:596
static SEGMENT_MIN_DISTANCE minDistancePointFromOnSegment(const double *w_ptr, const double *v_ptr, const double *p_ptr, double *const t_ptr=nullptr)
Find closet point on the segment from the point.
Definition: Tools.cpp:319
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
map< EntityHandle, TreeData > edgesToCut

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

1537  {
1538  CoreInterface &m_field = cOre;
1539  moab::Interface &moab = m_field.get_moab();
1541 
1542  // takes body skin
1543  Skinner skin(&moab);
1544  Range tets_skin;
1545  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
1546 
1547  // vertices on the skin
1548  Range tets_skin_verts;
1549  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1550  // edges on the skin
1551  Range tets_skin_edges;
1552  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1553  moab::Interface::UNION);
1554  // get edges on new surface
1555  Range cut_surface_edges;
1556  CHKERR moab.get_adjacencies(cutNewSurfaces, 1, false, cut_surface_edges,
1557  moab::Interface::UNION);
1558  Range cut_surface_verts;
1559  CHKERR moab.get_connectivity(cutNewSurfaces, cut_surface_verts, true);
1560 
1561  Range corners;
1562  if (corner_nodes)
1563  corners = *corner_nodes;
1564 
1565  Range fix_edges;
1566  if (fixed_edges)
1567  fix_edges = *fixed_edges;
1568 
1569  Range fixed_edges_verts;
1570  if (fixed_edges)
1571  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts, true);
1572 
1573  Range surface_skin;
1574  if (fRont.empty())
1575  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
1576  else
1577  surface_skin = fRont;
1578 
1579  auto get_point_coords = [&](EntityHandle v) {
1580  VectorDouble3 point(3);
1581  if (th)
1582  CHKERR moab.tag_get_data(th, &v, 1, &point[0]);
1583  else
1584  CHKERR moab.get_coords(&v, 1, &point[0]);
1585  return point;
1586  };
1587 
1588  struct VertMap {
1589  const double d;
1590  const EntityHandle v;
1591  const EntityHandle e;
1592  VertMap(const double d, const EntityHandle v, const EntityHandle e)
1593  : d(d), v(v), e(e) {}
1594  };
1595 
1596  typedef multi_index_container<
1597  VertMap,
1598  indexed_by<
1599  ordered_non_unique<member<VertMap, const double, &VertMap::d>>,
1600  ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::v>>,
1601  ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::e>>
1602 
1603  >>
1604  VertMapMultiIndex;
1605 
1606  VertMapMultiIndex verts_map;
1607 
1608  auto add_vert = [&](const EntityHandle v, const EntityHandle e,
1609  const double dist) {
1610  verts_map.insert(VertMap(dist, v, e));
1611  };
1612 
1613  // clear data ranges
1614  trimEdges.clear();
1615  edgesToTrim.clear();
1616  verticesOnTrimEdges.clear();
1617  trimNewVertices.clear();
1618 
1619  auto cut_this_edge = [&](const EntityHandle e, const double length, auto &ray,
1620  auto &ray_point) {
1621  trimEdges.insert(e);
1622  edgesToTrim[e].dIst = 1;
1623  edgesToTrim[e].lEngth = length;
1624  edgesToTrim[e].unitRayDir.resize(3, false);
1625  edgesToTrim[e].rayPoint.resize(3, false);
1626  for (int ii = 0; ii != 3; ++ii)
1627  edgesToTrim[e].unitRayDir[ii] = ray(ii);
1628  for (int ii = 0; ii != 3; ++ii)
1629  edgesToTrim[e].rayPoint[ii] = ray_point(ii);
1630  };
1631 
1633  int num_nodes;
1634 
1635  auto get_tag_data = [&](auto th, auto conn) {
1637  CHKERR moab.tag_get_data(th, &conn, 1, &t(0));
1638  return t;
1639  };
1640 
1641  double max_edge_length = 0;
1642 
1643  if (!fRont.empty()) {
1644  // Calculate distances
1645  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
1646  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
1648 
1649  for (auto s : surface_skin) {
1650 
1651  auto project_on_surface = [&](auto &t) {
1653 
1654  EntityHandle facets_out;
1655  FTensor::Tensor1<double, 3> t_point_on_cutting_surface;
1656  CHKERR treeSurfPtr->closest_to_location(&t(0), rootSetSurf, &t_p(0),
1657  facets_out);
1658 
1659  FTensor::Tensor1<double,3> t_normal;
1660  CHKERR m_field.getInterface<Tools>()->getTriNormal(facets_out,
1661  &t_normal(0));
1662  t_normal(i) /= sqrt(t_normal(i) * t_normal(i));
1663  const double dot = t_normal(i) * (t_p(i) - t(i));
1664  t_p(i) = t(i) + dot * t_normal(i);
1665 
1666  return t_p;
1667  };
1668 
1669  const EntityHandle *conn;
1670  CHKERR moab.get_connectivity(s, conn, num_nodes, true);
1671 
1672  VectorDouble6 coords_front(6);
1673  CHKERR moab.get_coords(conn, num_nodes, &coords_front[0]);
1674 
1675  FTensor::Tensor1<double *, 3> t_f0(&coords_front[0], &coords_front[1],
1676  &coords_front[2]);
1677  FTensor::Tensor1<double *, 3> t_f1(&coords_front[3], &coords_front[4],
1678  &coords_front[5]);
1679 
1680  auto t_p_f0 = project_on_surface(t_f0);
1681  auto t_p_f1 = project_on_surface(t_f1);
1682 
1683  CHKERR moab.set_coords(&conn[0], 1, &t_p_f0(0));
1684  CHKERR moab.set_coords(&conn[1], 1, &t_p_f1(0));
1685  }
1686  }
1687 
1688  if (debug)
1689  CHKERR SaveData(moab)("surface_skin_to_trim.vtk", surface_skin);
1690 
1692  Tag th_dist_front_vec;
1693  CHKERR moab.tag_get_handle("DIST_FRONT_VECTOR", th_dist_front_vec);
1694 
1695  if (debug)
1696  CHKERR SaveData(moab)("edges_potentially_to_trim.vtk", cut_surface_edges);
1697 
1698  // iterate over edges on cut surface
1699  for (auto e : cut_surface_edges) {
1700 
1701  // Get edge connectivity and coords
1702  const EntityHandle *conn_edge;
1703  CHKERR moab.get_connectivity(e, conn_edge, num_nodes, true);
1704 
1705  auto t_dist0 = get_tag_data(th_dist_front_vec, conn_edge[0]);
1706  auto t_dist1 = get_tag_data(th_dist_front_vec, conn_edge[1]);
1707 
1708  double coords_edge[3 * num_nodes];
1709  CHKERR moab.get_coords(conn_edge, num_nodes, coords_edge);
1710 
1711  FTensor::Tensor1<double *, 3> t_e0(&coords_edge[0], &coords_edge[1],
1712  &coords_edge[2]);
1713  FTensor::Tensor1<double *, 3> t_e1(&coords_edge[3], &coords_edge[4],
1714  &coords_edge[5]);
1715 
1716  FTensor::Tensor1<double, 3> t_edge_delta;
1717  t_edge_delta(i) = t_e1(i) - t_e0(i);
1718  const double edge_length2 = t_edge_delta(i) * t_edge_delta(i);
1719  const double edge_length = sqrt(edge_length2);
1720  if (edge_length == 0)
1721  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Zero edge length");
1722 
1723  max_edge_length = std::max(max_edge_length, edge_length);
1724  const double s0 = t_dist0(i) * t_edge_delta(i) / edge_length;
1725  const double s1 = t_dist1(i) * t_edge_delta(i) / edge_length;
1726  const double dot = t_dist0(i) * t_dist1(i);
1727 
1728  // iterate entities surface front to find cross to trim
1729  for (auto s : surface_skin) {
1730 
1731  auto get_edge_coors = [&](const auto e) {
1732  const EntityHandle *conn;
1733  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1734  VectorDouble6 coords(6);
1735  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1736  return coords;
1737  };
1738 
1739  auto coords_front = get_edge_coors(s);
1740 
1741  FTensor::Tensor1<double *, 3> t_f0(&coords_front[0], &coords_front[1],
1742  &coords_front[2]);
1743  FTensor::Tensor1<double *, 3> t_f1(&coords_front[3], &coords_front[4],
1744  &coords_front[5]);
1745 
1746  // find point of minilam distance between front and cut surface edge
1747  double t_edge = -1, t_front = -1;
1748  auto res = Tools::minDistanceFromSegments(&t_e0(0), &t_e1(0), &t_f0(0),
1749  &t_f1(0), &t_edge, &t_front);
1750 
1751  if (res != Tools::NO_SOLUTION) {
1752  // check if edges crossing each other in the middle (it not imply that
1753  // have common point)
1754  const double overlap_tol = 1e-2;
1755  if (
1756 
1757  (t_edge > -std::numeric_limits<float>::epsilon() &&
1758  t_edge < 1 + std::numeric_limits<float>::epsilon())
1759 
1760  &&
1761 
1762  (t_front >= -overlap_tol && t_front <= 1 + overlap_tol)
1763 
1764  ) {
1765 
1766  t_front = std::max(0., std::min(t_front, 1.));
1767 
1768  FTensor::Tensor1<double, 3> t_front_delta;
1769  t_front_delta(i) = t_f1(i) - t_f0(i);
1770  FTensor::Tensor1<double, 3> t_edge_point, t_front_point;
1771  t_edge_point(i) = t_e0(i) + t_edge * t_edge_delta(i);
1772  t_front_point(i) = t_f0(i) + t_front * t_front_delta(i);
1773 
1775  t_ray(i) = t_front_point(i) - t_edge_point(i);
1776  const double dist = sqrt(t_ray(i) * t_ray(i));
1777 
1778  // that imply that edges have common point
1779  if ((dist < 0.1 * edge_length) ||
1780  (s0 >= 0 && s1 < -std::numeric_limits<double>::epsilon() &&
1781  dot <= 0)) {
1782 
1783  auto check_to_add_edge = [&](const EntityHandle e,
1784  const double dist) {
1785  auto eit = edgesToTrim.find(e);
1786  if (eit != edgesToTrim.end())
1787  if (eit->second.dIst < dist)
1788  return false;
1789  return true;
1790  };
1791 
1792  if (t_edge < 0.5) {
1793  t_ray(i) = t_edge * t_edge_delta(i);
1794  const double ray_length = sqrt(t_ray(i) * t_ray(i));
1795  if (check_to_add_edge(e, dist)) {
1796  add_vert(conn_edge[0], e, fabs(t_edge));
1797  add_vert(conn_edge[1], e, fabs(t_edge - 1));
1798  cut_this_edge(e, edge_length, t_ray, t_e0);
1799  }
1800  } else {
1801  FTensor::Tensor1<double, 3> t_edge_point;
1802  t_edge_point(i) = t_e0(i) + t_edge * t_edge_delta(i);
1803  t_ray(i) = t_edge_point(i) - t_e1(i);
1804  const double ray_length = sqrt(t_ray(i) * t_ray(i));
1805  if (check_to_add_edge(e, dist)) {
1806  add_vert(conn_edge[0], e, fabs(t_edge));
1807  add_vert(conn_edge[1], e, fabs(t_edge - 1));
1808  cut_this_edge(e, edge_length, t_ray, t_e1);
1809  }
1810  }
1811  }
1812  }
1813  }
1814  }
1815  }
1816 
1817  if (debug)
1818  CHKERR SaveData(moab)("edges_selected_to_trim.vtk", trimEdges);
1819 
1820  auto get_quality_change = [&](const Range &adj_tets, const EntityHandle &v,
1821  const VectorDouble3 &new_pos) {
1822  double q0 = 1;
1823  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0, th);
1824 
1825  double min_q = 1;
1826  for (auto t : adj_tets) {
1827  int num_nodes;
1828  const EntityHandle *conn;
1829  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1830  VectorDouble12 coords(12);
1831  if (th)
1832  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1833  else
1834  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1835 
1836  for (int n = 0; n != 4; ++n) {
1837  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1838  if (conn[n] == v) {
1839  noalias(n_coords) = new_pos;
1840  } else {
1841  auto m = verticesOnTrimEdges.find(conn[n]);
1842  if (m != verticesOnTrimEdges.end())
1843  noalias(n_coords) =
1844  m->second.rayPoint + m->second.dIst * m->second.unitRayDir;
1845  }
1846  }
1847 
1848  const double q = Tools::volumeLengthQuality(&coords[0]);
1849  if (!std::isnormal(q)) {
1850  min_q = -2;
1851  break;
1852  }
1853  min_q = std::min(min_q, q);
1854  }
1855  return min_q / q0;
1856  };
1857 
1858  Range checked_verts;
1859  auto add_trim_vert = [&](const EntityHandle v, const EntityHandle e) {
1860  if (!(trimNewVertices.find(v) != trimNewVertices.end())) {
1861  auto &r = edgesToTrim.at(e);
1862  VectorDouble3 ray_point = get_point_coords(v);
1863  VectorDouble3 new_pos = r.rayPoint + r.dIst * r.unitRayDir;
1864  VectorDouble3 unit_ray_dir = ray_point - new_pos;
1865  Range adj_tets;
1866  CHKERR moab.get_adjacencies(&v, 1, 3, false, adj_tets);
1867  adj_tets = intersect(adj_tets, cutNewVolumes);
1868  double q = get_quality_change(adj_tets, v, new_pos);
1870  VectorDouble3 ray_dir = new_pos - ray_point;
1871  double dist = norm_2(ray_dir);
1872  verticesOnTrimEdges[v].dIst = 1;
1873  verticesOnTrimEdges[v].lEngth = dist;
1874  verticesOnTrimEdges[v].unitRayDir = ray_dir;
1875  verticesOnTrimEdges[v].rayPoint = ray_point;
1876  trimNewVertices.insert(v);
1877  }
1878  checked_verts.insert(v);
1879  }
1880  };
1881 
1882  auto add_no_move_trim = [&](const EntityHandle v, const EntityHandle e) {
1883  if (!(trimNewVertices.find(v) != trimNewVertices.end())) {
1884  auto &m = edgesToTrim.at(e);
1885  verticesOnTrimEdges[v] = m;
1886  verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
1887  verticesOnTrimEdges[v].lEngth = 0;
1888  verticesOnTrimEdges[v].dIst = 0;
1889  trimNewVertices.insert(v);
1890  checked_verts.insert(v);
1891  }
1892  };
1893 
1894  // Iterate over all vertives close to surface front and check if those can
1895  // be moved
1896 
1897  // filer nodes which distance is in given tolerance
1898  auto hi = verts_map.get<0>().upper_bound(tol);
1899  verts_map.get<0>().erase(hi, verts_map.end());
1900 
1901  auto remove_verts = [&](Range nodes) {
1902  for (auto n : nodes) {
1903  auto r = verts_map.get<1>().equal_range(n);
1904  verts_map.get<1>().erase(r.first, r.second);
1905  }
1906  };
1907 
1908  auto trim_verts = [&](const Range verts, const bool move) {
1909  VertMapMultiIndex verts_map_tmp;
1910  for (auto p = corners.pair_begin(); p != corners.pair_end(); ++p) {
1911  auto lo = verts_map.get<1>().lower_bound(p->first);
1912  auto hi = verts_map.get<1>().upper_bound(p->second);
1913  verts_map_tmp.insert(lo, hi);
1914  }
1915  if (move) {
1916  for (auto &m : verts_map_tmp.get<0>())
1917  add_trim_vert(m.v, m.e);
1918  } else {
1919  for (auto &m : verts_map_tmp.get<0>())
1920  add_no_move_trim(m.v, m.e);
1921  }
1922  };
1923 
1924  auto trim_edges = [&](const Range ents, const bool move) {
1925  VertMapMultiIndex verts_map_tmp;
1926  for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
1927  auto lo = verts_map.get<2>().lower_bound(p->first);
1928  auto hi = verts_map.get<2>().upper_bound(p->second);
1929  verts_map_tmp.insert(lo, hi);
1930  verts_map.get<2>().erase(lo, hi);
1931  }
1932  if (move) {
1933  for (auto &m : verts_map_tmp.get<0>())
1934  add_trim_vert(m.v, m.e);
1935  } else {
1936  for (auto &m : verts_map_tmp.get<0>())
1937  add_no_move_trim(m.v, m.e);
1938  }
1939  };
1940 
1941  auto intersect_self = [&](Range &a, const Range b) { a = intersect(a, b); };
1942 
1943  map<std::string, Range> range_maps;
1944  CHKERR skin.find_skin(0, cutNewSurfaces, false, range_maps["surface_skin"]);
1945  intersect_self(range_maps["surface_skin"], trimEdges);
1946  range_maps["fixed_edges_on_surface_skin"] =
1947  intersect(range_maps["surface_skin"], fix_edges);
1948  CHKERR moab.get_adjacencies(range_maps["fixed_edges_verts"], 1, false,
1949  range_maps["fixed_edges_verts_edges"],
1950  moab::Interface::UNION);
1951  intersect_self(range_maps["fixed_edges_verts_edges"], trimEdges);
1952  CHKERR moab.get_connectivity(
1953  range_maps["fixed_edges_verts_edges"],
1954  range_maps["fixed_edges_verts_edges_verts_on_the_skin"], false);
1955  intersect_self(range_maps["fixed_edges_verts_edges_verts_on_the_skin"],
1956  tets_skin_verts);
1957 
1958  // do not move nodes at the corners
1959  trim_verts(corners, false);
1960  remove_verts(corners);
1961  trim_edges(range_maps["fixed_edges_on_surface_skin"], true);
1962  remove_verts(range_maps["fixed_edges_on_surface_skin_verts"]);
1963  trim_verts(range_maps["fixed_edges_verts_edges_verts_on_the_skin"], false);
1964  remove_verts(range_maps["fixed_edges_verts_edges_verts_on_the_skin"]);
1965  trim_edges(range_maps["surface_skin"], true);
1966  trim_verts(tets_skin_verts, false);
1967  remove_verts(tets_skin_verts);
1968 
1969  for (auto &m : verts_map.get<0>())
1970  add_trim_vert(m.v, m.e);
1971 
1972  for (auto v : subtract(cut_surface_verts, checked_verts)) {
1973 
1974  if (!(trimNewVertices.find(v) != trimNewVertices.end())) {
1975 
1976  auto get_tag_data = [&](auto th, auto conn) {
1978  CHKERR moab.tag_get_data(th, &conn, 1, &t(0));
1979  return t;
1980  };
1981 
1982  auto t_dist = get_tag_data(th_dist_front_vec, v);
1983  const double d = sqrt(t_dist(i) * t_dist(i));
1984  if (d < tol * max_edge_length) {
1985 
1986  Range adj;
1987  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj);
1988  adj = intersect(adj, cut_surface_edges);
1989  double min_length = max_edge_length;
1990  for (auto e : adj) {
1991 
1992  // Get edge connectivity and coords
1993  const EntityHandle *conn_edge;
1994  CHKERR moab.get_connectivity(e, conn_edge, num_nodes, true);
1995  double coords_edge[3 * num_nodes];
1996  CHKERR moab.get_coords(conn_edge, num_nodes, coords_edge);
1997  FTensor::Tensor1<double *, 3> t_e0(&coords_edge[0], &coords_edge[1],
1998  &coords_edge[2]);
1999  FTensor::Tensor1<double *, 3> t_e1(&coords_edge[3], &coords_edge[4],
2000  &coords_edge[5]);
2001  FTensor::Tensor1<double, 3> t_edge_delta;
2002  t_edge_delta(i) = t_e1(i) - t_e0(i);
2003 
2004  const double length = sqrt(t_edge_delta(i) * t_edge_delta(i));
2005  min_length = std::min(min_length, length);
2006  }
2007 
2008  if (d < (tol * tol * min_length)) {
2009  verticesOnTrimEdges[v].dIst = 0;
2010  verticesOnTrimEdges[v].lEngth = 0;
2011  verticesOnTrimEdges[v].unitRayDir.resize(3);
2012  verticesOnTrimEdges[v].unitRayDir.clear();
2013  verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
2014  trimNewVertices.insert(v);
2015  }
2016  }
2017  }
2018  }
2019 
2020  for (auto m : verticesOnTrimEdges) {
2021  EntityHandle v = m.first;
2022  Range adj;
2023  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj);
2024  for (auto e : adj) {
2025  edgesToTrim.erase(e);
2026  trimEdges.erase(e);
2027  }
2028  }
2029 
2030  if (debug)
2031  if (!trimNewVertices.empty())
2032  CHKERR SaveData(moab)("trim_close_vertices.vtk", trimNewVertices);
2033 
2034  if (debug)
2035  if (!trimEdges.empty())
2036  CHKERR SaveData(moab)("trim_edges.vtk", trimEdges);
2037 
2039 }
map< EntityHandle, TreeData > edgesToTrim
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:32
moab::Interface & get_moab()
Definition: Core.hpp:266
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:342
map< EntityHandle, TreeData > verticesOnTrimEdges
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:88
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Types.hpp:90
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
double tol
static const bool debug
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:43
#define CHKERR
Inline error check.
Definition: definitions.h:596
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
MoFEMErrorCode createFrontLevelSets(Range vol, Tag th=nullptr, int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to surface front.

◆ getCutEdges()

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

Definition at line 470 of file CutMeshInterface.hpp.

470 { return cutEdges; }

◆ getCutFrontVolumes()

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

Definition at line 487 of file CutMeshInterface.hpp.

◆ getCutSurfaceVolumes()

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

Definition at line 486 of file CutMeshInterface.hpp.

◆ getFront()

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

Definition at line 468 of file CutMeshInterface.hpp.

468 { return fRont; }

◆ getMergedSurfaces()

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

Definition at line 484 of file CutMeshInterface.hpp.

◆ getMergedVolumes()

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

Definition at line 483 of file CutMeshInterface.hpp.

483 { return mergedVolumes; }

◆ getNewCutSurfaces()

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

Definition at line 472 of file CutMeshInterface.hpp.

◆ getNewCutVertices()

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

Definition at line 473 of file CutMeshInterface.hpp.

◆ getNewCutVolumes()

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

Definition at line 471 of file CutMeshInterface.hpp.

471 { return cutNewVolumes; }

◆ getNewTrimSurfaces()

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

Definition at line 480 of file CutMeshInterface.hpp.

◆ getNewTrimVertices()

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

Definition at line 481 of file CutMeshInterface.hpp.

◆ getNewTrimVolumes()

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

Definition at line 479 of file CutMeshInterface.hpp.

◆ getOptions()

MoFEMErrorCode MoFEM::CutMeshInterface::getOptions ( )

Get options from command line.

Returns
error code

Definition at line 51 of file CutMeshInterface.hpp.

51  {
53  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "cut_", "MOFEM Cut mesh options",
54  "none");
55 
56  CHKERR PetscOptionsInt("-linesearch_steps",
57  "number of bisection steps which line search do to "
58  "find optimal merged nodes position",
59  "", lineSearchSteps, &lineSearchSteps, PETSC_NULL);
60 
61  CHKERR PetscOptionsInt("-max_merging_cycles",
62  "number of maximal merging cycles", "",
64 
65  CHKERR PetscOptionsScalar("-project_entities_quality_trashold",
66  "project entities quality trashold", "",
68  &projectEntitiesQualityTrashold, PETSC_NULL);
69 
70  ierr = PetscOptionsEnd();
71  CHKERRG(ierr);
73  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ getSubInterfaceOptions()

MoFEMErrorCode MoFEM::CutMeshInterface::getSubInterfaceOptions ( )

Definition at line 45 of file CutMeshInterface.hpp.

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

◆ getSurface()

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

Definition at line 467 of file CutMeshInterface.hpp.

467 { return sUrface; }

◆ getTreeSurfPtr()

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

Definition at line 495 of file CutMeshInterface.hpp.

495  {
496  return treeSurfPtr;
497  }
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr

◆ getTrimEdges()

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

Definition at line 478 of file CutMeshInterface.hpp.

478 { return trimEdges; }

◆ getVolume()

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

Definition at line 466 of file CutMeshInterface.hpp.

466 { 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 445 of file CutMeshInterface.cpp.

445  {
446  CoreInterface &m_field = cOre;
447  moab::Interface &moab = m_field.get_moab();
449  Skinner skin(&moab);
450  Range tets_skin;
451  CHKERR skin.find_skin(0, vOlume, false, tets_skin);
452  Range tets_skin_edges;
453  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
454  moab::Interface::UNION);
455  Range surface_skin;
456  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
457  fRont = subtract(surface_skin, tets_skin_edges);
458  if (debug)
459  CHKERR SaveData(moab)("front.vtk", fRont);
461 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ mergeBadEdges() [1/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 2298 of file CutMeshInterface.cpp.

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

◆ mergeBadEdges() [2/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 3044 of file CutMeshInterface.cpp.

3048  {
3049  CoreInterface &m_field = cOre;
3051  Range tets_level;
3052  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3053  trim_bit, BitRefLevel().set(), MBTET, tets_level);
3054 
3055  Range edges_to_merge;
3056  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
3057  trim_bit, cut_bit | trim_bit, edges_to_merge);
3058 
3059  // get all entities not in database
3060  Range all_ents_not_in_database_before;
3061  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
3062  all_ents_not_in_database_before);
3063 
3064  edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
3065  if (debug)
3066  CHKERR SaveData(m_field.get_moab())("edges_to_merge.vtk", edges_to_merge);
3067 
3068  Range out_new_tets, out_new_surf;
3069  CHKERR mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
3070  corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
3071  th, update_meshsets, &bit, debug);
3072 
3073  // get all entities not in database after merge
3074  Range all_ents_not_in_database_after;
3075  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
3076  all_ents_not_in_database_after);
3077 
3078  // delete hanging entities
3079  all_ents_not_in_database_after =
3080  subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
3081  Range meshsets;
3082  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
3083  false);
3084  for (auto m : meshsets)
3085  CHKERR m_field.get_moab().remove_entities(m,
3086  all_ents_not_in_database_after);
3087 
3088  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
3089 
3090  mergedVolumes.swap(out_new_tets);
3091  mergedSurfaces.swap(out_new_surf);
3093 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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 getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ mergeSurface()

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

merge surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 125 of file CutMeshInterface.cpp.

125  {
127  sUrface.merge(surface);
129 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508

◆ mergeVolumes()

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

merge volume entities

Parameters
volumeentities which going to be added
Returns
error code

Definition at line 131 of file CutMeshInterface.cpp.

131  {
133  vOlume.merge(volume);
135 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508

◆ moveMidNodesOnCutEdges()

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

projecting of mid edge nodes on new mesh on surface

Returns
error code

Definition at line 1499 of file CutMeshInterface.cpp.

1499  {
1501 
1502  CoreInterface &m_field = cOre;
1503  moab::Interface &moab = m_field.get_moab();
1505 
1506  // Range out_side_vertices;
1507  for (auto m : verticesOnCutEdges) {
1508  double dist = m.second.dIst;
1509  VectorDouble3 new_coors = m.second.rayPoint + dist * m.second.unitRayDir;
1510  if (th)
1511  CHKERR moab.tag_set_data(th, &m.first, 1, &new_coors[0]);
1512  else
1513  CHKERR moab.set_coords(&m.first, 1, &new_coors[0]);
1514  }
1515 
1517 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ moveMidNodesOnTrimmedEdges()

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

move trimmed edges mid nodes

Returns
error code

Definition at line 1519 of file CutMeshInterface.cpp.

1519  {
1520  CoreInterface &m_field = cOre;
1521  moab::Interface &moab = m_field.get_moab();
1523  for (auto &v : verticesOnTrimEdges) {
1524  double dist = v.second.dIst;
1525  VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1526  if (th)
1527  CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1528  else
1529  CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1530  }
1532 }
moab::Interface & get_moab()
Definition: Core.hpp:266
map< EntityHandle, TreeData > verticesOnTrimEdges
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ projectZeroDistanceEnts() [1/2]

MoFEMErrorCode MoFEM::CutMeshInterface::projectZeroDistanceEnts ( Range *  fixed_edges,
Range *  corner_nodes,
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 1034 of file CutMeshInterface.cpp.

1038  {
1039  CoreInterface &m_field = cOre;
1040  moab::Interface &moab = m_field.get_moab();
1042 
1043  // Get entities on body skin
1044  Skinner skin(&moab);
1045  Range tets_skin;
1046  CHKERR skin.find_skin(0, vOlume, false, tets_skin);
1047  Range tets_skin_edges;
1048  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1049  moab::Interface::UNION);
1050  Range tets_skin_verts;
1051  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1052 
1053  // Get entities in volume
1054  Range vol_faces, vol_edges, vol_nodes;
1055  CHKERR moab.get_adjacencies(vOlume, 2, false, vol_faces,
1056  moab::Interface::UNION);
1057  CHKERR moab.get_adjacencies(vOlume, 1, false, vol_edges,
1058  moab::Interface::UNION);
1059  CHKERR moab.get_connectivity(vOlume, vol_nodes, true);
1060 
1061  // Get nodes on cut edges
1062  Range cut_edge_verts;
1063  CHKERR moab.get_connectivity(cutEdges, cut_edge_verts, true);
1064 
1065  // Get faces and edges
1066  Range cut_edges_faces;
1067  CHKERR moab.get_adjacencies(cut_edge_verts, 2, true, cut_edges_faces,
1068  moab::Interface::UNION);
1069  cut_edges_faces = intersect(cut_edges_faces, vol_faces);
1070  Range cut_edges_faces_verts;
1071  CHKERR moab.get_connectivity(cut_edges_faces, cut_edges_faces_verts, true);
1072  cut_edges_faces_verts = subtract(cut_edges_faces_verts, cut_edge_verts);
1073  Range to_remove_cut_edges_faces;
1074  CHKERR moab.get_adjacencies(cut_edges_faces_verts, 2, true,
1075  to_remove_cut_edges_faces,
1076  moab::Interface::UNION);
1077  // Those are faces which have vertices adjacent to cut edges vertices without
1078  // hanging nodes (nodes not adjacent to cutting edge)
1079  cut_edges_faces = subtract(cut_edges_faces, to_remove_cut_edges_faces);
1080  if (debug)
1081  CHKERR SaveData(moab)("cut_edges_faces.vtk", cut_edges_faces);
1082  cut_edges_faces.merge(cutEdges);
1083 
1084  Range fixed_edges_nodes;
1085  if (fixed_edges)
1086  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_nodes, true);
1087 
1088  Tag th_dist_normal;
1089  CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_VECTOR", th_dist_normal);
1090 
1091  // Create map of closes points to the surface
1092  enum TYPE { FREE = 0, SKIN, FIXED, CORNER, NOT_MOVE };
1093  map<EntityHandle, std::pair<std::pair<TYPE, EntityHandle>, TreeData>>
1094  min_dist_map;
1095  double ave_cut_edge_length = 0;
1096  for (auto e : cutEdges) {
1097 
1098  auto eit = edgesToCut.find(e);
1099  if (eit != edgesToCut.end()) {
1100 
1101  TYPE edge_type = FREE;
1102  if (tets_skin_edges.find(e) != tets_skin_edges.end())
1103  edge_type = SKIN;
1104  if (fixed_edges)
1105  if (fixed_edges->find(e) != fixed_edges->end())
1106  edge_type = FIXED;
1107 
1108  int num_nodes;
1109  const EntityHandle *conn;
1110  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1111  VectorDouble6 pos(6);
1112  CHKERR moab.get_coords(conn, num_nodes, &pos[0]);
1113  VectorDouble3 p[2];
1114  p[0] = getVectorAdaptor(&pos[0], 3);
1115  p[1] = getVectorAdaptor(&pos[3], 3);
1116  ave_cut_edge_length += norm_2(p[0] - p[1]);
1117 
1118  auto &d = eit->second;
1119  VectorDouble3 new_pos = d.rayPoint + d.dIst * d.unitRayDir;
1120  for (int nn = 0; nn != 2; ++nn) {
1121 
1122  VectorDouble3 ray = new_pos - p[nn];
1123  const double dist = norm_2(ray);
1124  const double length = dist;
1125 
1126  bool add_node = true;
1127  auto vit = min_dist_map.find(conn[nn]);
1128  if (vit != min_dist_map.end()) {
1129  if (vit->second.second.dIst < dist)
1130  add_node = false;
1131  }
1132 
1133  if (add_node) {
1134  min_dist_map[conn[nn]].first.first = edge_type;
1135  min_dist_map[conn[nn]].first.second = e;
1136  auto &data = min_dist_map[conn[nn]].second;
1137  data.lEngth = length;
1138  data.rayPoint = p[nn];
1139  data.dIst = dist;
1140  if (dist > 0)
1141  data.unitRayDir = ray / dist;
1142  else {
1143  data.unitRayDir.resize(3);
1144  data.unitRayDir.clear();
1145  }
1146  }
1147  }
1148 
1149  } else
1150  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Edge not found");
1151  }
1152 
1153  ave_cut_edge_length /= cutEdges.size();
1154 
1155  auto get_min_quality =
1156  [&](const Range &adj_tets,
1157  const map<EntityHandle, TreeData> &vertices_on_cut_edges) {
1158  double min_q = 1;
1159  for (auto t : adj_tets) {
1160  int num_nodes;
1161  const EntityHandle *conn;
1162  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1163  VectorDouble12 coords(12);
1164  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1165  for (auto n : {0, 1, 2, 3}) {
1166  auto mit = vertices_on_cut_edges.find(conn[n]);
1167  if (mit != vertices_on_cut_edges.end()) {
1168  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1169  const double dist = mit->second.dIst;
1170  noalias(n_coords) =
1171  mit->second.rayPoint + dist * mit->second.unitRayDir;
1172  }
1173  }
1174  min_q = std::min(min_q, Tools::volumeLengthQuality(&coords[0]));
1175  }
1176  return min_q;
1177  };
1178 
1179  auto get_quality_change =
1180  [&](const Range &adj_tets,
1181  map<EntityHandle, TreeData> vertices_on_cut_edges) {
1182  double q0 = get_min_quality(adj_tets, verticesOnCutEdges);
1183  vertices_on_cut_edges.insert(verticesOnCutEdges.begin(),
1184  verticesOnCutEdges.end());
1185  double q = get_min_quality(adj_tets, vertices_on_cut_edges);
1186  return q / q0;
1187  };
1188 
1189  auto get_conn = [&moab](const EntityHandle &e, const EntityHandle *&conn,
1190  int &num_nodes) {
1192  EntityType type = moab.type_from_handle(e);
1193  if (type == MBVERTEX) {
1194  conn = &e;
1195  num_nodes = 1;
1196  } else {
1197  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1198  }
1200  };
1201 
1202  auto project_node = [&](const EntityHandle v,
1203  map<EntityHandle, TreeData> &vertices_on_cut_edges) {
1205 
1206  vertices_on_cut_edges[v].dIst = min_dist_map[v].second.dIst;
1207  vertices_on_cut_edges[v].lEngth = min_dist_map[v].second.lEngth;
1208  vertices_on_cut_edges[v].unitRayDir = min_dist_map[v].second.unitRayDir;
1209  vertices_on_cut_edges[v].rayPoint = min_dist_map[v].second.rayPoint;
1210 
1212  };
1213 
1214  auto remove_surface_tets = [&](Range &zero_dist_tets,
1215  Range zero_distance_ents,
1216  Range zero_distance_verts) {
1218  Range zero_dist_all_verts;
1219  CHKERR moab.get_connectivity(zero_distance_ents, zero_dist_all_verts, true);
1220  zero_dist_all_verts.merge(zero_distance_verts);
1221  CHKERR moab.get_adjacencies(zero_dist_all_verts, 3, false, zero_dist_tets,
1222  moab::Interface::UNION);
1223  zero_dist_tets = intersect(zero_dist_tets, vOlume);
1224  Range zero_tets_verts;
1225  CHKERR moab.get_connectivity(zero_dist_tets, zero_tets_verts, true);
1226  zero_tets_verts = subtract(zero_tets_verts, zero_dist_all_verts);
1227  Range free_tets;
1228  CHKERR moab.get_adjacencies(zero_tets_verts, 3, false, free_tets,
1229  moab::Interface::UNION);
1230  zero_dist_tets = subtract(zero_dist_tets, free_tets);
1231 
1233  };
1234 
1235  for (int d = 2; d >= 0; --d) {
1236 
1237  Range ents;
1238  if (d > 0)
1239  ents = cut_edges_faces.subset_by_dimension(d);
1240  else
1241  ents = cut_edge_verts;
1242 
1243  // make list of entities
1244  multimap<double, EntityHandle> ents_to_check;
1245  for (auto f : ents) {
1246  int num_nodes;
1247  const EntityHandle *conn;
1248  CHKERR get_conn(f, conn, num_nodes);
1249  VectorDouble3 dist(3);
1250 
1251  for (int n = 0; n != num_nodes; ++n) {
1252  auto &d = min_dist_map[conn[n]];
1253  dist[n] = d.second.lEngth;
1254  }
1255 
1256  double max_dist = 0;
1257  for (int n = 0; n != num_nodes; ++n)
1258  max_dist = std::max(max_dist, fabs(dist[n]));
1259  if (max_dist < close_tol * ave_cut_edge_length)
1260  ents_to_check.insert(std::pair<double, EntityHandle>(max_dist, f));
1261  }
1262 
1263  if (debug) {
1264  Range ents;
1265  for (auto m : ents_to_check)
1266  ents.insert(m.second);
1267  CHKERR SaveData(moab)("ents_to_check_to_project.vtk", ents);
1268  }
1269 
1270  double ray_point[3], unit_ray_dir[3];
1271  VectorAdaptor vec_unit_ray_dir(
1272  3, ublas::shallow_array_adaptor<double>(3, unit_ray_dir));
1273  VectorAdaptor vec_ray_point(
1274  3, ublas::shallow_array_adaptor<double>(3, ray_point));
1275 
1276  for (auto m : ents_to_check) {
1277 
1278  EntityHandle f = m.second;
1279 
1280  int num_nodes;
1281  const EntityHandle *conn;
1282  CHKERR get_conn(f, conn, num_nodes);
1283  VectorDouble9 coords(9);
1284  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1285 
1286  Range adj_tets;
1287  CHKERR moab.get_adjacencies(conn, num_nodes, 3, false, adj_tets,
1288  moab::Interface::UNION);
1289  adj_tets = intersect(adj_tets, vOlume);
1290 
1291  map<EntityHandle, TreeData> vertices_on_cut_edges;
1292  for (int n = 0; n != num_nodes; ++n)
1293  CHKERR project_node(conn[n], vertices_on_cut_edges);
1294 
1295  const double q = get_quality_change(adj_tets, vertices_on_cut_edges);
1296  if (std::isnormal(q) && q > projectEntitiesQualityTrashold) {
1297  EntityHandle type = moab.type_from_handle(f);
1298 
1299  Range zero_dist_tets;
1300  if (type == MBVERTEX) {
1301  Range zero_distance_verts_test = zeroDistanceVerts;
1302  zero_distance_verts_test.insert(f);
1303  CHKERR remove_surface_tets(zero_dist_tets, zeroDistanceEnts,
1304  zero_distance_verts_test);
1305  } else {
1306  Range zero_distance_ents_test = zeroDistanceEnts;
1307  zero_distance_ents_test.insert(f);
1308  CHKERR remove_surface_tets(zero_dist_tets, zero_distance_ents_test,
1310  }
1311 
1312  if (zero_dist_tets.empty()) {
1313  verticesOnCutEdges.insert(vertices_on_cut_edges.begin(),
1314  vertices_on_cut_edges.end());
1315  if (type == MBVERTEX) {
1316  zeroDistanceVerts.insert(f);
1317  } else {
1318  zeroDistanceEnts.insert(f);
1319  }
1320  }
1321  }
1322  }
1323  }
1324 
1325  for (auto &v : verticesOnCutEdges) {
1326 
1327  TYPE node_type;
1328 
1329  if (corner_nodes->find(v.first) != corner_nodes->end())
1330  node_type = CORNER;
1331  else if (fixed_edges_nodes.find(v.first) != fixed_edges_nodes.end())
1332  node_type = FIXED;
1333  else if (tets_skin_verts.find(v.first) != tets_skin_verts.end())
1334  node_type = SKIN;
1335  else
1336  node_type = FREE;
1337 
1338  if (node_type > min_dist_map[v.first].first.first)
1339  v.second.unitRayDir.clear();
1340  }
1341 
1342  for (auto f : unite(zeroDistanceEnts, zeroDistanceVerts)) {
1343  int num_nodes;
1344  const EntityHandle *conn;
1345  CHKERR get_conn(f, conn, num_nodes);
1346  Range adj_edges;
1347  CHKERR moab.get_adjacencies(conn, num_nodes, 1, false, adj_edges,
1348  moab::Interface::UNION);
1349  for (auto e : adj_edges) {
1350  cutEdges.erase(e);
1351  edgesToCut.erase(e);
1352  }
1353  }
1354 
1355  if (debug)
1356  SaveData(m_field.get_moab())("projected_zero_distance_ents.vtk",
1358 
1360 }
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:32
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:88
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Types.hpp:90
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
static const bool debug
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:43
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:101
#define CHKERR
Inline error check.
Definition: definitions.h:596
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
VectorBoundedArray< double, 9 > VectorDouble9
Definition: Types.hpp:89
map< EntityHandle, TreeData > edgesToCut

◆ projectZeroDistanceEnts() [2/2]

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

Definition at line 474 of file CutMeshInterface.hpp.

474  {
475  return zeroDistanceEnts;
476  }

◆ query_interface()

MoFEMErrorCode MoFEM::CutMeshInterface::query_interface ( const MOFEMuuid uuid,
UnknownInterface **  iface 
) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 21 of file CutMeshInterface.cpp.

22  {
24  *iface = NULL;
25  if (uuid == IDD_MOFEMCutMesh) {
26  *iface = const_cast<CutMeshInterface *>(this);
28  }
29  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
31 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static const MOFEMuuid IDD_MOFEMCutMesh

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

718  {
719  CoreInterface &m_field = cOre;
720  moab::Interface &moab = m_field.get_moab();
721  MeshRefinement *refiner;
722  BitRefManager *bit_ref_manager;
724  CHKERR m_field.getInterface(refiner);
725  CHKERR m_field.getInterface(bit_ref_manager);
726 
727  auto add_bit = [&](const int bit) {
729  CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
730  verb);
731  Range adj_ents;
732  for (auto d : {2, 1, 0})
733  CHKERR moab.get_adjacencies(vOlume, d, true, adj_ents,
734  moab::Interface::UNION);
735  CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
736  verb);
738  };
739  CHKERR add_bit(init_bit_level);
740 
741  auto update_range = [&](Range *r_ptr) {
743  if (r_ptr) {
744  Range childs;
745  CHKERR bit_ref_manager->updateRange(*r_ptr, childs);
746  r_ptr->merge(childs);
747  }
749  };
750 
751  auto refine = [&](const BitRefLevel bit, const Range tets) {
753  Range verts;
754  CHKERR moab.get_connectivity(tets, verts, true);
755  Range ref_edges;
756  CHKERR moab.get_adjacencies(verts, 1, true, ref_edges,
757  moab::Interface::UNION);
758 
759  CHKERR refiner->add_vertices_in_the_middle_of_edges(ref_edges, bit);
760  CHKERR refiner->refine_TET(vOlume, bit, false, verb);
761 
762  CHKERR update_range(fixed_edges);
763  CHKERR update_range(&vOlume);
764  CHKERR m_field.getInterface<MeshsetsManager>()
765  ->updateAllMeshsetsByEntitiesChildren(bit);
766 
767  Range bit_tets;
768  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
769  bit, BitRefLevel().set(), MBTET, bit_tets);
770  vOlume = intersect(vOlume, bit_tets);
771 
772  Range bit_edges;
773  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
774  bit, BitRefLevel().set(), MBEDGE, bit_edges);
775  if (fixed_edges)
776  *fixed_edges = intersect(*fixed_edges, bit_edges);
777 
779  };
780 
781  for (int ll = init_bit_level; ll != init_bit_level + surf_levels; ++ll) {
783  CHKERR refine(BitRefLevel().set(ll + 1),
785  }
786 
787  for (int ll = init_bit_level + surf_levels;
788  ll != init_bit_level + surf_levels + front_levels; ++ll) {
790  CHKERR refine(BitRefLevel().set(ll + 1), cutFrontVolumes);
791  }
792 
794 
795  if (debug)
796  CHKERR SaveData(m_field.get_moab())("refinedVolume.vtk", vOlume);
797 
799 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
MoFEMErrorCode createLevelSets(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())
Create a level sets, i.e. distances from surface and surface front.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

2232  {
2233  CoreInterface &m_field = cOre;
2234  moab::Interface &moab = m_field.get_moab();
2235  PrismInterface *interface;
2237  CHKERR m_field.getInterface(interface);
2238  // Remove tris on surface front
2239  {
2240  Range front_tris;
2241  EntityHandle meshset;
2242  CHKERR moab.create_meshset(MESHSET_SET, meshset);
2243  CHKERR moab.add_entities(meshset, ents);
2244  CHKERR interface->findIfTringleHasThreeNodesOnInternalSurfaceSkin(
2245  meshset, split_bit, true, front_tris);
2246  CHKERR moab.delete_entities(&meshset, 1);
2247  ents = subtract(ents, front_tris);
2248  }
2249  // Remove entities on skin
2250  Range tets;
2251  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2252  split_bit, BitRefLevel().set(), MBTET, tets);
2253  // Remove entities on skin
2254  Skinner skin(&moab);
2255  Range tets_skin;
2256  CHKERR skin.find_skin(0, tets, false, tets_skin);
2257  ents = subtract(ents, tets_skin);
2258 
2260 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ saveCutEdges()

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

Definition at line 3128 of file CutMeshInterface.cpp.

3128  {
3129  CoreInterface &m_field = cOre;
3130  moab::Interface &moab = m_field.get_moab();
3132  CHKERR SaveData(moab)(prefix + "out_vol.vtk", vOlume);
3133  CHKERR SaveData(moab)(prefix + "out_surface.vtk", sUrface);
3134  CHKERR SaveData(moab)(prefix + "out_cut_edges.vtk", cutEdges);
3135  CHKERR SaveData(moab)(prefix + "out_cut_new_volumes.vtk", cutNewVolumes);
3136  CHKERR SaveData(moab)(prefix + "out_cut_new_surfaces.vtk", cutNewSurfaces);
3137  CHKERR SaveData(moab)(prefix + "out_cut_zero_distance_ents.vtk",
3140 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ saveTrimEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveTrimEdges ( )

Definition at line 3142 of file CutMeshInterface.cpp.

3142  {
3143  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
3145  CHKERR SaveData(moab)("out_trim_surface.vtk", sUrface);
3146  CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
3147  CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
3148  CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
3150 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

3112  {
3113  CoreInterface &m_field = cOre;
3114  moab::Interface &moab = m_field.get_moab();
3116  Range nodes;
3117  if (bit.none())
3118  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3119  else
3120  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3121  bit, mask, MBVERTEX, nodes);
3122  std::vector<double> coords(3 * nodes.size());
3123  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
3124  CHKERR moab.set_coords(nodes, &coords[0]);
3126 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ setFront()

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

Definition at line 47 of file CutMeshInterface.cpp.

47  {
49  fRont = front;
51 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508

◆ setSurface()

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

set surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 53 of file CutMeshInterface.cpp.

53  {
55  sUrface = surface;
57 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508

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

3095  {
3096  CoreInterface &m_field = cOre;
3097  moab::Interface &moab = m_field.get_moab();
3099  Range nodes;
3100  if (bit.none())
3101  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3102  else
3103  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3104  bit, BitRefLevel().set(), MBVERTEX, nodes);
3105  std::vector<double> coords(3 * nodes.size());
3106  CHKERR moab.get_coords(nodes, &coords[0]);
3107  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
3109 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ setTrimFixedEdges()

void MoFEM::CutMeshInterface::setTrimFixedEdges ( const bool  b)

Definition at line 489 of file CutMeshInterface.hpp.

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

119  {
121  vOlume = volume;
123 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508

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

139  {
140  CoreInterface &m_field = cOre;
141  moab::Interface &moab = m_field.get_moab();
143 
144  // Get cutting surface skin
145  Skinner skin(&moab);
146  Range surface_skin;
147  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
148 
149  CHKERR snapSurfaceToEdges(surface_skin, fixed_edges, rel_tol, abs_tol, th,
150  debug);
151 
153 }
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)
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

159  {
160  CoreInterface &m_field = cOre;
161  moab::Interface &moab = m_field.get_moab();
164 
165  map<EntityHandle, double> map_verts_length;
166 
167  for (auto f : surface_edges) {
168  int num_nodes;
169  const EntityHandle *conn_skin;
170  CHKERR moab.get_connectivity(f, conn_skin, num_nodes, true);
171  VectorDouble6 coords_skin(6);
172  if (th)
173  CHKERR moab.tag_get_data(th, conn_skin, num_nodes, &coords_skin[0]);
174  else
175  CHKERR moab.get_coords(conn_skin, num_nodes, &coords_skin[0]);
177  &coords_skin[0], &coords_skin[1], &coords_skin[2]);
179  &coords_skin[3], &coords_skin[4], &coords_skin[5]);
180  FTensor::Tensor1<double, 3> t_skin_delta;
181  t_skin_delta(i) = t_n1(i) - t_n0(i);
182  const double skin_edge_length = sqrt(t_skin_delta(i) * t_skin_delta(i));
183  for (int nn = 0; nn != 2; ++nn) {
184  auto it = map_verts_length.find(conn_skin[nn]);
185  if (it != map_verts_length.end())
186  it->second = std::min(it->second, skin_edge_length);
187  else
188  map_verts_length[conn_skin[nn]] = skin_edge_length;
189  }
190  }
191 
192  for (auto m : map_verts_length) {
193 
195  if (th)
196  CHKERR moab.tag_get_data(th, &m.first, 1, &t_n(0));
197  else
198  CHKERR moab.get_coords(&m.first, 1, &t_n(0));
199 
200  double min_dist = rel_tol * m.second;
201  FTensor::Tensor1<double, 3> t_min_coords;
202  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
203  &t_n(0), 1, fixed_edges, &min_dist, &t_min_coords(0));
204 
205  if (min_dist < rel_tol * m.second || min_dist < abs_tol) {
206  if(debug)
207  cerr << "Snap " << min_dist << endl;
208  if (th)
209  CHKERR moab.tag_set_data(th, &m.first, 1, &t_min_coords(0));
210  else
211  CHKERR moab.set_coords(&m.first, 1, &t_min_coords(0));
212  }
213  }
214 
216 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:88
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

2264  {
2265  CoreInterface &m_field = cOre;
2266  moab::Interface &moab = m_field.get_moab();
2267  PrismInterface *interface;
2269  CHKERR m_field.getInterface(interface);
2270  EntityHandle meshset_volume;
2271  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
2272  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2273  split_bit, BitRefLevel().set(), MBTET, meshset_volume);
2274  EntityHandle meshset_surface;
2275  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
2276  CHKERR moab.add_entities(meshset_surface, ents);
2277  CHKERR interface->getSides(meshset_surface, split_bit, true);
2278  CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
2279  true);
2280  CHKERR moab.delete_entities(&meshset_volume, 1);
2281  CHKERR moab.delete_entities(&meshset_surface, 1);
2282  if (th) {
2283  Range prisms;
2284  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2285  bit, BitRefLevel().set(), MBPRISM, prisms);
2286  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
2287  int num_nodes;
2288  const EntityHandle *conn;
2289  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
2290  MatrixDouble data(3, 3);
2291  CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
2292  CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
2293  }
2294  }
2296 }
moab::Interface & get_moab()
Definition: Core.hpp:266
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

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

2042  {
2043  CoreInterface &m_field = cOre;
2044  moab::Interface &moab = m_field.get_moab();
2045  MeshRefinement *refiner;
2046  const RefEntity_multiIndex *ref_ents_ptr;
2048 
2049  CHKERR m_field.getInterface(refiner);
2050  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
2051  CHKERR refiner->add_vertices_in_the_middle_of_edges(trimEdges, bit);
2052  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
2053  debug ? &trimEdges : NULL);
2054 
2055  trimNewVolumes.clear();
2056  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2057  bit, BitRefLevel().set(), MBTET, trimNewVolumes);
2058 
2059  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
2060  mit != edgesToTrim.end(); mit++) {
2061  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
2062  boost::make_tuple(MBVERTEX, mit->first));
2063  if (vit ==
2064  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end())
2065  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2066  "No vertex on trim edges, that make no sense");
2067 
2068  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
2069  if ((ref_ent->getBitRefLevel() & bit).any()) {
2070  EntityHandle vert = ref_ent->getRefEnt();
2071  trimNewVertices.insert(vert);
2072  verticesOnTrimEdges[vert] = mit->second;
2073  }
2074  }
2075 
2076  // Get faces which are trimmed
2077  trimNewSurfaces.clear();
2078  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2079  bit, bit, MBTRI, trimNewSurfaces);
2080 
2081  Range trim_new_surfaces_nodes;
2082  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
2083  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
2084  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
2085  Range faces_not_on_surface;
2086  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
2087  faces_not_on_surface, moab::Interface::UNION);
2088  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
2089 
2090  // Get surfaces which are not trimmed and add them to surface
2091  Range all_surfaces_on_bit_level;
2092  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2093  bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
2094  all_surfaces_on_bit_level =
2095  intersect(all_surfaces_on_bit_level, cutNewSurfaces);
2096 
2097  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
2098 
2100 }
map< EntityHandle, TreeData > edgesToTrim
moab::Interface & get_moab()
Definition: Core.hpp:266
map< EntityHandle, TreeData > verticesOnTrimEdges
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ trimOnly()

MoFEMErrorCode MoFEM::CutMeshInterface::trimOnly ( const BitRefLevel  trim_bit,
Tag  th,
const double  tol_cut_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);
275  CHKERR trimEdgesInTheMiddle(trim_bit, debug);
276  if (fixed_edges)
277  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
278  *fixed_edges);
279 
280  if (corner_nodes)
281  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
282  *corner_nodes);
283 
284  if (update_meshsets)
285  CHKERR m_field.getInterface<MeshsetsManager>()
286  ->updateAllMeshsetsByEntitiesChildren(trim_bit);
287 
288  // move nodes
290 
291  // remove faces
292  CHKERR trimSurface(fixed_edges, corner_nodes, debug);
293 
294  if (debug) {
296  Range bit2_edges;
297  CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
298  trim_bit, BitRefLevel().set(), MBEDGE, bit2_edges);
299  CHKERR SaveData(moab)("trim_fixed_edges.vtk",
300  intersect(*fixed_edges, bit2_edges));
301  }
302 
304 }
MoFEMErrorCode trimSurface(Range *fixed_edge, Range *corner_nodes, const bool debug=false)
Trim surface from faces beyond front.
MoFEMErrorCode trimEdgesInTheMiddle(const BitRefLevel bit, const bool debug=false)
trim edges
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, const bool debug=false)
Find edges to trimEdges.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
MoFEMErrorCode saveTrimEdges()
MoFEMErrorCode moveMidNodesOnTrimmedEdges(Tag th=NULL)
move trimmed edges mid nodes

◆ trimSurface()

MoFEMErrorCode MoFEM::CutMeshInterface::trimSurface ( Range *  fixed_edge,
Range *  corner_nodes,
const bool  debug = false 
)

Trim surface from faces beyond front.

Parameters
fixed_edge
corner_nodes
debug
Returns
MoFEMErrorCode

Definition at line 2102 of file CutMeshInterface.cpp.

2104  {
2105  CoreInterface &m_field = cOre;
2106  moab::Interface &moab = m_field.get_moab();
2108 
2109  Skinner skin(&moab);
2110  Range trim_tets_skin;
2111  CHKERR skin.find_skin(0, trimNewVolumes, false, trim_tets_skin);
2112  Range trim_tets_skin_edges;
2113  CHKERR moab.get_adjacencies(trim_tets_skin, 1, false, trim_tets_skin_edges,
2114  moab::Interface::UNION);
2115 
2116  Range barrier_vertices(trimNewVertices);
2117 
2118  if (fixed_edges && trimFixedEdges) {
2119 
2120  // get all vertices on fixed edges and surface
2121  Range trim_surface_edges;
2122  CHKERR moab.get_adjacencies(trimNewSurfaces, 1, false, trim_surface_edges,
2123  moab::Interface::UNION);
2124  Range fixed_edges_on_trim_surface;
2125  fixed_edges_on_trim_surface = intersect(trim_surface_edges, *fixed_edges);
2126  Range fixed_edges_on_trim_surface_verts;
2127  CHKERR moab.get_connectivity(fixed_edges_on_trim_surface,
2128  fixed_edges_on_trim_surface_verts, false);
2129 
2130  // get faces adjacent to barrier_vertices
2131  Range barrier_vertices_faces;
2132  CHKERR moab.get_adjacencies(barrier_vertices, 2, false,
2133  barrier_vertices_faces, moab::Interface::UNION);
2134  barrier_vertices_faces = intersect(barrier_vertices_faces, trimNewSurfaces);
2135 
2136  // get vertices on fixed edges
2137  Range fixed_edges_vertices;
2138  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_vertices, false);
2139  fixed_edges_vertices = intersect(barrier_vertices, fixed_edges_vertices);
2140  fixed_edges_vertices =
2141  subtract(fixed_edges_vertices, fixed_edges_on_trim_surface_verts);
2142  if (corner_nodes)
2143  fixed_edges_vertices.merge(intersect(barrier_vertices, *corner_nodes));
2144 
2145  // get faces adjacent to vertices on fixed edges
2146  Range fixed_edges_faces;
2147  CHKERR moab.get_adjacencies(fixed_edges_vertices, 2, false,
2148  fixed_edges_faces, moab::Interface::UNION);
2149  fixed_edges_faces = intersect(fixed_edges_faces, barrier_vertices_faces);
2150 
2151  if (debug && !fixed_edges_faces.empty())
2152  CHKERR SaveData(m_field.get_moab())("fixed_edges_faces.vtk",
2153  fixed_edges_faces);
2154 
2155  // get nodes on faces
2156  Range fixed_edges_faces_vertices;
2157  CHKERR moab.get_connectivity(fixed_edges_faces, fixed_edges_faces_vertices,
2158  false);
2159  barrier_vertices.merge(fixed_edges_faces_vertices);
2160  }
2161 
2162  auto remove_faces_on_skin = [&]() {
2164  Range skin_faces = intersect(trimNewSurfaces, trim_tets_skin);
2165  if (!skin_faces.empty()) {
2166  Range add_to_barrier;
2167  CHKERR moab.get_connectivity(skin_faces, add_to_barrier, false);
2168  barrier_vertices.merge(add_to_barrier);
2169  for (auto f : skin_faces)
2170  trimNewSurfaces.erase(f);
2171  }
2173  };
2174 
2175  auto get_trim_free_edges = [&]() {
2176  // get current surface skin
2177  Range trim_surf_skin;
2178  CHKERR skin.find_skin(0, trimNewSurfaces, false, trim_surf_skin);
2179  trim_surf_skin = subtract(trim_surf_skin, trim_tets_skin_edges);
2180 
2181  Range trim_surf_skin_verts;
2182  CHKERR moab.get_connectivity(trim_surf_skin, trim_surf_skin_verts, true);
2183  for (auto e : barrier_vertices)
2184  trim_surf_skin_verts.erase(e);
2185 
2186  Range free_edges;
2187  CHKERR moab.get_adjacencies(trim_surf_skin_verts, 1, false, free_edges,
2188  moab::Interface::UNION);
2189  free_edges = intersect(free_edges, trim_surf_skin);
2190 
2191  return free_edges;
2192  };
2193 
2194  CHKERR remove_faces_on_skin();
2195 
2196  if (debug && !barrier_vertices.empty())
2197  CHKERR SaveData(m_field.get_moab())("barrier_vertices.vtk",
2198  barrier_vertices);
2199 
2200  int nn = 0;
2201  Range out_edges;
2202  while (!(out_edges = get_trim_free_edges()).empty()) {
2203 
2204  if (debug && !trimNewSurfaces.empty())
2205  CHKERR SaveData(m_field.get_moab())(
2206  "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
2207  trimNewSurfaces);
2208 
2209  if (debug && !out_edges.empty())
2210  CHKERR SaveData(m_field.get_moab())(
2211  "trimNewSurfacesOutsideVerts_" +
2212  boost::lexical_cast<std::string>(nn) + ".vtk",
2213  out_edges);
2214 
2215  Range outside_faces;
2216  CHKERR moab.get_adjacencies(out_edges, 2, false, outside_faces,
2217  moab::Interface::UNION);
2218  trimNewSurfaces = subtract(trimNewSurfaces, outside_faces);
2219  ++nn;
2220  }
2221 
2222  if (debug && !trimNewSurfaces.empty())
2223  CHKERR SaveData(m_field.get_moab())(
2224  "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
2225  trimNewSurfaces);
2226 
2228 }
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

Member Data Documentation

◆ aveLength

double MoFEM::CutMeshInterface::aveLength
private

Average edge length.

Definition at line 577 of file CutMeshInterface.hpp.

◆ cOre

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

Definition at line 37 of file CutMeshInterface.hpp.

◆ cutEdges

Range MoFEM::CutMeshInterface::cutEdges
private

Definition at line 550 of file CutMeshInterface.hpp.

◆ cutFrontVolumes

Range MoFEM::CutMeshInterface::cutFrontVolumes
private

Definition at line 581 of file CutMeshInterface.hpp.

◆ cutNewSurfaces

Range MoFEM::CutMeshInterface::cutNewSurfaces
private

Definition at line 552 of file CutMeshInterface.hpp.

◆ cutNewVertices

Range MoFEM::CutMeshInterface::cutNewVertices
private

Definition at line 555 of file CutMeshInterface.hpp.

◆ cutNewVolumes

Range MoFEM::CutMeshInterface::cutNewVolumes
private

Definition at line 551 of file CutMeshInterface.hpp.

◆ cutSurfaceVolumes

Range MoFEM::CutMeshInterface::cutSurfaceVolumes
private

Definition at line 580 of file CutMeshInterface.hpp.

◆ edgesToCut

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

Definition at line 572 of file CutMeshInterface.hpp.

◆ edgesToTrim

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

Definition at line 574 of file CutMeshInterface.hpp.

◆ fRont

Range MoFEM::CutMeshInterface::fRont
private

Definition at line 541 of file CutMeshInterface.hpp.

◆ lineSearchSteps

int MoFEM::CutMeshInterface::lineSearchSteps

Definition at line 41 of file CutMeshInterface.hpp.

◆ maxLength

double MoFEM::CutMeshInterface::maxLength
private

Maximal edge length.

Definition at line 578 of file CutMeshInterface.hpp.

◆ mergedSurfaces

Range MoFEM::CutMeshInterface::mergedSurfaces
private

Definition at line 563 of file CutMeshInterface.hpp.

◆ mergedVolumes

Range MoFEM::CutMeshInterface::mergedVolumes
private

Definition at line 562 of file CutMeshInterface.hpp.

◆ nbMaxMergingCycles

int MoFEM::CutMeshInterface::nbMaxMergingCycles

Definition at line 42 of file CutMeshInterface.hpp.

◆ projectEntitiesQualityTrashold

double MoFEM::CutMeshInterface::projectEntitiesQualityTrashold

Definition at line 43 of file CutMeshInterface.hpp.

◆ rootSetSurf

EntityHandle MoFEM::CutMeshInterface::rootSetSurf
private

Definition at line 548 of file CutMeshInterface.hpp.

◆ sUrface

Range MoFEM::CutMeshInterface::sUrface
private

Definition at line 542 of file CutMeshInterface.hpp.

◆ treeSurfPtr

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

Definition at line 547 of file CutMeshInterface.hpp.

◆ trimEdges

Range MoFEM::CutMeshInterface::trimEdges
private

Definition at line 560 of file CutMeshInterface.hpp.

◆ trimFixedEdges

bool MoFEM::CutMeshInterface::trimFixedEdges
private

Definition at line 545 of file CutMeshInterface.hpp.

◆ trimNewSurfaces

Range MoFEM::CutMeshInterface::trimNewSurfaces
private

Definition at line 559 of file CutMeshInterface.hpp.

◆ trimNewVertices

Range MoFEM::CutMeshInterface::trimNewVertices
private

Definition at line 558 of file CutMeshInterface.hpp.

◆ trimNewVolumes

Range MoFEM::CutMeshInterface::trimNewVolumes
private

Definition at line 557 of file CutMeshInterface.hpp.

◆ verticesOnCutEdges

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

Definition at line 573 of file CutMeshInterface.hpp.

◆ verticesOnTrimEdges

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

Definition at line 575 of file CutMeshInterface.hpp.

◆ vOlume

Range MoFEM::CutMeshInterface::vOlume
private

Definition at line 543 of file CutMeshInterface.hpp.

◆ zeroDistanceEnts

Range MoFEM::CutMeshInterface::zeroDistanceEnts
private

Definition at line 553 of file CutMeshInterface.hpp.

◆ zeroDistanceVerts

Range MoFEM::CutMeshInterface::zeroDistanceVerts
private

Definition at line 554 of file CutMeshInterface.hpp.


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