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 (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 low_tol, int verb=QUIET, const bool debug=false)
 find edges to cut More...
 
MoFEMErrorCode projectZeroDistanceEnts (Range *fixed_edges, Range *corner_nodes, const double low_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 & getCutVolumes () 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 cutVolumes
 
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 529 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 229 of file CutMeshInterface.cpp.

229  {
230  CoreInterface &m_field = cOre;
231  moab::Interface &moab = m_field.get_moab();
233  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
234  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
237 }
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:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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:476
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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:77
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:121
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:104
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ createFrontLevelSets()

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

Calculate distance from mesh nodes to surface front.

Parameters
verb
debug
Returns
MoFEMErrorCode

Definition at line 547 of file CutMeshInterface.cpp.

548  {
549  CoreInterface &m_field = cOre;
550  moab::Interface &moab = m_field.get_moab();
552 
553  auto create_tag = [&](const std::string name, const int dim) {
554  Tag th;
555  rval = moab.tag_get_handle(name.c_str(), th);
556  if (rval == MB_SUCCESS)
557  return th;
558  std::vector<double> def_val(dim, 0);
559  CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
560  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
561  return th;
562  };
563 
564  Range vol_vertices;
565  CHKERR moab.get_connectivity(vOlume, vol_vertices, true);
566 
567  std::vector<double> min_distances_from_front(vol_vertices.size(), -1);
568  std::vector<double> points_on_edges(3 * vol_vertices.size(), 0);
569  std::vector<EntityHandle> closest_edges(vol_vertices.size(), 0);
570 
571  std::vector<double> coords(3 * vol_vertices.size());
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:476
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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 < 0.5 * 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:476
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:91
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:104
#define CHKERR
Inline error check.
Definition: definitions.h:595
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:88
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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 
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 }
MoFEMErrorCode createFrontLevelSets(int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to surface front.
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:476
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:595
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:406

◆ createSurfaceLevelSets()

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

Calculate distance from mesh nodes to cut surface.

Parameters
intersect_vol
verb
debug
Returns
MoFEMErrorCode

Definition at line 477 of file CutMeshInterface.cpp.

478  {
479  CoreInterface &m_field = cOre;
480  moab::Interface &moab = m_field.get_moab();
482 
483  auto create_tag = [&](const std::string name, const int dim) {
484  Tag th;
485  rval = moab.tag_get_handle(name.c_str(), th);
486  if (rval == MB_SUCCESS)
487  return th;
488  std::vector<double> def_val(dim, 0);
489  CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
490  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
491 
492  return th;
493  };
494 
495  auto set_vol = [&](const Range &vol_verts, std::vector<double> &coords,
496  std::vector<double> &dist_surface_vec,
497  std::vector<double> &dist_surface_normal_vec) {
499 
500  coords.resize(3 * vol_verts.size());
501  dist_surface_vec.resize(3 * vol_verts.size());
502  dist_surface_normal_vec.resize(3 * vol_verts.size());
503  CHKERR moab.get_coords(vol_verts, &*coords.begin());
504  for (auto v : vol_verts) {
505 
506  const int index = vol_verts.index(v);
507  auto point_in = getVectorAdaptor(&coords[3 * index], 3);
508  VectorDouble3 point_out(3);
509  EntityHandle facets_out;
510  CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
511  &point_out[0], facets_out);
512 
513  VectorDouble3 delta = point_out - point_in;
514  auto dist_vec = getVectorAdaptor(&dist_surface_vec[3 * index], 3);
515  noalias(dist_vec) = delta;
516 
517  VectorDouble3 n(3);
518  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
519  n /= norm_2(n);
520  auto dist_normal_vec =
521  getVectorAdaptor(&dist_surface_normal_vec[3 * index], 3);
522  noalias(dist_normal_vec) = inner_prod(delta, n) * n;
523 
524  }
525 
527  };
528 
529  std::vector<double> coords;
530  std::vector<double> dist_surface_vec;
531  std::vector<double> dist_surface_normal_vec;
532  Range vol_verts;
533  CHKERR moab.get_connectivity(vOlume, vol_verts, true);
534 
535  CHKERR set_vol(vol_verts, coords, dist_surface_vec, dist_surface_normal_vec);
536 
537  auto th_dist_surface_vec = create_tag("DIST_SURFACE_VECTOR", 3);
538  auto th_dist_surface_normal_vec = create_tag("DIST_SURFACE_NORMAL_VECTOR", 3);
539  CHKERR moab.tag_set_data(th_dist_surface_vec, vol_verts,
540  &*dist_surface_vec.begin());
541  CHKERR moab.tag_set_data(th_dist_surface_normal_vec, vol_verts,
542  &*dist_surface_normal_vec.begin());
543 
545 }
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:476
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:595
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:88
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

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

321  {
322  CoreInterface &m_field = cOre;
323  moab::Interface &moab = m_field.get_moab();
325 
326  std::vector<BitRefLevel> bit_levels;
327 
328  auto get_back_bit_levels = [&]() {
329  bit_levels.push_back(BitRefLevel());
330  bit_levels.back().set(first_bit + bit_levels.size() - 1);
331  return bit_levels.back();
332  };
333 
334  BitRefLevel cut_bit = get_back_bit_levels();
335 
336  CHKERR cutOnly(unite(cutSurfaceVolumes, cutFrontVolumes), cut_bit, th,
337  tol_cut, tol_cut_close, fixed_edges, corner_nodes,
338  update_meshsets, debug);
339 
340  auto get_min_quality = [&m_field](const BitRefLevel bit, Tag th) {
341  Range tets_level; // test at level
342  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
343  bit, BitRefLevel().set(), MBTET, tets_level);
344  double min_q = 1;
345  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
346  return min_q;
347  };
348 
349  PetscPrintf(PETSC_COMM_WORLD, "Min quality cut %6.4g\n",
350  get_min_quality(cut_bit, th));
351 
352  Range starting_volume = cutNewVolumes;
353  Range starting_volume_nodes;
354  CHKERR m_field.get_moab().get_connectivity(starting_volume,
355  starting_volume_nodes, true);
356  std::vector<double> staring_volume_coords(3 * starting_volume_nodes.size());
357  CHKERR m_field.get_moab().get_coords(starting_volume_nodes,
358  &*staring_volume_coords.begin());
359 
360  if (th) {
361  std::vector<double> staring_volume_th_coords(3 *
362  starting_volume_nodes.size());
363  CHKERR m_field.get_moab().tag_get_data(th, starting_volume_nodes,
364  &*staring_volume_th_coords.begin());
365  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
366  &*staring_volume_th_coords.begin());
367  }
368 
369  if (th)
370  CHKERR setTagData(th);
371 
372  BitRefLevel trim_bit = get_back_bit_levels();
373 
374  CHKERR trimOnly(trim_bit, th, tol_trim_close, fixed_edges, corner_nodes,
375  update_meshsets, debug);
376 
377  PetscPrintf(PETSC_COMM_WORLD, "Min quality trim %3.2g\n",
378  get_min_quality(trim_bit, th));
379 
380  first_bit += bit_levels.size() - 1;
381 
382  if (th)
383  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
384  &*staring_volume_coords.begin());
385 
387 }
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
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:476
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:595
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:406
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 1306 of file CutMeshInterface.cpp.

1310  {
1311  CoreInterface &m_field = cOre;
1312  moab::Interface &moab = m_field.get_moab();
1313  MeshRefinement *refiner;
1314  const RefEntity_multiIndex *ref_ents_ptr;
1316 
1317  if (cutEdges.size() != edgesToCut.size())
1318  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1319 
1320  CHKERR m_field.getInterface(refiner);
1321  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1322  CHKERR refiner->add_vertices_in_the_middel_of_edges(cutEdges, bit);
1323  CHKERR refiner->refine_TET(vOlume, bit, false, QUIET,
1324  debug ? &cutEdges : NULL);
1325 
1326  if (debug)
1327  if (cutEdges.size() != edgesToCut.size())
1328  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1329 
1330  cut_vols.clear();
1331  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1332  bit, BitRefLevel().set(), MBTET, cut_vols);
1333  cut_surf.clear();
1334  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1335  bit, BitRefLevel().set(), MBTRI, cut_surf);
1336 
1337  // Find new vertices on cut edges
1338  cut_verts.clear();
1339  CHKERR moab.get_connectivity(zeroDistanceEnts, cut_verts, true);
1340  cut_verts.merge(zeroDistanceVerts);
1341  for (auto &m : edgesToCut) {
1342  RefEntity_multiIndex::index<
1343  Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vit =
1344  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1345  boost::make_tuple(MBVERTEX, m.first));
1346  if (vit ==
1347  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1348  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1349  "No vertex on cut edges, that make no sense");
1350  }
1351  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1352  if ((ref_ent->getBitRefLevel() & bit).any()) {
1353  EntityHandle vert = ref_ent->getRefEnt();
1354  cut_verts.insert(vert);
1355  verticesOnCutEdges[vert] = m.second;
1356  } else {
1357  SETERRQ1(
1358  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1359  "Vertex has wrong bit ref level %s",
1360  boost::lexical_cast<std::string>(ref_ent->getBitRefLevel()).c_str());
1361  }
1362  }
1363 
1364  // Add zero distance entities faces
1365  Range tets_skin;
1366  Skinner skin(&moab);
1367  CHKERR skin.find_skin(0, cut_vols, false, tets_skin);
1368  cut_surf.merge(zeroDistanceEnts.subset_by_type(MBTRI));
1369 
1370  // At that point cut_surf has all newly created faces, now take all
1371  // nodes on those faces and subtract nodes on cut edges. Faces adjacent to
1372  // nodes which left are not part of surface.
1373  Range diff_verts;
1374  CHKERR moab.get_connectivity(unite(cut_surf, zeroDistanceEnts), diff_verts,
1375  true);
1376  diff_verts = subtract(diff_verts, cut_verts);
1377  Range subtract_faces;
1378  CHKERR moab.get_adjacencies(diff_verts, 2, false, subtract_faces,
1379  moab::Interface::UNION);
1380  cut_surf = subtract(cut_surf, unite(subtract_faces, tets_skin));
1381  cut_verts.clear();
1382  CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1383 
1385 }
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:476
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:595
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:406
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 240 of file CutMeshInterface.cpp.

243  {
244  CoreInterface &m_field = cOre;
245  moab::Interface &moab = m_field.get_moab();
247 
248  // cut mesh
249  CHKERR findEdgesToCut(vol, fixed_edges, corner_nodes, tol_cut, QUIET, debug);
250  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_cut_close,
251  QUIET, debug);
254  if (fixed_edges)
255  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
256  *fixed_edges);
257  if (corner_nodes)
258  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
259  *corner_nodes);
260  if (update_meshsets)
261  CHKERR m_field.getInterface<MeshsetsManager>()
262  ->updateAllMeshsetsByEntitiesChildren(cut_bit);
264 
265  if (debug) {
267  if(fixed_edges)
268  CHKERR SaveData(moab)("out_cut_new_fixed_edges.vtk", *fixed_edges);
269  }
270 
272 }
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:476
MoFEMErrorCode cutEdgesInMiddle(const BitRefLevel bit, Range &cut_vols, Range &cut_surf, Range &cut_verts, const bool debug=false)
cut edges
MoFEMErrorCode findEdgesToCut(Range vol, Range *fixed_edges, Range *corner_nodes, const double low_tol, int verb=QUIET, const bool debug=false)
find edges to cut
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

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

392  {
393  CoreInterface &m_field = cOre;
395 
396  std::vector<BitRefLevel> bit_levels;
397 
398  auto get_back_bit_levels = [&]() {
399  bit_levels.push_back(BitRefLevel());
400  bit_levels.back().set(first_bit + bit_levels.size() - 1);
401  return bit_levels.back();
402  };
403 
404  if (debug) {
405  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
406  "ents_not_in_database.vtk", "VTK", "");
407  }
408 
409  CHKERR cutAndTrim(first_bit, th, tol_cut, tol_cut_close, tol_trim_close,
410  &fixed_edges, &corner_nodes, update_meshsets, debug);
411  if (debug) {
412  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
413  "cut_trim_ents_not_in_database.vtk", "VTK", "");
414  }
415 
416  BitRefLevel bit_level1 = BitRefLevel().set(first_bit - 1);
417  BitRefLevel bit_level2 = get_back_bit_levels();
418  BitRefLevel bit_level3 = get_back_bit_levels();
419 
420  CHKERR mergeBadEdges(fraction_level, bit_level1, bit_level2, bit_level3,
421  getNewTrimSurfaces(), fixed_edges, corner_nodes, th,
422  update_meshsets, debug);
423 
424  auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
425  Range tets_level; // test at level
426  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
427  bit, BitRefLevel().set(), MBTET, tets_level);
428  double min_q = 1;
429  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
430  if (min_q < 0 && debug) {
431  CHKERR m_field.getInterface<Tools>()->writeTetsWithQuality(
432  "negative_tets.vtk", "VTK", "", tets_level, th);
433  }
434  return min_q;
435  };
436 
437  PetscPrintf(PETSC_COMM_WORLD, "Min quality node merge %6.4g\n",
438  get_min_quality(bit_level3, th));
439 
440  CHKERR cOre.getInterface<BitRefManager>()->updateRange(fixed_edges,
441  fixed_edges);
442  CHKERR cOre.getInterface<BitRefManager>()->updateRange(corner_nodes,
443  corner_nodes);
444 
445  first_bit += bit_levels.size() - 1;
446 
447  if (debug) {
448  CHKERR cOre.getInterface<BitRefManager>()->writeBitLevelByType(
449  bit_level3, BitRefLevel().set(), MBTET, "out_tets_merged.vtk", "VTK",
450  "");
451  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
452  "cut_trim_merge_ents_not_in_database.vtk", "VTK", "");
453  CHKERR SaveData(m_field.get_moab())("merged_surface.vtk", mergedSurfaces);
454  }
455 
457 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:595
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:406

◆ findEdgesToCut()

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

find edges to cut

Parameters
verbverbosity level
Returns
error code

Definition at line 802 of file CutMeshInterface.cpp.

805  {
806  CoreInterface &m_field = cOre;
807  moab::Interface &moab = m_field.get_moab();
809 
810  edgesToCut.clear();
811  cutEdges.clear();
812 
813  zeroDistanceVerts.clear();
814  zeroDistanceEnts.clear();
815  verticesOnCutEdges.clear();
816 
817  Tag th_dist_normal;
818  CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_VECTOR", th_dist_normal);
819 
820  auto not_project_node = [this, &moab](const EntityHandle v) {
822  VectorDouble3 s0(3);
823  CHKERR moab.get_coords(&v, 1, &s0[0]);
824  verticesOnCutEdges[v].dIst = 0;
825  verticesOnCutEdges[v].lEngth = 0;
826  verticesOnCutEdges[v].unitRayDir.resize(3, false);
827  verticesOnCutEdges[v].unitRayDir.clear();
828  verticesOnCutEdges[v].rayPoint = s0;
830  };
831 
832  auto get_ave_edge_length = [&](const EntityHandle ent,
833  const Range &vol_edges) {
834  Range adj_edges;
835  if (moab.type_from_handle(ent) == MBVERTEX) {
836  CHKERR moab.get_adjacencies(&ent, 1, 1, false, adj_edges);
837  } else {
838  adj_edges.insert(ent);
839  }
840  adj_edges = intersect(adj_edges, vol_edges);
841  double ave_l = 0;
842  for (auto e : adj_edges) {
843  int num_nodes;
844  const EntityHandle *conn;
845  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
846  VectorDouble6 coords(6);
847  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
849  &coords[0], &coords[1], &coords[2]);
851  &coords[3], &coords[4], &coords[5]);
853  t_n0(i) -= t_n1(i);
854  ave_l += sqrt(t_n0(i) * t_n0(i));
855  }
856  return ave_l / adj_edges.size();
857  };
858 
859  auto get_tag_data = [&](auto th, auto conn) {
860  const void *ptr;
861  CHKERR moab.tag_get_by_ptr(th, &conn, 1, &ptr);
862  return getVectorAdaptor(
863  const_cast<double *>(static_cast<const double *>(ptr)), 3);
864  };
865 
866  Range vol_edges;
867  CHKERR moab.get_adjacencies(vol, 1, true, vol_edges, moab::Interface::UNION);
868 
869  aveLength = 0;
870  maxLength = 0;
871  int nb_ave_length = 0;
872  for (auto e : vol_edges) {
873 
874  int num_nodes;
875  const EntityHandle *conn;
876  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
877  const double tol = get_ave_edge_length(e, vol_edges) * low_tol;
878 
879  VectorDouble6 coords(6);
880  CHKERR moab.get_coords(conn, 2, &coords[0]);
881  VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
882  VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
883  VectorDouble3 ray = n1 - n0;
884  const double ray_length = norm_2(ray);
885  ray /= ray_length;
886 
887  auto dist_vec0 = get_tag_data(th_dist_normal, conn[0]);
888  auto dist_vec1 = get_tag_data(th_dist_normal, conn[1]);
889 
890  const double dist_normal0 = norm_2(dist_vec0);
891  const double dist_normal1 = norm_2(dist_vec1);
892 
893  // check if cutting point is close to the end of the edges
894  if (dist_normal0 < tol && dist_normal1 < tol) {
895 
896  aveLength += norm_2(ray);
897  maxLength = fmax(maxLength, norm_2(ray));
898  ++nb_ave_length;
899 
900  for (int n : {0, 1})
901  CHKERR not_project_node(conn[n]);
902  zeroDistanceEnts.insert(e);
903 
904  } else if (inner_prod(dist_vec0, dist_vec1) < 0) {
905 
906  // Edges is on two sides of the surface
907 
908  const double s = dist_normal0 / (dist_normal1 + dist_normal0);
909  edgesToCut[e].dIst = s * ray_length;
910  edgesToCut[e].lEngth = ray_length;
911  edgesToCut[e].unitRayDir = ray;
912  edgesToCut[e].rayPoint = n0;
913  cutEdges.insert(e);
914 
915  aveLength += norm_2(ray);
916  maxLength = fmax(maxLength, norm_2(ray));
917  ++nb_ave_length;
918  }
919  }
920  aveLength /= nb_ave_length;
921 
922  Range vol_vertices;
923  CHKERR moab.get_connectivity(vol, vol_vertices, true);
924  for (auto v : vol_vertices) {
925 
926  VectorDouble3 dist_normal(3);
927  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, &*dist_normal.begin());
928  const double dist = norm_2(dist_normal);
929 
930  const double tol = get_ave_edge_length(v, vol_edges) * low_tol;
931  if (dist < tol) {
932  CHKERR not_project_node(v);
933  zeroDistanceVerts.insert(v);
934  }
935  }
936 
937  cutVolumes.clear();
938  CHKERR moab.get_adjacencies(unite(cutEdges, zeroDistanceVerts), 3, false,
939  cutVolumes, moab::Interface::UNION);
940  cutVolumes = intersect(cutVolumes, vOlume);
941 
942  for (auto v : zeroDistanceVerts) {
943  Range adj_edges;
944  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges,
945  moab::Interface::UNION);
946  for (auto e : adj_edges) {
947  cutEdges.erase(e);
948  edgesToCut.erase(e);
949  }
950  }
951 
952  if (debug)
953  CHKERR SaveData(m_field.get_moab())("cut_edges.vtk", cutEdges);
954 
955  if (debug)
956  CHKERR SaveData(m_field.get_moab())(
957  "cut_zero_dist.vtk", unite(zeroDistanceVerts, zeroDistanceEnts));
958 
960 }
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:500
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
double aveLength
Average edge length.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:91
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:104
#define CHKERR
Inline error check.
Definition: definitions.h:595
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:88
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:406
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 1422 of file CutMeshInterface.cpp.

1425  {
1426  CoreInterface &m_field = cOre;
1427  moab::Interface &moab = m_field.get_moab();
1429 
1430  // takes body skin
1431  Skinner skin(&moab);
1432  Range tets_skin;
1433  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
1434 
1435  // vertices on the skin
1436  Range tets_skin_verts;
1437  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1438  // edges on the skin
1439  Range tets_skin_edges;
1440  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1441  moab::Interface::UNION);
1442  // get edges on new surface
1443  Range edges;
1444  CHKERR moab.get_adjacencies(cutNewSurfaces, 1, false, edges,
1445  moab::Interface::UNION);
1446  Range corners;
1447  if (corner_nodes)
1448  corners = *corner_nodes;
1449 
1450  Range fix_edges;
1451  if (fixed_edges)
1452  fix_edges = *fixed_edges;
1453 
1454  Range fixed_edges_verts;
1455  if (fixed_edges)
1456  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts, true);
1457 
1458  Range surface_skin;
1459  if (fRont.empty())
1460  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
1461  else
1462  surface_skin = fRont;
1463 
1464  auto get_point_coords = [&](EntityHandle v) {
1465  VectorDouble3 point(3);
1466  if (th)
1467  CHKERR moab.tag_get_data(th, &v, 1, &point[0]);
1468  else
1469  CHKERR moab.get_coords(&v, 1, &point[0]);
1470  return point;
1471  };
1472 
1473  struct VertMap {
1474  const double d;
1475  const EntityHandle v;
1476  const EntityHandle e;
1477  VertMap(const double d, const EntityHandle v, const EntityHandle e)
1478  : d(d), v(v), e(e) {}
1479  };
1480 
1481  typedef multi_index_container<
1482  VertMap,
1483  indexed_by<
1484  ordered_non_unique<member<VertMap, const double, &VertMap::d>>,
1485  ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::v>>,
1486  ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::e>>
1487 
1488  >>
1489  VertMapMultiIndex;
1490 
1491  VertMapMultiIndex verts_map;
1492 
1493  auto add_vert = [&](const EntityHandle v, const EntityHandle e,
1494  const double dist) {
1495  verts_map.insert(VertMap(dist, v, e));
1496  };
1497 
1498  // clear data ranges
1499  trimEdges.clear();
1500  edgesToTrim.clear();
1501  verticesOnTrimEdges.clear();
1502  trimNewVertices.clear();
1503 
1504  if (debug)
1505  CHKERR SaveData(moab)("edges_potentially_to_trim.vtk", edges);
1506 
1507  if (debug)
1508  CHKERR SaveData(moab)("surface_skin_to_trim.vtk", surface_skin);
1509 
1510  auto cut_this_edge = [&](const EntityHandle e, const double length, auto &ray,
1511  auto &ray_point) {
1512  trimEdges.insert(e);
1513  edgesToTrim[e].dIst = 1;
1514  edgesToTrim[e].lEngth = length;
1515  edgesToTrim[e].unitRayDir.resize(3, false);
1516  edgesToTrim[e].rayPoint.resize(3, false);
1517  for (int ii = 0; ii != 3; ++ii)
1518  edgesToTrim[e].unitRayDir[ii] = ray(ii);
1519  for (int ii = 0; ii != 3; ++ii)
1520  edgesToTrim[e].rayPoint[ii] = ray_point(ii);
1521  };
1522 
1524  int num_nodes;
1525 
1526  // iterate over edges on cut surface
1527  for (auto e : edges) {
1528 
1529  // Get edge connectivity and coords
1530  const EntityHandle *conn_edge;
1531  CHKERR moab.get_connectivity(e, conn_edge, num_nodes, true);
1532  double coords_edge[3 * num_nodes];
1533  CHKERR moab.get_coords(conn_edge, num_nodes, coords_edge);
1534 
1535  FTensor::Tensor1<double *, 3> t_e0(&coords_edge[0], &coords_edge[1],
1536  &coords_edge[2]);
1537  FTensor::Tensor1<double *, 3> t_e1(&coords_edge[3], &coords_edge[4],
1538  &coords_edge[5]);
1539 
1540  FTensor::Tensor1<double, 3> t_edge_delta;
1541  t_edge_delta(i) = t_e1(i) - t_e0(i);
1542  const double edge_length2 = t_edge_delta(i) * t_edge_delta(i);
1543  const double edge_length = sqrt(edge_length2);
1544  if (edge_length == 0)
1545  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Zero edge length");
1546 
1547  auto get_edge_coors = [&](const auto e) {
1548  const EntityHandle *conn;
1549  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1550  VectorDouble6 coords(6);
1551  if (th)
1552  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1553  else
1554  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1555  return coords;
1556  };
1557 
1558  // iterate entities on inner surface font to find that edge cross
1559  for (auto s : surface_skin) {
1560 
1561  auto coords_front = get_edge_coors(s);
1562 
1563  FTensor::Tensor1<double *, 3> t_f0(&coords_front[0], &coords_front[1],
1564  &coords_front[2]);
1565  FTensor::Tensor1<double *, 3> t_f1(&coords_front[3], &coords_front[4],
1566  &coords_front[5]);
1567 
1568  // find point of minilam distance between front and cut surface edge
1569  double t_edge = -1, t_front = -1;
1570  auto res = Tools::minDistanceFromSegments(&t_e0(0), &t_e1(0), &t_f0(0),
1571  &t_f1(0), &t_edge, &t_front);
1572 
1573  if (res != Tools::NO_SOLUTION) {
1574 
1575  // check if edges crossing each other in the middle (it not imply that
1576  // have common point)
1577  const double overlap_tol = 1e-2;
1578  if (t_edge > -std::numeric_limits<float>::epsilon() &&
1579  t_edge < 1 + std::numeric_limits<float>::epsilon() &&
1580  t_front >= -overlap_tol && t_front <= 1 + overlap_tol) {
1581 
1582  FTensor::Tensor1<double, 3> t_front_delta;
1583  t_front_delta(i) = t_f1(i) - t_f0(i);
1584  FTensor::Tensor1<double, 3> t_edge_delta;
1585  t_edge_delta(i) = t_e1(i) - t_e0(i);
1586  FTensor::Tensor1<double, 3> t_edge_point, t_front_point;
1587  t_edge_point(i) = t_e0(i) + t_edge * t_edge_delta(i);
1588  t_front_point(i) = t_f0(i) + t_front * t_front_delta(i);
1590  t_ray(i) = t_front_point(i) - t_edge_point(i);
1591  const double dist = sqrt(t_ray(i) * t_ray(i));
1592 
1593  // that imply that edges have common point
1594  if ((dist / edge_length) < 1.) {
1595 
1596  auto check_to_add_edge = [&](const EntityHandle e,
1597  const double dist) {
1598  auto eit = edgesToTrim.find(e);
1599  if (eit != edgesToTrim.end())
1600  if (eit->second.dIst < dist)
1601  return false;
1602  return true;
1603  };
1604 
1606  if (t_edge < 0.5) {
1607  t_ray(i) = t_edge * t_edge_delta(i);
1608  const double ray_length = sqrt(t_ray(i) * t_ray(i));
1609  if (check_to_add_edge(e, ray_length)) {
1610  add_vert(conn_edge[0], e, fabs(t_edge));
1611  add_vert(conn_edge[1], e, fabs(t_edge - 1));
1612  cut_this_edge(e, edge_length, t_ray, t_e0);
1613  }
1614  } else {
1615  FTensor::Tensor1<double, 3> t_edge_point;
1616  t_edge_point(i) = t_e0(i) + t_edge * t_edge_delta(i);
1617  t_ray(i) = t_edge_point(i) - t_e1(i);
1618  const double ray_length = sqrt(t_ray(i) * t_ray(i));
1619  if (check_to_add_edge(e, ray_length)) {
1620  add_vert(conn_edge[0], e, fabs(t_edge));
1621  add_vert(conn_edge[1], e, fabs(t_edge - 1));
1622  cut_this_edge(e, edge_length, t_ray, t_e1);
1623  }
1624  }
1625  }
1626  }
1627  }
1628 
1629  }
1630  }
1631 
1632  if (debug)
1633  CHKERR SaveData(moab)("edges_selected_to_trim.vtk", trimEdges);
1634 
1635  auto get_quality_change = [&](const Range &adj_tets, const EntityHandle &v,
1636  const VectorDouble3 &new_pos) {
1637  double q0 = 1;
1638  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0, th);
1639 
1640  double min_q = 1;
1641  for (auto t : adj_tets) {
1642  int num_nodes;
1643  const EntityHandle *conn;
1644  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1645  VectorDouble12 coords(12);
1646  if (th)
1647  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1648  else
1649  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1650 
1651  for (int n = 0; n != 4; ++n) {
1652  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1653  if (conn[n] == v) {
1654  noalias(n_coords) = new_pos;
1655  } else {
1656  auto m = verticesOnTrimEdges.find(conn[n]);
1657  if (m != verticesOnTrimEdges.end())
1658  noalias(n_coords) =
1659  m->second.rayPoint + m->second.dIst * m->second.unitRayDir;
1660  }
1661  }
1662 
1663  const double q = Tools::volumeLengthQuality(&coords[0]);
1664  if (!std::isnormal(q)) {
1665  min_q = -2;
1666  break;
1667  }
1668  min_q = std::min(min_q, q);
1669  }
1670  return min_q / q0;
1671  };
1672 
1673  auto add_trim_vert = [&](const EntityHandle v, const EntityHandle e) {
1674  if (!(trimNewVertices.find(v) != trimNewVertices.end())) {
1675  auto &r = edgesToTrim.at(e);
1676  VectorDouble3 ray_point = get_point_coords(v);
1677  VectorDouble3 new_pos = r.rayPoint + r.dIst * r.unitRayDir;
1678  VectorDouble3 unit_ray_dir = ray_point - new_pos;
1679  Range adj_tets;
1680  CHKERR moab.get_adjacencies(&v, 1, 3, false, adj_tets);
1681  adj_tets = intersect(adj_tets, cutNewVolumes);
1682  double q = get_quality_change(adj_tets, v, new_pos);
1684  VectorDouble3 ray_dir = new_pos - ray_point;
1685  double dist = norm_2(ray_dir);
1686  verticesOnTrimEdges[v].dIst = 1;
1687  verticesOnTrimEdges[v].lEngth = dist;
1688  verticesOnTrimEdges[v].unitRayDir = ray_dir;
1689  verticesOnTrimEdges[v].rayPoint = ray_point;
1690  trimNewVertices.insert(v);
1691  }
1692  }
1693  };
1694 
1695  auto add_no_move_trim = [&](const EntityHandle v, const EntityHandle e) {
1696  if (!(trimNewVertices.find(v) != trimNewVertices.end())) {
1697  auto &m = edgesToTrim.at(e);
1698  verticesOnTrimEdges[v] = m;
1699  verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
1700  verticesOnTrimEdges[v].lEngth = 0;
1701  verticesOnTrimEdges[v].dIst = 0;
1702  trimNewVertices.insert(v);
1703  }
1704  };
1705 
1706  // Iterate over all vertives close to surface front and check if those can be
1707  // moved
1708 
1709  // filer nodes which distance is in given tolerance
1710  auto hi = verts_map.get<0>().upper_bound(tol);
1711  verts_map.get<0>().erase(hi, verts_map.end());
1712 
1713  auto remove_verts = [&](Range nodes) {
1714  for (auto n : nodes) {
1715  auto r = verts_map.get<1>().equal_range(n);
1716  verts_map.get<1>().erase(r.first, r.second);
1717  }
1718  };
1719 
1720  auto remove_edges = [&](Range edges) {
1721  for (auto e : edges) {
1722  auto r = verts_map.get<2>().equal_range(e);
1723  verts_map.get<2>().erase(r.first, r.second);
1724  }
1725  };
1726 
1727  auto trim_verts = [&](const Range verts, const bool move) {
1728  VertMapMultiIndex verts_map_tmp;
1729  for (auto p = corners.pair_begin(); p != corners.pair_end(); ++p) {
1730  auto lo = verts_map.get<1>().lower_bound(p->first);
1731  auto hi = verts_map.get<1>().upper_bound(p->second);
1732  verts_map_tmp.insert(lo, hi);
1733  }
1734  if(move) {
1735  for (auto &m : verts_map_tmp.get<0>())
1736  add_trim_vert(m.v, m.e);
1737  } else {
1738  for (auto &m : verts_map_tmp.get<0>())
1739  add_no_move_trim(m.v, m.e);
1740  }
1741  };
1742 
1743  auto trim_edges = [&](const Range ents, const bool move) {
1744  VertMapMultiIndex verts_map_tmp;
1745  for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
1746  auto lo = verts_map.get<2>().lower_bound(p->first);
1747  auto hi = verts_map.get<2>().upper_bound(p->second);
1748  verts_map_tmp.insert(lo, hi);
1749  verts_map.get<2>().erase(lo, hi);
1750  }
1751  if (move) {
1752  for (auto &m : verts_map_tmp.get<0>())
1753  add_trim_vert(m.v, m.e);
1754  } else {
1755  for (auto &m : verts_map_tmp.get<0>())
1756  add_no_move_trim(m.v, m.e);
1757  }
1758  };
1759 
1760  auto intersect_self = [&](Range &a, const Range b) { a = intersect(a, b); };
1761 
1762  map<std::string, Range> range_maps;
1763  CHKERR skin.find_skin(0, cutNewSurfaces, false, range_maps["surface_skin"]);
1764  intersect_self(range_maps["surface_skin"], trimEdges);
1765  range_maps["fixed_edges_on_surface_skin"] =
1766  intersect(range_maps["surface_skin"], fix_edges);
1767  CHKERR moab.get_adjacencies(range_maps["fixed_edges_verts"], 1, false,
1768  range_maps["fixed_edges_verts_edges"],
1769  moab::Interface::UNION);
1770  intersect_self(range_maps["fixed_edges_verts_edges"], trimEdges);
1771  CHKERR moab.get_connectivity(
1772  range_maps["fixed_edges_verts_edges"],
1773  range_maps["fixed_edges_verts_edges_verts_on_the_skin"], false);
1774  intersect_self(range_maps["fixed_edges_verts_edges_verts_on_the_skin"],
1775  tets_skin_verts);
1776 
1777  // do not move nodes at the corners
1778  trim_verts(corners, false);
1779  remove_verts(corners);
1780  trim_edges(range_maps["fixed_edges_on_surface_skin"], true);
1781  remove_verts(range_maps["fixed_edges_on_surface_skin_verts"]);
1782  trim_verts(range_maps["fixed_edges_verts_edges_verts_on_the_skin"], false);
1783  remove_verts(range_maps["fixed_edges_verts_edges_verts_on_the_skin"]);
1784  trim_edges(range_maps["surface_skin"], true);
1785  trim_verts(tets_skin_verts, false);
1786  remove_verts(tets_skin_verts);
1787 
1788  for (auto &m : verts_map.get<0>())
1789  add_trim_vert(m.v, m.e);
1790 
1791  for (auto m : verticesOnTrimEdges) {
1792  EntityHandle v = m.first;
1793  Range adj;
1794  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj);
1795  for (auto e : adj) {
1796  edgesToTrim.erase(e);
1797  trimEdges.erase(e);
1798  }
1799  }
1800 
1801  if (debug)
1802  if (!trimNewVertices.empty())
1803  CHKERR SaveData(moab)("trim_close_vertices.vtk", trimNewVertices);
1804 
1805  if (debug)
1806  if (!trimEdges.empty())
1807  CHKERR SaveData(moab)("trim_edges.vtk", trimEdges);
1808 
1810 }
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:341
map< EntityHandle, TreeData > verticesOnTrimEdges
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:91
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Types.hpp:93
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:595
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:88
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ getCutEdges()

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

Definition at line 460 of file CutMeshInterface.hpp.

460 { return cutEdges; }

◆ getCutFrontVolumes()

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

Definition at line 478 of file CutMeshInterface.hpp.

◆ getCutSurfaceVolumes()

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

Definition at line 477 of file CutMeshInterface.hpp.

◆ getCutVolumes()

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

Definition at line 461 of file CutMeshInterface.hpp.

461 { return cutVolumes; }

◆ getFront()

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

Definition at line 458 of file CutMeshInterface.hpp.

458 { return fRont; }

◆ getMergedSurfaces()

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

Definition at line 475 of file CutMeshInterface.hpp.

◆ getMergedVolumes()

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

Definition at line 474 of file CutMeshInterface.hpp.

474 { return mergedVolumes; }

◆ getNewCutSurfaces()

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

Definition at line 463 of file CutMeshInterface.hpp.

◆ getNewCutVertices()

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

Definition at line 464 of file CutMeshInterface.hpp.

◆ getNewCutVolumes()

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

Definition at line 462 of file CutMeshInterface.hpp.

462 { return cutNewVolumes; }

◆ getNewTrimSurfaces()

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

Definition at line 471 of file CutMeshInterface.hpp.

◆ getNewTrimVertices()

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

Definition at line 472 of file CutMeshInterface.hpp.

◆ getNewTrimVolumes()

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

Definition at line 470 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 
57  CHKERR PetscOptionsInt("-linesearch_steps",
58  "number of bisection steps which line search do to "
59  "find optimal merged nodes position",
60  "", lineSearchSteps, &lineSearchSteps, PETSC_NULL);
61 
62  CHKERR PetscOptionsInt("-max_merging_cycles",
63  "number of maximal merging cycles", "",
65 
66  CHKERR PetscOptionsScalar("-project_entities_quality_trashold",
67  "project entities quality trashold", "",
69  &projectEntitiesQualityTrashold, PETSC_NULL);
70 
71  ierr = PetscOptionsEnd();
72  CHKERRG(ierr);
74  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:543
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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 457 of file CutMeshInterface.hpp.

457 { return sUrface; }

◆ getTreeSurfPtr()

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

Definition at line 486 of file CutMeshInterface.hpp.

486  {
487  return treeSurfPtr;
488  }
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr

◆ getTrimEdges()

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

Definition at line 469 of file CutMeshInterface.hpp.

469 { return trimEdges; }

◆ getVolume()

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

Definition at line 456 of file CutMeshInterface.hpp.

456 { 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 459 of file CutMeshInterface.cpp.

459  {
460  CoreInterface &m_field = cOre;
461  moab::Interface &moab = m_field.get_moab();
463  Skinner skin(&moab);
464  Range tets_skin;
465  CHKERR skin.find_skin(0, vOlume, false, tets_skin);
466  Range tets_skin_edges;
467  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
468  moab::Interface::UNION);
469  Range surface_skin;
470  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
471  fRont = subtract(surface_skin, tets_skin_edges);
472  if (debug)
473  CHKERR SaveData(moab)("front.vtk", fRont);
475 }
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:476
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

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

2044  {
2045  CoreInterface &m_field = cOre;
2046  moab::Interface &moab = m_field.get_moab();
2048 
2049  /**
2050  * \brief Merge nodes
2051  */
2052  struct MergeNodes {
2053  CoreInterface &mField;
2054  const bool onlyIfImproveQuality;
2055  Tag tH;
2056  bool updateMehsets;
2057 
2058  MergeNodes(CoreInterface &m_field, const bool only_if_improve_quality,
2059  Tag th, bool update_mehsets)
2060  : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
2061  tH(th), updateMehsets(update_mehsets) {
2062  mField.getInterface(nodeMergerPtr);
2063  }
2064  NodeMergerInterface *nodeMergerPtr;
2065  MoFEMErrorCode mergeNodes(int line_search, EntityHandle father,
2066  EntityHandle mother, Range &proc_tets,
2067  bool add_child = true) {
2068  moab::Interface &moab = mField.get_moab();
2070  const EntityHandle conn[] = {father, mother};
2071  Range vert_tets;
2072  CHKERR moab.get_adjacencies(conn, 2, 3, false, vert_tets,
2073  moab::Interface::UNION);
2074  vert_tets = intersect(vert_tets, proc_tets);
2075  Range out_tets;
2076  CHKERR nodeMergerPtr->mergeNodes(father, mother, out_tets, &vert_tets,
2077  onlyIfImproveQuality, 0, line_search,
2078  tH);
2079 
2080  if (add_child && nodeMergerPtr->getSuccessMerge()) {
2081 
2082  Range::iterator lo, hi = proc_tets.begin();
2083  for (auto pt = vert_tets.pair_begin(); pt != vert_tets.pair_end();
2084  ++pt) {
2085  lo = proc_tets.lower_bound(hi, proc_tets.end(), pt->first);
2086  if (lo != proc_tets.end()) {
2087  hi = proc_tets.upper_bound(lo, proc_tets.end(), pt->second);
2088  proc_tets.erase(lo, hi);
2089  } else
2090  break;
2091  }
2092  proc_tets.merge(out_tets);
2093 
2094  auto &parent_child_map = nodeMergerPtr->getParentChildMap();
2095 
2096  struct ChangeChild {
2097  EntityHandle child;
2098  ChangeChild(const EntityHandle child) : child(child) {}
2099  void operator()(NodeMergerInterface::ParentChild &p) {
2100  p.cHild = child;
2101  }
2102  };
2103 
2104  std::vector<decltype(parentsChildMap.get<0>().begin())> it_vec;
2105  it_vec.reserve(parentsChildMap.size());
2106 
2107  for (auto &p : parent_child_map) {
2108 
2109  it_vec.clear();
2110  for (auto it = parentsChildMap.get<0>().equal_range(p.pArent);
2111  it.first != it.second; ++it.first)
2112  it_vec.emplace_back(it.first);
2113 
2114  for (auto it = parentsChildMap.get<1>().equal_range(p.pArent);
2115  it.first != it.second; ++it.first)
2116  it_vec.emplace_back(parentsChildMap.project<0>(it.first));
2117 
2118  for (auto &it : it_vec)
2119  parentsChildMap.modify(it, ChangeChild(p.cHild));
2120  }
2121 
2122  parentsChildMap.insert(parent_child_map.begin(),
2123  parent_child_map.end());
2124  }
2126  }
2127 
2128  MoFEMErrorCode updateRangeByChilds(Range &new_surf, Range &edges_to_merge,
2129  Range &not_merged_edges,
2130  bool add_child) {
2131  moab::Interface &moab = mField.get_moab();
2133  if (add_child) {
2134 
2135  std::vector<EntityHandle> parents_ents_vec(parentsChildMap.size());
2136  for (auto &it : parentsChildMap)
2137  parents_ents_vec.emplace_back(it.pArent);
2138  Range parent_ents;
2139  parent_ents.insert_list(parents_ents_vec.begin(),
2140  parents_ents_vec.end());
2141 
2142  Range surf_parent_ents = intersect(new_surf, parent_ents);
2143  new_surf = subtract(new_surf, surf_parent_ents);
2144  Range child_surf_ents;
2145  CHKERR updateRangeByChilds(parentsChildMap, surf_parent_ents,
2146  child_surf_ents);
2147  new_surf.merge(child_surf_ents);
2148 
2149  Range edges_to_merge_parent_ents =
2150  intersect(edges_to_merge, parent_ents);
2151  edges_to_merge = subtract(edges_to_merge, edges_to_merge_parent_ents);
2152  Range merged_child_edge_ents;
2153  CHKERR updateRangeByChilds(parentsChildMap, edges_to_merge_parent_ents,
2154  merged_child_edge_ents);
2155 
2156  Range not_merged_edges_child_ents =
2157  intersect(not_merged_edges, parent_ents);
2158  not_merged_edges =
2159  subtract(not_merged_edges, not_merged_edges_child_ents);
2160  Range not_merged_child_edge_ents;
2161  CHKERR updateRangeByChilds(parentsChildMap, not_merged_edges_child_ents,
2162  not_merged_child_edge_ents);
2163 
2164  merged_child_edge_ents =
2165  subtract(merged_child_edge_ents, not_merged_child_edge_ents);
2166  edges_to_merge.merge(merged_child_edge_ents);
2167  not_merged_edges.merge(not_merged_child_edge_ents);
2168 
2169  if (updateMehsets) {
2171  (*mField.getInterface<MeshsetsManager>()), cubit_it)) {
2172  EntityHandle cubit_meshset = cubit_it->meshset;
2173  Range meshset_ents;
2174  CHKERR moab.get_entities_by_handle(cubit_meshset, meshset_ents,
2175  true);
2176  Range child_ents;
2177  CHKERR updateRangeByChilds(parentsChildMap, meshset_ents,
2178  child_ents);
2179  CHKERR moab.add_entities(cubit_meshset, child_ents);
2180  }
2181  }
2182  }
2183 
2185  };
2186 
2187  private:
2188  NodeMergerInterface::ParentChildMap parentsChildMap;
2189  std::vector<EntityHandle> childsVec;
2190 
2191  inline MoFEMErrorCode updateRangeByChilds(
2192  const NodeMergerInterface::ParentChildMap &parent_child_map,
2193  const Range &parents, Range &childs) {
2195  childsVec.clear();
2196  childsVec.reserve(parents.size());
2197  for (auto pit = parents.pair_begin(); pit != parents.pair_end(); pit++) {
2198  auto it = parent_child_map.lower_bound(pit->first);
2199  if (it != parent_child_map.end()) {
2200  for (auto hi_it = parent_child_map.upper_bound(pit->second);
2201  it != hi_it; ++it)
2202  childsVec.emplace_back(it->cHild);
2203  }
2204  }
2205  childs.insert_list(childsVec.begin(), childsVec.end());
2207  }
2208  };
2209 
2210  /**
2211  * \brief Calculate edge length
2212  */
2213  struct LengthMap {
2214  Tag tH;
2215  CoreInterface &mField;
2216  moab::Interface &moab;
2217  const double maxLength;
2218  LengthMap(CoreInterface &m_field, Tag th, double max_length)
2219  : tH(th), mField(m_field), moab(m_field.get_moab()),
2220  maxLength(max_length) {
2221  gammaL = 1.;
2222  gammaQ = 3.;
2223  acceptedThrasholdMergeQuality = 0.5;
2224  }
2225 
2226  double
2227  gammaL; ///< Controls importance of length when ranking edges for merge
2228  double
2229  gammaQ; ///< Controls importance of quality when ranking edges for merge
2230  double acceptedThrasholdMergeQuality; ///< Do not merge quality if quality
2231  ///< above accepted thrashold
2232 
2233  MoFEMErrorCode operator()(const Range &tets, const Range &edges,
2234  LengthMapData_multi_index &length_map,
2235  double &ave) const {
2236  int num_nodes;
2237  const EntityHandle *conn;
2238  std::array<double, 12> coords;
2240  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
2241  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
2242  VectorDouble3 delta(3);
2243 
2244  struct NodeQuality {
2245  EntityHandle node;
2246  double quality;
2247  NodeQuality(const EntityHandle node) : node(node), quality(1) {}
2248  };
2249 
2250  typedef multi_index_container<
2251  NodeQuality, indexed_by<
2252 
2253  sequenced<>,
2254 
2255  hashed_non_unique<tag<Ent_mi_tag>,
2256  member<NodeQuality, EntityHandle,
2257  &NodeQuality::node>>
2258 
2259  >>
2260  NodeQuality_sequence;
2261 
2262  NodeQuality_sequence node_quality_sequence;
2263 
2264  Range edges_nodes;
2265  CHKERR moab.get_connectivity(tets, edges_nodes, false);
2266  Range edges_tets;
2267  CHKERR moab.get_adjacencies(edges, 3, false, edges_tets,
2268  moab::Interface::UNION);
2269  edges_tets = intersect(edges_tets, tets);
2270 
2271  for (auto node : edges_nodes)
2272  node_quality_sequence.emplace_back(node);
2273 
2274  for (auto tet : edges_tets) {
2275 
2276  CHKERR moab.get_connectivity(tet, conn, num_nodes, true);
2277  if (tH)
2278  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2279  else
2280  CHKERR moab.get_coords(conn, num_nodes, coords.data());
2281 
2282  const double q = Tools::volumeLengthQuality(coords.data());
2283 
2284  for (auto n : {0, 1, 2, 3}) {
2285  auto it = node_quality_sequence.get<1>().find(conn[n]);
2286  if (it != node_quality_sequence.get<1>().end())
2287  const_cast<double &>(it->quality) = std::min(q, it->quality);
2288  }
2289  }
2290 
2291  for (auto edge : edges) {
2292 
2293  CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
2294 
2295  if (tH)
2296  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2297  else
2298  CHKERR moab.get_coords(conn, num_nodes, coords.data());
2299 
2300  double q = 1;
2301  for (auto n : {0, 1}) {
2302  auto it = node_quality_sequence.get<1>().find(conn[n]);
2303  if (it != node_quality_sequence.get<1>().end())
2304  q = std::min(q, it->quality);
2305  }
2306 
2307  if (q < acceptedThrasholdMergeQuality) {
2308  noalias(delta) = (s0 - s1) / maxLength;
2309  double dot = inner_prod(delta, delta);
2310  double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
2311  length_map.insert(LengthMapData(val, q, edge));
2312  }
2313  }
2314 
2315  ave = 0;
2316  for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
2317  length_map.get<0>().begin();
2318  mit != length_map.get<0>().end(); mit++) {
2319  ave += mit->qUality;
2320  }
2321  ave /= length_map.size();
2323  }
2324  };
2325 
2326  /**
2327  * \brief Topological relations
2328  */
2329  struct Toplogy {
2330 
2331  CoreInterface &mField;
2332  Tag tH;
2333  const double tOL;
2334  Toplogy(CoreInterface &m_field, Tag th, const double tol)
2335  : mField(m_field), tH(th), tOL(tol) {}
2336 
2337  enum TYPE {
2338  FREE = 0,
2339  SKIN = 1 << 0,
2340  SURFACE = 1 << 1,
2341  SURFACE_SKIN = 1 << 2,
2342  FRONT_ENDS = 1 << 3,
2343  FIX_EDGES = 1 << 4,
2344  FIX_CORNERS = 1 << 5
2345  };
2346 
2347  typedef map<int, Range> SetsMap;
2348 
2349  MoFEMErrorCode classifyVerts(const Range &surface, const Range &tets,
2350  const Range &fixed_edges,
2351  const Range &corner_nodes,
2352  SetsMap &sets_map) const {
2353  moab::Interface &moab(mField.get_moab());
2354  Skinner skin(&moab);
2356 
2357  sets_map[FIX_CORNERS].merge(corner_nodes);
2358  Range fixed_verts;
2359  CHKERR moab.get_connectivity(fixed_edges, fixed_verts, true);
2360  sets_map[FIX_EDGES].swap(fixed_verts);
2361 
2362  Range tets_skin;
2363  CHKERR skin.find_skin(0, tets, false, tets_skin);
2364  Range tets_skin_edges;
2365  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2366  moab::Interface::UNION);
2367 
2368  // surface skin
2369  Range surface_skin;
2370  CHKERR skin.find_skin(0, surface, false, surface_skin);
2371  Range front_in_the_body;
2372  front_in_the_body = subtract(surface_skin, tets_skin_edges);
2373  Range front_ends;
2374  CHKERR skin.find_skin(0, front_in_the_body, false, front_ends);
2375  sets_map[FRONT_ENDS].swap(front_ends);
2376 
2377  Range surface_skin_verts;
2378  CHKERR moab.get_connectivity(surface_skin, surface_skin_verts, true);
2379  sets_map[SURFACE_SKIN].swap(surface_skin_verts);
2380 
2381  // surface
2382  Range surface_verts;
2383  CHKERR moab.get_connectivity(surface, surface_verts, true);
2384  sets_map[SURFACE].swap(surface_verts);
2385 
2386  // skin
2387  Range tets_skin_verts;
2388  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
2389  sets_map[SKIN].swap(tets_skin_verts);
2390 
2391  Range tets_verts;
2392  CHKERR moab.get_connectivity(tets, tets_verts, true);
2393  sets_map[FREE].swap(tets_verts);
2394 
2396  }
2397 
2398  MoFEMErrorCode getProcTets(const Range &tets, const Range &edges_to_merge,
2399  Range &proc_tets) const {
2400  moab::Interface &moab(mField.get_moab());
2402  Range edges_to_merge_verts;
2403  CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts, true);
2404  Range edges_to_merge_verts_tets;
2405  CHKERR moab.get_adjacencies(edges_to_merge_verts, 3, false,
2406  edges_to_merge_verts_tets,
2407  moab::Interface::UNION);
2408  edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
2409  proc_tets.swap(edges_to_merge_verts_tets);
2411  }
2412 
2413  MoFEMErrorCode removeBadEdges(const Range &surface, const Range &tets,
2414  const Range &fixed_edges,
2415  const Range &corner_nodes,
2416  Range &edges_to_merge,
2417  Range &not_merged_edges) {
2418  moab::Interface &moab(mField.get_moab());
2420 
2421  // find skin
2422  Skinner skin(&moab);
2423  Range tets_skin;
2424  CHKERR skin.find_skin(0, tets, false, tets_skin);
2425  Range surface_skin;
2426  CHKERR skin.find_skin(0, surface, false, surface_skin);
2427 
2428  // end nodes
2429  Range tets_skin_edges;
2430  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2431  moab::Interface::UNION);
2432  Range surface_front;
2433  surface_front = subtract(surface_skin, tets_skin_edges);
2434  Range surface_front_nodes;
2435  CHKERR moab.get_connectivity(surface_front, surface_front_nodes, true);
2436  Range ends_nodes;
2437  CHKERR skin.find_skin(0, surface_front, false, ends_nodes);
2438 
2439  // remove bad merges
2440 
2441  // get surface and body skin verts
2442  Range surface_edges;
2443  CHKERR moab.get_adjacencies(surface, 1, false, surface_edges,
2444  moab::Interface::UNION);
2445  // get nodes on the surface
2446  Range surface_edges_verts;
2447  CHKERR moab.get_connectivity(surface_edges, surface_edges_verts, true);
2448  // get vertices on the body skin
2449  Range tets_skin_edges_verts;
2450  CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts,
2451  true);
2452 
2453  Range edges_to_remove;
2454 
2455  // remove edges self connected to body skin
2456  {
2457  Range ents_nodes_and_edges;
2458  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2459  ents_nodes_and_edges.merge(tets_skin_edges);
2460  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2461  0, false);
2462  }
2463  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2464  not_merged_edges.merge(edges_to_remove);
2465 
2466  // remove edges self connected to surface
2467  {
2468  Range ents_nodes_and_edges;
2469  ents_nodes_and_edges.merge(surface_edges_verts);
2470  ents_nodes_and_edges.merge(surface_edges);
2471  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2472  ents_nodes_and_edges.merge(tets_skin_edges);
2473  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2474  0, false);
2475  }
2476  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2477  not_merged_edges.merge(edges_to_remove);
2478 
2479  // remove edges adjacent corner_nodes execpt those on fixed edges
2480  Range fixed_edges_nodes;
2481  CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes, true);
2482  {
2483  Range ents_nodes_and_edges;
2484  ents_nodes_and_edges.merge(fixed_edges_nodes);
2485  ents_nodes_and_edges.merge(ends_nodes);
2486  ents_nodes_and_edges.merge(corner_nodes);
2487  ents_nodes_and_edges.merge(fixed_edges);
2488  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2489  0, false);
2490  }
2491  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2492  not_merged_edges.merge(edges_to_remove);
2493 
2494  // remove edges self connected to surface
2495  CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove, 0, false);
2496  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2497  not_merged_edges.merge(edges_to_remove);
2498 
2499  // remove edges self contented on surface skin
2500  {
2501  Range ents_nodes_and_edges;
2502  ents_nodes_and_edges.merge(surface_skin);
2503  ents_nodes_and_edges.merge(fixed_edges_nodes);
2504  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2505  0, false);
2506  }
2507  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2508  not_merged_edges.merge(edges_to_remove);
2509 
2510  // remove edges connecting crack front and fixed edges, except those short
2511  {
2512  Range ents_nodes_and_edges;
2513  ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
2514  ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2515  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2516  0, false);
2517  }
2518  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2519  not_merged_edges.merge(edges_to_remove);
2520 
2521  // remove crack front nodes connected to the surface, except those short
2522  {
2523  Range ents_nodes_and_edges;
2524  ents_nodes_and_edges.merge(surface_front_nodes);
2525  ents_nodes_and_edges.merge(surface_front);
2526  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2527  ents_nodes_and_edges.merge(tets_skin_edges);
2528  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2529  tOL, false);
2530  }
2531  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2532  not_merged_edges.merge(edges_to_remove);
2533 
2535  }
2536 
2537  private:
2538  MoFEMErrorCode removeSelfConectingEdges(const Range &ents,
2539  Range &edges_to_remove,
2540  const bool length,
2541  bool debug) const {
2542  moab::Interface &moab(mField.get_moab());
2544  // get nodes
2545  Range ents_nodes = ents.subset_by_type(MBVERTEX);
2546  if (ents_nodes.empty()) {
2547  CHKERR moab.get_connectivity(ents, ents_nodes, true);
2548  }
2549  // edges adj. to nodes
2550  Range ents_nodes_edges;
2551  CHKERR moab.get_adjacencies(ents_nodes, 1, false, ents_nodes_edges,
2552  moab::Interface::UNION);
2553  // nodes of adj. edges
2554  Range ents_nodes_edges_nodes;
2555  CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2556  true);
2557  // hanging nodes
2558  ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2559  Range ents_nodes_edges_nodes_edges;
2560  CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1, false,
2561  ents_nodes_edges_nodes_edges,
2562  moab::Interface::UNION);
2563  // remove edges adj. to hanging edges
2564  ents_nodes_edges =
2565  subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2566  ents_nodes_edges =
2567  subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2568  if (length > 0) {
2569  Range::iterator eit = ents_nodes_edges.begin();
2570  for (; eit != ents_nodes_edges.end();) {
2571 
2572  int num_nodes;
2573  const EntityHandle *conn;
2574  CHKERR moab.get_connectivity(*eit, conn, num_nodes, true);
2575  double coords[6];
2576  if (tH)
2577  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
2578  else
2579  CHKERR moab.get_coords(conn, num_nodes, coords);
2580 
2581  auto get_edge_length = [coords]() {
2583  &coords[0], &coords[1], &coords[2]);
2586  t_delta(i) = t_coords(i);
2587  ++t_coords;
2588  t_delta(i) -= t_coords(i);
2589  return sqrt(t_delta(i) * t_delta(i));
2590  };
2591 
2592  if (get_edge_length() < tOL) {
2593  eit = ents_nodes_edges.erase(eit);
2594  } else {
2595  eit++;
2596  }
2597  }
2598  }
2599  edges_to_remove.swap(ents_nodes_edges);
2600  if (debug) {
2601  CHKERR SaveData(moab)("edges_to_remove.vtk", edges_to_remove);
2602  }
2604  }
2605  };
2606 
2607  Range not_merged_edges;
2608  const double tol = 1e-3;
2609  CHKERR Toplogy(m_field, th, tol * aveLength)
2610  .removeBadEdges(surface, tets, fixed_edges, corner_nodes, edges_to_merge,
2611  not_merged_edges);
2612  Toplogy::SetsMap sets_map;
2613  CHKERR Toplogy(m_field, th, tol * aveLength)
2614  .classifyVerts(surface, tets, fixed_edges, corner_nodes, sets_map);
2615  if (debug) {
2616  for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2617  sit != sets_map.rend(); sit++) {
2618  std::string name = "classification_verts_" +
2619  boost::lexical_cast<std::string>(sit->first) + ".vtk";
2620  if (!sit->second.empty())
2621  CHKERR SaveData(moab)(name, sit->second);
2622  }
2623  }
2624  Range proc_tets;
2625  CHKERR Toplogy(m_field, th, tol * aveLength)
2626  .getProcTets(tets, edges_to_merge, proc_tets);
2627  out_tets = subtract(tets, proc_tets);
2628 
2629  if (bit_ptr) {
2630  Range all_out_ents = out_tets;
2631  for (int dd = 2; dd >= 0; dd--) {
2632  CHKERR moab.get_adjacencies(out_tets, dd, false, all_out_ents,
2633  moab::Interface::UNION);
2634  }
2635  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(all_out_ents,
2636  *bit_ptr);
2637  }
2638 
2639  int nb_nodes_merged = 0;
2640  LengthMapData_multi_index length_map;
2641  new_surf = surface;
2642 
2643  auto save_merge_step = [&](const int pp, const Range collapsed_edges) {
2645  Range adj_faces;
2646  CHKERR moab.get_adjacencies(proc_tets, 2, false, adj_faces,
2647  moab::Interface::UNION);
2648  std::string name;
2649  name = "node_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2650  CHKERR SaveData(moab)(
2651  name, unite(intersect(new_surf, adj_faces), collapsed_edges));
2652  name =
2653  "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2654  CHKERR SaveData(moab)(
2655  name, unite(intersect(new_surf, adj_faces), edges_to_merge));
2657  };
2658 
2659  if (debug)
2660  CHKERR save_merge_step(0, Range());
2661 
2662  double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2663  for (int pp = 0; pp != nbMaxMergingCycles; ++pp) {
2664 
2665  int nb_nodes_merged_p = nb_nodes_merged;
2666  length_map.clear();
2667  min_pp = min_p;
2668  min_p = min;
2669  CHKERR LengthMap(m_field, th, aveLength)(proc_tets, edges_to_merge,
2670  length_map, ave);
2671  min = length_map.get<2>().begin()->qUality;
2672  if (pp == 0) {
2673  ave0 = ave;
2674  }
2675 
2676  int nn = 0;
2677  Range collapsed_edges;
2678  MergeNodes merge_nodes(m_field, true, th, update_meshsets);
2679 
2680  for (auto mit = length_map.get<0>().begin();
2681  mit != length_map.get<0>().end(); mit++, nn++) {
2682 
2683  if (!mit->skip) {
2684 
2685  auto get_father_and_mother =
2686  [&](EntityHandle &father, EntityHandle &mother, int &line_search) {
2688  int num_nodes;
2689  const EntityHandle *conn;
2690  CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes, true);
2691  std::array<int, 2> conn_type = {0, 0};
2692  for (int nn = 0; nn != 2; nn++) {
2693  for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2694  sit != sets_map.rend(); sit++) {
2695  if (sit->second.find(conn[nn]) != sit->second.end()) {
2696  conn_type[nn] |= sit->first;
2697  }
2698  }
2699  }
2700  int type_father, type_mother;
2701  if (conn_type[0] > conn_type[1]) {
2702  father = conn[0];
2703  mother = conn[1];
2704  type_father = conn_type[0];
2705  type_mother = conn_type[1];
2706  } else {
2707  father = conn[1];
2708  mother = conn[0];
2709  type_father = conn_type[1];
2710  type_mother = conn_type[0];
2711  }
2712  if (type_father == type_mother) {
2713  line_search = lineSearchSteps;
2714  }
2716  };
2717 
2718  int line_search = 0;
2719  EntityHandle father, mother;
2720  CHKERR get_father_and_mother(father, mother, line_search);
2721  CHKERR merge_nodes.mergeNodes(line_search, father, mother, proc_tets);
2722  if (m_field.getInterface<NodeMergerInterface>()->getSuccessMerge()) {
2723  Range adj_mother_tets;
2724  CHKERR moab.get_adjacencies(&mother, 1, 3, false, adj_mother_tets);
2725  Range adj_mother_tets_nodes;
2726  CHKERR moab.get_connectivity(adj_mother_tets, adj_mother_tets_nodes,
2727  true);
2728  Range adj_edges;
2729  CHKERR moab.get_adjacencies(adj_mother_tets_nodes, 1, false,
2730  adj_edges, moab::Interface::UNION);
2731  CHKERR moab.get_adjacencies(&father, 1, 1, false, adj_edges,
2732  moab::Interface::UNION);
2733  Range proc_edges;
2734  CHKERR moab.get_adjacencies(proc_tets, 1, true, proc_edges);
2735  adj_edges = intersect(proc_edges, adj_edges);
2736  for (auto ait : adj_edges) {
2737  auto miit = length_map.get<1>().find(ait);
2738  if (miit != length_map.get<1>().end())
2739  (const_cast<LengthMapData &>(*miit)).skip = true;
2740  }
2741  nb_nodes_merged++;
2742  collapsed_edges.insert(mit->eDge);
2743  }
2744 
2745  if (nn > static_cast<int>(length_map.size() / fraction_level))
2746  break;
2747  if (mit->qUality > ave)
2748  break;
2749  }
2750  }
2751 
2752  CHKERR merge_nodes.updateRangeByChilds(new_surf, edges_to_merge,
2753  not_merged_edges, true);
2754 
2755  PetscPrintf(m_field.get_comm(),
2756  "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e\n", pp,
2757  nb_nodes_merged, ave, min);
2758 
2759  if (debug)
2760  CHKERR save_merge_step(pp + 1, collapsed_edges);
2761 
2762  if (nb_nodes_merged == nb_nodes_merged_p)
2763  break;
2764  if (min > 1e-2 && min == min_pp)
2765  break;
2766  if (min > ave0)
2767  break;
2768 
2769  Range adj_edges;
2770  CHKERR moab.get_adjacencies(proc_tets, 1, false, adj_edges,
2771  moab::Interface::UNION);
2772  edges_to_merge = intersect(edges_to_merge, adj_edges);
2773  CHKERR Toplogy(m_field, th, tol * aveLength)
2774  .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
2775  edges_to_merge, not_merged_edges);
2776  }
2777 
2778  if (bit_ptr)
2779  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(proc_tets,
2780  *bit_ptr);
2781 
2782  out_tets.merge(proc_tets);
2783  Range adj_faces;
2784  CHKERR moab.get_adjacencies(out_tets, 2, false, adj_faces,
2785  moab::Interface::UNION);
2786  new_surf = intersect(new_surf, adj_faces);
2787 
2789 }
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:500
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
double aveLength
Average edge length.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
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:104
#define CHKERR
Inline error check.
Definition: definitions.h:595
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:88
#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:406

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

2795  {
2796  CoreInterface &m_field = cOre;
2798  Range tets_level;
2799  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2800  trim_bit, BitRefLevel().set(), MBTET, tets_level);
2801 
2802  Range edges_to_merge;
2803  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
2804  trim_bit, cut_bit | trim_bit, edges_to_merge);
2805 
2806  // get all entities not in database
2807  Range all_ents_not_in_database_before;
2808  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2809  all_ents_not_in_database_before);
2810 
2811  edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
2812  if(debug)
2813  CHKERR SaveData(m_field.get_moab())("edges_to_merge.vtk", edges_to_merge);
2814 
2815  Range out_new_tets, out_new_surf;
2816  CHKERR mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
2817  corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
2818  th, update_meshsets, &bit, debug);
2819 
2820  // get all entities not in database after merge
2821  Range all_ents_not_in_database_after;
2822  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2823  all_ents_not_in_database_after);
2824 
2825  // delete hanging entities
2826  all_ents_not_in_database_after =
2827  subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
2828  Range meshsets;
2829  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
2830  false);
2831  for (auto m : meshsets)
2832  CHKERR m_field.get_moab().remove_entities(m,
2833  all_ents_not_in_database_after);
2834 
2835  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
2836 
2837  mergedVolumes.swap(out_new_tets);
2838  mergedSurfaces.swap(out_new_surf);
2840 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:595
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:406

◆ 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:500
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507

◆ 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:500
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507

◆ moveMidNodesOnCutEdges()

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

projecting of mid edge nodes on new mesh on surface

Returns
error code

Definition at line 1387 of file CutMeshInterface.cpp.

1387  {
1389 
1390  CoreInterface &m_field = cOre;
1391  moab::Interface &moab = m_field.get_moab();
1393 
1394  // Range out_side_vertices;
1395  for (auto m : verticesOnCutEdges) {
1396  double dist = m.second.dIst;
1397  VectorDouble3 new_coors = m.second.rayPoint + dist * m.second.unitRayDir;
1398  if (th)
1399  CHKERR moab.tag_set_data(th, &m.first, 1, &new_coors[0]);
1400  else
1401  CHKERR moab.set_coords(&m.first, 1, &new_coors[0]);
1402  }
1403 
1405 }
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:500
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:88
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:406

◆ moveMidNodesOnTrimmedEdges()

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

move trimmed edges mid nodes

Returns
error code

Definition at line 1407 of file CutMeshInterface.cpp.

1407  {
1408  CoreInterface &m_field = cOre;
1409  moab::Interface &moab = m_field.get_moab();
1411  for (auto &v : verticesOnTrimEdges) {
1412  double dist = v.second.dIst;
1413  VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1414  if (th)
1415  CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1416  else
1417  CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1418  }
1420 }
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:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:88
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ projectZeroDistanceEnts() [1/2]

MoFEMErrorCode MoFEM::CutMeshInterface::projectZeroDistanceEnts ( Range *  fixed_edges,
Range *  corner_nodes,
const double  low_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
low_tois tolerance how close entities has to be
verbverbosity level
debugtrue for debuging purposes

Definition at line 962 of file CutMeshInterface.cpp.

966  {
967  CoreInterface &m_field = cOre;
968  moab::Interface &moab = m_field.get_moab();
970 
971  // Get entities on body skin
972  Skinner skin(&moab);
973  Range tets_skin;
974  CHKERR skin.find_skin(0, vOlume, false, tets_skin);
975  Range tets_skin_edges;
976  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
977  moab::Interface::UNION);
978  Range tets_skin_verts;
979  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
980 
981  // Get entities in volume
982  Range vol_faces, vol_edges, vol_nodes;
983  CHKERR moab.get_adjacencies(vOlume, 2, false, vol_faces,
984  moab::Interface::UNION);
985  CHKERR moab.get_adjacencies(vOlume, 1, false, vol_edges,
986  moab::Interface::UNION);
987  CHKERR moab.get_connectivity(vOlume, vol_nodes, true);
988 
989  // Get nodes on cut edges
990  Range cut_edge_verts;
991  CHKERR moab.get_connectivity(cutEdges, cut_edge_verts, true);
992 
993  // Get faces and edges
994  Range cut_edges_faces;
995  CHKERR moab.get_adjacencies(cut_edge_verts, 2, true, cut_edges_faces,
996  moab::Interface::UNION);
997  cut_edges_faces = intersect(cut_edges_faces, vol_faces);
998  Range cut_edges_faces_verts;
999  CHKERR moab.get_connectivity(cut_edges_faces, cut_edges_faces_verts, true);
1000  cut_edges_faces_verts = subtract(cut_edges_faces_verts, cut_edge_verts);
1001  Range to_remove_cut_edges_faces;
1002  CHKERR moab.get_adjacencies(cut_edges_faces_verts, 2, true,
1003  to_remove_cut_edges_faces,
1004  moab::Interface::UNION);
1005  // Those are faces which have vertices adjacent to cut edges vertices without
1006  // hanging nodes (nodes not adjacent to cutting edge)
1007  cut_edges_faces = subtract(cut_edges_faces, to_remove_cut_edges_faces);
1008  if (debug)
1009  CHKERR SaveData(moab)("cut_edges_faces.vtk", cut_edges_faces);
1010  cut_edges_faces.merge(cutEdges);
1011 
1012  Range fixed_edges_nodes;
1013  if (fixed_edges)
1014  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_nodes, true);
1015 
1016  Tag th_dist_normal;
1017  CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_VECTOR", th_dist_normal);
1018 
1019  // Create map of closes points to the surface
1020  enum TYPE { FREE = 0, SKIN, FIXED, CORNER, NOT_MOVE };
1021  map<EntityHandle, std::pair<std::pair<TYPE, EntityHandle>, TreeData>>
1022  min_dist_map;
1023  double ave_cut_edge_length = 0;
1024  for (auto e : cutEdges) {
1025 
1026  TYPE edge_type = FREE;
1027  if (tets_skin_edges.find(e) != tets_skin_edges.end())
1028  edge_type = SKIN;
1029  if (fixed_edges)
1030  if (fixed_edges->find(e) != fixed_edges->end())
1031  edge_type = FIXED;
1032 
1033  auto eit = edgesToCut.find(e);
1034  if (eit != edgesToCut.end()) {
1035 
1036  int num_nodes;
1037  const EntityHandle *conn;
1038  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1039  VectorDouble6 pos(6);
1040  CHKERR moab.get_coords(conn, num_nodes, &pos[0]);
1041  VectorDouble3 p[2];
1042  p[0] = getVectorAdaptor(&pos[0], 3);
1043  p[1] = getVectorAdaptor(&pos[3], 3);
1044  ave_cut_edge_length += norm_2(p[0] - p[1]);
1045 
1046  VectorDouble6 dist_normal(6);
1047  CHKERR moab.tag_get_data(th_dist_normal, conn, num_nodes,
1048  &dist_normal[0]);
1049  auto get_normal_dist = [](const double *normal) {
1050  FTensor::Tensor1<const double *, 3> t_n(normal, &normal[1], &normal[2]);
1052  return sqrt(t_n(i) * t_n(i));
1053  };
1054 
1055  auto &d = eit->second;
1056  VectorDouble3 new_pos = d.rayPoint + d.dIst * d.unitRayDir;
1057  for (int nn = 0; nn != 2; ++nn) {
1058 
1059  VectorDouble3 ray = new_pos - p[nn];
1060  const double dist = norm_2(ray);
1061  const double length = get_normal_dist(&dist_normal[3 * nn]);
1062 
1063  bool add_node = true;
1064  auto vit = min_dist_map.find(conn[nn]);
1065  if (vit != min_dist_map.end()) {
1066  if (vit->second.first.first > edge_type)
1067  add_node = false;
1068  else if (vit->second.first.first == edge_type) {
1069  if (vit->second.second.dIst < dist)
1070  add_node = false;
1071  }
1072  }
1073 
1074  if (add_node) {
1075  min_dist_map[conn[nn]].first.first = edge_type;
1076  min_dist_map[conn[nn]].first.second = e;
1077  auto &data = min_dist_map[conn[nn]].second;
1078  data.lEngth = length;
1079  data.rayPoint = p[nn];
1080  data.dIst = dist;
1081  if (dist > 0)
1082  data.unitRayDir = ray / dist;
1083  else {
1084  data.unitRayDir.resize(3);
1085  data.unitRayDir.clear();
1086  }
1087  }
1088  }
1089 
1090  } else
1091  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Edge not found");
1092  }
1093 
1094  ave_cut_edge_length /= cutEdges.size();
1095 
1096  auto get_quality_change =
1097  [this, &m_field,
1098  &moab](const Range &adj_tets,
1099  map<EntityHandle, TreeData> vertices_on_cut_edges) {
1100  vertices_on_cut_edges.insert(verticesOnCutEdges.begin(),
1101  verticesOnCutEdges.end());
1102  double q0 = 1;
1103  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0);
1104  double min_q = 1;
1105  for (auto t : adj_tets) {
1106  int num_nodes;
1107  const EntityHandle *conn;
1108  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1109  VectorDouble12 coords(12);
1110  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1111  for (int n = 0; n != 4; ++n) {
1112  bool ray_found = false;
1113  auto mit = vertices_on_cut_edges.find(conn[n]);
1114  if (mit != vertices_on_cut_edges.end()) {
1115  ray_found = true;
1116  }
1117  if (ray_found) {
1118  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1119  double dist = mit->second.dIst;
1120  noalias(n_coords) =
1121  mit->second.rayPoint + dist * mit->second.unitRayDir;
1122  }
1123  }
1124  const double q = Tools::volumeLengthQuality(&coords[0]);
1125  if (!std::isnormal(q)) {
1126  min_q = -2;
1127  break;
1128  }
1129  min_q = std::min(min_q, q);
1130  }
1131  return min_q / q0;
1132  };
1133 
1134  auto get_conn = [&moab](const EntityHandle &e, const EntityHandle *&conn,
1135  int &num_nodes) {
1137  EntityType type = moab.type_from_handle(e);
1138  if (type == MBVERTEX) {
1139  conn = &e;
1140  num_nodes = 1;
1141  } else {
1142  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1143  }
1145  };
1146 
1147  auto project_node = [&](const EntityHandle v,
1148  map<EntityHandle, TreeData> &vertices_on_cut_edges) {
1150  vertices_on_cut_edges[v].dIst = min_dist_map[v].second.dIst;
1151  vertices_on_cut_edges[v].lEngth = min_dist_map[v].second.lEngth;
1152  vertices_on_cut_edges[v].unitRayDir = min_dist_map[v].second.unitRayDir;
1153  vertices_on_cut_edges[v].rayPoint = min_dist_map[v].second.rayPoint;
1155  };
1156 
1157  auto remove_surface_tets = [&](Range &zero_dist_tets,
1158  Range zero_distance_ents,
1159  Range zero_distance_verts) {
1161  Range zero_dist_all_verts;
1162  CHKERR moab.get_connectivity(zero_distance_ents, zero_dist_all_verts, true);
1163  zero_dist_all_verts.merge(zero_distance_verts);
1164  CHKERR moab.get_adjacencies(zero_dist_all_verts, 3, false, zero_dist_tets,
1165  moab::Interface::UNION);
1166  Range zero_tets_verts;
1167  CHKERR moab.get_connectivity(zero_dist_tets, zero_tets_verts, true);
1168  zero_tets_verts = subtract(zero_tets_verts, zero_dist_all_verts);
1169  Range free_tets;
1170  CHKERR moab.get_adjacencies(zero_tets_verts, 3, false, free_tets,
1171  moab::Interface::UNION);
1172  zero_dist_tets = subtract(zero_dist_tets, free_tets);
1173 
1175  };
1176 
1177  for (int d = 2; d >= 0; --d) {
1178 
1179  Range ents;
1180  if (d > 0)
1181  ents = cut_edges_faces.subset_by_dimension(d);
1182  else
1183  ents = cut_edge_verts;
1184 
1185  // make list of entities
1186  multimap<double, EntityHandle> ents_to_check;
1187  for (auto f : ents) {
1188  int num_nodes;
1189  const EntityHandle *conn;
1190  CHKERR get_conn(f, conn, num_nodes);
1191  VectorDouble3 dist(3);
1192 
1193  bool move = true;
1194  for (int n = 0; n != num_nodes; ++n) {
1195  auto &d = min_dist_map[conn[n]];
1196  if (d.first.first == NOT_MOVE) {
1197  move = false;
1198  break;
1199  }
1200  dist[n] = d.second.lEngth;
1201  }
1202 
1203  if (move) {
1204  double max_dist = 0;
1205  for (int n = 0; n != num_nodes; ++n)
1206  max_dist = std::max(max_dist, fabs(dist[n]));
1207  if (max_dist < low_tol * ave_cut_edge_length)
1208  ents_to_check.insert(std::pair<double, EntityHandle>(max_dist, f));
1209  }
1210  }
1211 
1212  double ray_point[3], unit_ray_dir[3];
1213  VectorAdaptor vec_unit_ray_dir(
1214  3, ublas::shallow_array_adaptor<double>(3, unit_ray_dir));
1215  VectorAdaptor vec_ray_point(
1216  3, ublas::shallow_array_adaptor<double>(3, ray_point));
1217 
1218  for (auto m : ents_to_check) {
1219 
1220  EntityHandle f = m.second;
1221 
1222  int num_nodes;
1223  const EntityHandle *conn;
1224  CHKERR get_conn(f, conn, num_nodes);
1225  VectorDouble9 coords(9);
1226  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1227 
1228  Range adj_tets;
1229  CHKERR moab.get_adjacencies(conn, num_nodes, 3, false, adj_tets,
1230  moab::Interface::UNION);
1231  adj_tets = intersect(adj_tets, vOlume);
1232 
1233  map<EntityHandle, TreeData> vertices_on_cut_edges;
1234  for (int n = 0; n != num_nodes; ++n)
1235  CHKERR project_node(conn[n], vertices_on_cut_edges);
1236 
1237  const double q = get_quality_change(adj_tets, vertices_on_cut_edges);
1238 
1240  EntityHandle type = moab.type_from_handle(f);
1241 
1242  Range zero_dist_tets;
1243  if (type == MBVERTEX) {
1244  Range zero_distance_verts_test = zeroDistanceVerts;
1245  zero_distance_verts_test.insert(f);
1246  CHKERR remove_surface_tets(zero_dist_tets, zeroDistanceEnts,
1247  zero_distance_verts_test);
1248  } else {
1249  Range zero_distance_ents_test = zeroDistanceEnts;
1250  zero_distance_ents_test.insert(f);
1251  CHKERR remove_surface_tets(zero_dist_tets, zero_distance_ents_test,
1253  }
1254 
1255  if (zero_dist_tets.empty()) {
1256  verticesOnCutEdges.insert(vertices_on_cut_edges.begin(),
1257  vertices_on_cut_edges.end());
1258 
1259  if (type == MBVERTEX) {
1260  zeroDistanceVerts.insert(f);
1261  } else {
1262  zeroDistanceEnts.insert(f);
1263  }
1264  }
1265  }
1266  }
1267  }
1268 
1269  for (auto &v : verticesOnCutEdges) {
1270 
1271  TYPE node_type;
1272 
1273  if (corner_nodes->find(v.first) != corner_nodes->end())
1274  node_type = CORNER;
1275  else if (fixed_edges_nodes.find(v.first) != fixed_edges_nodes.end())
1276  node_type = FIXED;
1277  else if (tets_skin_verts.find(v.first) != tets_skin_verts.end())
1278  node_type = SKIN;
1279  else
1280  node_type = FREE;
1281 
1282  if (node_type > min_dist_map[v.first].first.first)
1283  v.second.unitRayDir.clear();
1284  }
1285 
1286  for (auto f : unite(zeroDistanceEnts, zeroDistanceVerts)) {
1287  int num_nodes;
1288  const EntityHandle *conn;
1289  CHKERR get_conn(f, conn, num_nodes);
1290  Range adj_edges;
1291  CHKERR moab.get_adjacencies(conn, num_nodes, 1, false, adj_edges,
1292  moab::Interface::UNION);
1293  for (auto e : adj_edges) {
1294  cutEdges.erase(e);
1295  edgesToCut.erase(e);
1296  }
1297  }
1298 
1299  if (debug)
1300  SaveData(m_field.get_moab())("projected_zero_distance_ents.vtk",
1302 
1304 }
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:500
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:91
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Types.hpp:93
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:104
#define CHKERR
Inline error check.
Definition: definitions.h:595
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:88
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:406
VectorBoundedArray< double, 9 > VectorDouble9
Definition: Types.hpp:92
map< EntityHandle, TreeData > edgesToCut

◆ projectZeroDistanceEnts() [2/2]

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

Definition at line 465 of file CutMeshInterface.hpp.

465  {
466  return zeroDistanceEnts;
467  }

◆ 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:500
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
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 715 of file CutMeshInterface.cpp.

717  {
718  CoreInterface &m_field = cOre;
719  moab::Interface &moab = m_field.get_moab();
720  MeshRefinement *refiner;
721  BitRefManager *bit_ref_manager;
723  CHKERR m_field.getInterface(refiner);
724  CHKERR m_field.getInterface(bit_ref_manager);
725 
726  auto add_bit = [&](const int bit) {
728  CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
729  verb);
730  Range adj_ents;
731  for (auto d : {2, 1, 0})
732  CHKERR moab.get_adjacencies(vOlume, d, true, adj_ents,
733  moab::Interface::UNION);
734  CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
735  verb);
737  };
738  CHKERR add_bit(init_bit_level);
739 
740  int very_last_bit = init_bit_level + surf_levels + front_levels + 2;
741 
742  auto update_range = [&](Range *r_ptr) {
744  if (r_ptr) {
745  Range childs;
746  CHKERR bit_ref_manager->updateRange(*r_ptr, childs);
747  r_ptr->merge(childs);
748  }
750  };
751 
752  auto refine = [&](const BitRefLevel bit, const Range tets) {
754  Range verts;
755  CHKERR moab.get_connectivity(tets, verts, true);
756  Range ref_edges;
757  CHKERR moab.get_adjacencies(verts, 1, true, ref_edges,
758  moab::Interface::UNION);
759 
760  CHKERR refiner->add_vertices_in_the_middel_of_edges(ref_edges, bit);
761  CHKERR refiner->refine_TET(vOlume, bit, false, verb);
762 
763  CHKERR update_range(fixed_edges);
764  CHKERR update_range(&vOlume);
765  CHKERR m_field.getInterface<MeshsetsManager>()
766  ->updateAllMeshsetsByEntitiesChildren(bit);
767 
768  Range bit_tets;
769  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
770  bit, BitRefLevel().set(), MBTET, bit_tets);
771  vOlume = intersect(vOlume, bit_tets);
772 
773  Range bit_edges;
774  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
775  bit, BitRefLevel().set(), MBEDGE, bit_edges);
776  if (fixed_edges)
777  *fixed_edges = intersect(*fixed_edges, bit_edges);
778 
780  };
781 
782  for (int ll = init_bit_level; ll != init_bit_level + surf_levels; ++ll) {
784  CHKERR refine(BitRefLevel().set(ll + 1),
786  }
787 
788  for (int ll = init_bit_level + surf_levels;
789  ll != init_bit_level + surf_levels + front_levels; ++ll) {
791  CHKERR refine(BitRefLevel().set(ll + 1), cutFrontVolumes);
792  }
793 
795 
796  if (debug)
797  CHKERR SaveData(m_field.get_moab())("refinedVolume.vtk", vOlume);
798 
800 }
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:476
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:595
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:406

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

1974  {
1975  CoreInterface &m_field = cOre;
1976  moab::Interface &moab = m_field.get_moab();
1977  PrismInterface *interface;
1979  CHKERR m_field.getInterface(interface);
1980  // Remove tris on surface front
1981  {
1982  Range front_tris;
1983  EntityHandle meshset;
1984  CHKERR moab.create_meshset(MESHSET_SET, meshset);
1985  CHKERR moab.add_entities(meshset, ents);
1986  CHKERR interface->findIfTringleHasThreeNodesOnInternalSurfaceSkin(
1987  meshset, split_bit, true, front_tris);
1988  CHKERR moab.delete_entities(&meshset, 1);
1989  ents = subtract(ents, front_tris);
1990  }
1991  // Remove entities on skin
1992  Range tets;
1993  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1994  split_bit, BitRefLevel().set(), MBTET, tets);
1995  // Remove entities on skin
1996  Skinner skin(&moab);
1997  Range tets_skin;
1998  CHKERR skin.find_skin(0, tets, false, tets_skin);
1999  ents = subtract(ents, tets_skin);
2000 
2002 }
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:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
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:406

◆ saveCutEdges()

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

Definition at line 2875 of file CutMeshInterface.cpp.

2875  {
2876  CoreInterface &m_field = cOre;
2877  moab::Interface &moab = m_field.get_moab();
2879  CHKERR SaveData(moab)(prefix + "out_vol.vtk", vOlume);
2880  CHKERR SaveData(moab)(prefix + "out_surface.vtk", sUrface);
2881  CHKERR SaveData(moab)(prefix + "out_cut_edges.vtk", cutEdges);
2882  CHKERR SaveData(moab)(prefix + "out_cut_volumes.vtk", cutVolumes);
2883  CHKERR SaveData(moab)(prefix + "out_cut_new_volumes.vtk", cutNewVolumes);
2884  CHKERR SaveData(moab)(prefix + "out_cut_new_surfaces.vtk", cutNewSurfaces);
2885  CHKERR SaveData(moab)(prefix + "out_cut_zero_distance_ents.vtk",
2888 }
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:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ saveTrimEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveTrimEdges ( )

Definition at line 2890 of file CutMeshInterface.cpp.

2890  {
2891  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
2893  CHKERR SaveData(moab)("out_trim_surface.vtk", sUrface);
2894  CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
2895  CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
2896  CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
2898 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

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

2859  {
2860  CoreInterface &m_field = cOre;
2861  moab::Interface &moab = m_field.get_moab();
2863  Range nodes;
2864  if (bit.none())
2865  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
2866  else
2867  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2868  bit, mask, MBVERTEX, nodes);
2869  std::vector<double> coords(3 * nodes.size());
2870  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
2871  CHKERR moab.set_coords(nodes, &coords[0]);
2873 }
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:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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:500
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507

◆ 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:500
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507

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

2842  {
2843  CoreInterface &m_field = cOre;
2844  moab::Interface &moab = m_field.get_moab();
2846  Range nodes;
2847  if (bit.none())
2848  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
2849  else
2850  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2851  bit, BitRefLevel().set(), MBVERTEX, nodes);
2852  std::vector<double> coords(3 * nodes.size());
2853  CHKERR moab.get_coords(nodes, &coords[0]);
2854  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
2856 }
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:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
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:406

◆ setTrimFixedEdges()

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

Definition at line 480 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:500
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507

◆ 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:476
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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  auto get_point = [i](auto &t_w, auto &t_delta, auto t) {
167  t = std::max(0., std::min(1., t));
168  t_p(i) = t_w(i) + t * t_delta(i);
169  return t_p;
170  };
171 
172  auto get_distance = [i](auto &t_p, auto &t_n) {
173  FTensor::Tensor1<double, 3> t_dist_vector;
174  t_dist_vector(i) = t_p(i) - t_n(i);
175  return sqrt(t_dist_vector(i) * t_dist_vector(i));
176  };
177 
178  map<EntityHandle, double> map_verts_length;
179 
180  for (auto f : surface_edges) {
181  int num_nodes;
182  const EntityHandle *conn_skin;
183  CHKERR moab.get_connectivity(f, conn_skin, num_nodes, true);
184  VectorDouble6 coords_skin(6);
185  if (th)
186  CHKERR moab.tag_get_data(th, conn_skin, num_nodes, &coords_skin[0]);
187  else
188  CHKERR moab.get_coords(conn_skin, num_nodes, &coords_skin[0]);
190  &coords_skin[0], &coords_skin[1], &coords_skin[2]);
192  &coords_skin[3], &coords_skin[4], &coords_skin[5]);
193  FTensor::Tensor1<double, 3> t_skin_delta;
194  t_skin_delta(i) = t_n1(i) - t_n0(i);
195  const double skin_edge_length = sqrt(t_skin_delta(i) * t_skin_delta(i));
196  for (int nn = 0; nn != 2; ++nn) {
197  auto it = map_verts_length.find(conn_skin[nn]);
198  if (it != map_verts_length.end())
199  it->second = std::min(it->second, skin_edge_length);
200  else
201  map_verts_length[conn_skin[nn]] = skin_edge_length;
202  }
203  }
204 
205  for (auto m : map_verts_length) {
206 
208  if (th)
209  CHKERR moab.tag_get_data(th, &m.first, 1, &t_n(0));
210  else
211  CHKERR moab.get_coords(&m.first, 1, &t_n(0));
212 
213  double min_dist = rel_tol * m.second;
214  FTensor::Tensor1<double, 3> t_min_coords;
215  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
216  &t_n(0), 1, fixed_edges, &min_dist, &t_min_coords(0));
217 
218  if (min_dist < rel_tol * m.second || min_dist < abs_tol) {
219  if (th)
220  CHKERR moab.tag_set_data(th, &m.first, 1, &t_min_coords(0));
221  else
222  CHKERR moab.set_coords(&m.first, 1, &t_min_coords(0));
223  }
224  }
225 
227 }
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:476
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:91
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:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

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

2006  {
2007  CoreInterface &m_field = cOre;
2008  moab::Interface &moab = m_field.get_moab();
2009  PrismInterface *interface;
2011  CHKERR m_field.getInterface(interface);
2012  EntityHandle meshset_volume;
2013  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
2014  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2015  split_bit, BitRefLevel().set(), MBTET, meshset_volume);
2016  EntityHandle meshset_surface;
2017  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
2018  CHKERR moab.add_entities(meshset_surface, ents);
2019  CHKERR interface->getSides(meshset_surface, split_bit, true);
2020  CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
2021  true);
2022  CHKERR moab.delete_entities(&meshset_volume, 1);
2023  CHKERR moab.delete_entities(&meshset_surface, 1);
2024  if (th) {
2025  Range prisms;
2026  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2027  bit, BitRefLevel().set(), MBPRISM, prisms);
2028  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
2029  int num_nodes;
2030  const EntityHandle *conn;
2031  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
2032  MatrixDouble data(3, 3);
2033  CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
2034  CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
2035  }
2036  }
2038 }
moab::Interface & get_moab()
Definition: Core.hpp:266
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
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:406

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

1813  {
1814  CoreInterface &m_field = cOre;
1815  moab::Interface &moab = m_field.get_moab();
1816  MeshRefinement *refiner;
1817  const RefEntity_multiIndex *ref_ents_ptr;
1819 
1820  CHKERR m_field.getInterface(refiner);
1821  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1822  CHKERR refiner->add_vertices_in_the_middel_of_edges(trimEdges, bit);
1823  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
1824  debug ? &trimEdges : NULL);
1825 
1826  trimNewVolumes.clear();
1827  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1828  bit, BitRefLevel().set(), MBTET, trimNewVolumes);
1829 
1830  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
1831  mit != edgesToTrim.end(); mit++) {
1832  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1833  boost::make_tuple(MBVERTEX, mit->first));
1834  if (vit ==
1835  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1836  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1837  "No vertex on trim edges, that make no sense");
1838  }
1839  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1840  if ((ref_ent->getBitRefLevel() & bit).any()) {
1841  EntityHandle vert = ref_ent->getRefEnt();
1842  trimNewVertices.insert(vert);
1843  verticesOnTrimEdges[vert] = mit->second;
1844  }
1845  }
1846 
1847  // Get faces which are trimmed
1848  trimNewSurfaces.clear();
1849  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1850  bit, bit, MBTRI, trimNewSurfaces);
1851 
1852  Range trim_new_surfaces_nodes;
1853  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
1854  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
1855  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
1856  Range faces_not_on_surface;
1857  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
1858  faces_not_on_surface, moab::Interface::UNION);
1859  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
1860 
1861  // Get surfaces which are not trimmed and add them to surface
1862  Range all_surfaces_on_bit_level;
1863  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1864  bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
1865  all_surfaces_on_bit_level =
1866  intersect(all_surfaces_on_bit_level, cutNewSurfaces);
1867 
1868  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
1869 
1870  Range trim_new_vertices_faces;
1871  CHKERR moab.get_adjacencies(trimNewVertices, 2, false,
1872  trim_new_vertices_faces, moab::Interface::UNION);
1873  trim_new_vertices_faces = intersect(trimNewSurfaces, trim_new_vertices_faces);
1874 
1876 }
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:476
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:595
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:406

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

279  {
280  CoreInterface &m_field = cOre;
281  moab::Interface &moab = m_field.get_moab();
283 
284  // trim mesh
285  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim_close, debug);
286  CHKERR trimEdgesInTheMiddle(trim_bit, debug);
287  if (fixed_edges) {
288  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
289  *fixed_edges);
290  }
291  if (corner_nodes) {
292  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
293  *corner_nodes);
294  }
295  if (update_meshsets) {
296  CHKERR m_field.getInterface<MeshsetsManager>()
297  ->updateAllMeshsetsByEntitiesChildren(trim_bit);
298  }
299 
300  // move nodes
302 
303  // remove faces
304  CHKERR trimSurface(fixed_edges, corner_nodes, debug);
305 
306  if (debug) {
308  Range bit2_edges;
309  CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
310  trim_bit, BitRefLevel().set(), MBEDGE, bit2_edges);
311  CHKERR SaveData(moab)("trim_fixed_edges.vtk",
312  intersect(*fixed_edges, bit2_edges));
313  }
314 
316 }
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:476
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:595
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:406
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 1878 of file CutMeshInterface.cpp.

1880  {
1881  CoreInterface &m_field = cOre;
1882  moab::Interface &moab = m_field.get_moab();
1884 
1885  Skinner skin(&moab);
1886  Range trim_tets_skin;
1887  CHKERR skin.find_skin(0, trimNewVolumes, false, trim_tets_skin);
1888  Range trim_tets_skin_edges;
1889  CHKERR moab.get_adjacencies(trim_tets_skin, 1, false, trim_tets_skin_edges,
1890  moab::Interface::UNION);
1891 
1892  Range barrier_vertices(trimNewVertices);
1893 
1894  if (fixed_edges && trimFixedEdges) {
1895 
1896  // get all vertices on fixed edges and surface
1897  Range trim_surface_edges;
1898  CHKERR moab.get_adjacencies(trimNewSurfaces, 1, false, trim_surface_edges,
1899  moab::Interface::UNION);
1900  Range fixed_edges_on_trim_surface;
1901  fixed_edges_on_trim_surface = intersect(trim_surface_edges, *fixed_edges);
1902  Range fixed_edges_on_trim_surface_verts;
1903  CHKERR moab.get_connectivity(fixed_edges_on_trim_surface,
1904  fixed_edges_on_trim_surface_verts, false);
1905 
1906  // get faces adjacent to barrier_vertices
1907  Range barrier_vertices_faces;
1908  CHKERR moab.get_adjacencies(barrier_vertices, 2, false,
1909  barrier_vertices_faces, moab::Interface::UNION);
1910  barrier_vertices_faces = intersect(barrier_vertices_faces, trimNewSurfaces);
1911 
1912  // get vertices on fixed edges
1913  Range fixed_edges_vertices;
1914  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_vertices, false);
1915  fixed_edges_vertices = intersect(barrier_vertices, fixed_edges_vertices);
1916  fixed_edges_vertices =
1917  subtract(fixed_edges_vertices, fixed_edges_on_trim_surface_verts);
1918  if (corner_nodes)
1919  fixed_edges_vertices.merge(intersect(barrier_vertices, *corner_nodes));
1920 
1921  // get faces adjacent to vertices on fixed edges
1922  Range fixed_edges_faces;
1923  CHKERR moab.get_adjacencies(fixed_edges_vertices, 2, false,
1924  fixed_edges_faces, moab::Interface::UNION);
1925  fixed_edges_faces = intersect(fixed_edges_faces, barrier_vertices_faces);
1926 
1927  if(debug && !fixed_edges_faces.empty())
1928  CHKERR SaveData(m_field.get_moab())("fixed_edges_faces.vtk",
1929  fixed_edges_faces);
1930 
1931  // get nodes on faces
1932  Range fixed_edges_faces_vertices;
1933  CHKERR moab.get_connectivity(fixed_edges_faces, fixed_edges_faces_vertices,
1934  false);
1935  barrier_vertices.merge(fixed_edges_faces_vertices);
1936  }
1937 
1938  if (debug && !barrier_vertices.empty())
1939  CHKERR SaveData(m_field.get_moab())("barrier_vertices.vtk",
1940  barrier_vertices);
1941 
1942  auto get_trim_skin_verts = [&]() {
1943  Range trim_surf_skin;
1944  CHKERR skin.find_skin(0, trimNewSurfaces, false, trim_surf_skin);
1945  trim_surf_skin = subtract(trim_surf_skin, trim_tets_skin_edges);
1946  Range trim_surf_skin_verts;
1947  CHKERR moab.get_connectivity(trim_surf_skin, trim_surf_skin_verts, true);
1948  trim_surf_skin_verts = subtract(trim_surf_skin_verts, barrier_vertices);
1949  return trim_surf_skin_verts;
1950  };
1951 
1952  int nn = 0;
1953 
1954  Range outside_verts;
1955  while (!(outside_verts = get_trim_skin_verts()).empty()) {
1956 
1957  if (debug && !trimNewSurfaces.empty())
1958  CHKERR SaveData(m_field.get_moab())(
1959  "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
1960  unite(trimNewSurfaces, outside_verts));
1961 
1962  Range outside_faces;
1963  CHKERR moab.get_adjacencies(outside_verts, 2, false, outside_faces,
1964  moab::Interface::UNION);
1965  trimNewSurfaces = subtract(trimNewSurfaces, outside_faces);
1966  ++nn;
1967  }
1968 
1970 }
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:476
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

Member Data Documentation

◆ aveLength

double MoFEM::CutMeshInterface::aveLength
private

Average edge length.

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

◆ cutFrontVolumes

Range MoFEM::CutMeshInterface::cutFrontVolumes
private

Definition at line 573 of file CutMeshInterface.hpp.

◆ cutNewSurfaces

Range MoFEM::CutMeshInterface::cutNewSurfaces
private

Definition at line 544 of file CutMeshInterface.hpp.

◆ cutNewVertices

Range MoFEM::CutMeshInterface::cutNewVertices
private

Definition at line 547 of file CutMeshInterface.hpp.

◆ cutNewVolumes

Range MoFEM::CutMeshInterface::cutNewVolumes
private

Definition at line 543 of file CutMeshInterface.hpp.

◆ cutSurfaceVolumes

Range MoFEM::CutMeshInterface::cutSurfaceVolumes
private

Definition at line 572 of file CutMeshInterface.hpp.

◆ cutVolumes

Range MoFEM::CutMeshInterface::cutVolumes
private

Definition at line 542 of file CutMeshInterface.hpp.

◆ edgesToCut

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

Definition at line 564 of file CutMeshInterface.hpp.

◆ edgesToTrim

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

Definition at line 566 of file CutMeshInterface.hpp.

◆ fRont

Range MoFEM::CutMeshInterface::fRont
private

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

◆ mergedSurfaces

Range MoFEM::CutMeshInterface::mergedSurfaces
private

Definition at line 555 of file CutMeshInterface.hpp.

◆ mergedVolumes

Range MoFEM::CutMeshInterface::mergedVolumes
private

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

◆ sUrface

Range MoFEM::CutMeshInterface::sUrface
private

Definition at line 533 of file CutMeshInterface.hpp.

◆ treeSurfPtr

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

Definition at line 538 of file CutMeshInterface.hpp.

◆ trimEdges

Range MoFEM::CutMeshInterface::trimEdges
private

Definition at line 552 of file CutMeshInterface.hpp.

◆ trimFixedEdges

bool MoFEM::CutMeshInterface::trimFixedEdges
private

Definition at line 536 of file CutMeshInterface.hpp.

◆ trimNewSurfaces

Range MoFEM::CutMeshInterface::trimNewSurfaces
private

Definition at line 551 of file CutMeshInterface.hpp.

◆ trimNewVertices

Range MoFEM::CutMeshInterface::trimNewVertices
private

Definition at line 550 of file CutMeshInterface.hpp.

◆ trimNewVolumes

Range MoFEM::CutMeshInterface::trimNewVolumes
private

Definition at line 549 of file CutMeshInterface.hpp.

◆ verticesOnCutEdges

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

Definition at line 565 of file CutMeshInterface.hpp.

◆ verticesOnTrimEdges

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

Definition at line 567 of file CutMeshInterface.hpp.

◆ vOlume

Range MoFEM::CutMeshInterface::vOlume
private

Definition at line 534 of file CutMeshInterface.hpp.

◆ zeroDistanceEnts

Range MoFEM::CutMeshInterface::zeroDistanceEnts
private

Definition at line 545 of file CutMeshInterface.hpp.

◆ zeroDistanceVerts

Range MoFEM::CutMeshInterface::zeroDistanceVerts
private

Definition at line 546 of file CutMeshInterface.hpp.


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