v0.9.1
Classes | Public Types | Public Member Functions | Public Attributes | Private Attributes | List of all members
MoFEM::CutMeshInterface Struct Reference

Interface to cut meshes. More...

#include <src/interfaces/CutMeshInterface.hpp>

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

Classes

struct  LengthMapData
 
struct  SaveData
 
struct  TreeData
 

Public Types

typedef multi_index_container< LengthMapData, indexed_by< ordered_non_unique< member< LengthMapData, double, &LengthMapData::lEngth > >, hashed_unique< member< LengthMapData, EntityHandle, &LengthMapData::eDge > >, ordered_non_unique< member< LengthMapData, double, &LengthMapData::qUality > > > > LengthMapData_multi_index
 

Public Member Functions

MoFEMErrorCode query_interface (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 setConstrainSurface (const Range surf)
 Set the constrain surface object. More...
 
MoFEMErrorCode mergeSurface (const Range surface)
 merge surface entities More...
 
MoFEMErrorCode mergeVolumes (const Range volume)
 merge volume entities More...
 
MoFEMErrorCode snapSurfaceSkinToEdges (const Range fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
 
MoFEMErrorCode snapSurfaceToEdges (const Range surface_edges, const Range fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
 
MoFEMErrorCode buildTree ()
 build tree More...
 
MoFEMErrorCode cutOnly (Range vol, const BitRefLevel cut_bit, Tag th, const double tol_geometry, const double tol_cut_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
 Cut mesh only. More...
 
MoFEMErrorCode trimOnly (const BitRefLevel trim_bit, Tag th, const double tol_geometry, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
 Trim mesh only. More...
 
MoFEMErrorCode cutAndTrim (int &first_bit, Tag th, const double tol_geometry, const double tol_cut_close, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
 Cut and trim. More...
 
MoFEMErrorCode cutTrimAndMerge (int &first_bit, const int fraction_level, Tag th, const double tol_geometry, const double tol_cut_close, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)
 Cut, trim and merge. More...
 
MoFEMErrorCode makeFront (const bool debug=false)
 Create front from the surface. More...
 
MoFEMErrorCode createSurfaceLevelSets (int verb=QUIET, const bool debug=false)
 Calculate distance from mesh nodes to cut surface. More...
 
MoFEMErrorCode createFrontLevelSets (Range vol, Tag th=nullptr, int verb=QUIET, const bool debug=false)
 Calculate distance from mesh nodes to surface front. More...
 
MoFEMErrorCode findLevelSetVolumes (Tag th, Range &vol_edges, const bool remove_adj_prims_edges, int verb=QUIET, const bool debug=false, const std::string edges_file_name=string())
 Find level set volumes. More...
 
MoFEMErrorCode findLevelSetVolumes (int verb=QUIET, const bool debug=false)
 Find level set volumes. More...
 
MoFEMErrorCode refineMesh (const int init_bit_level, const int surf_levels, const int front_levels, Range *fixed_edges=nullptr, int verb=QUIET, const bool debug=false)
 Refine and set level sets. More...
 
MoFEMErrorCode findEdgesToCut (Range vol, int verb=QUIET, const bool debug=false)
 find edges to cut More...
 
MoFEMErrorCode projectZeroDistanceEnts (Range *fixed_edges, Range *corner_nodes, const double geometry_tol=0, const double close_tol=0, const int verb=QUIET, const bool debug=false)
 Find entities on cut surface which can be projected. More...
 
MoFEMErrorCode cutEdgesInMiddle (const BitRefLevel bit, Range &cut_vols, Range &cut_surf, Range &cut_verts, const bool debug=false)
 cut edges More...
 
MoFEMErrorCode moveMidNodesOnCutEdges (Tag th=NULL)
 projecting of mid edge nodes on new mesh on surface More...
 
MoFEMErrorCode findEdgesToTrim (Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol_geometry=1e-4, const double tol=1e-4, const bool debug=false)
 Find edges to trimEdges. More...
 
MoFEMErrorCode trimEdgesInTheMiddle (const BitRefLevel bit, const bool debug=false)
 trim edges More...
 
MoFEMErrorCode moveMidNodesOnTrimmedEdges (Tag th=NULL)
 move trimmed edges mid nodes More...
 
MoFEMErrorCode trimSurface (Range *fixed_edge, Range *corner_nodes, const bool debug=false)
 Trim surface from faces beyond front. More...
 
MoFEMErrorCode removePathologicalFrontTris (const BitRefLevel split_bit, Range &ents)
 Remove pathological elements on surface internal front. More...
 
MoFEMErrorCode splitSides (const BitRefLevel split_bit, const BitRefLevel bit, const Range &ents, Tag th=NULL)
 split sides More...
 
MoFEMErrorCode mergeBadEdges (const int fraction_level, const Range &tets, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Range &merged_nodes, Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets=false, const BitRefLevel *bit_ptr=NULL, const bool debug=false)
 Merge edges. More...
 
MoFEMErrorCode mergeBadEdges (const int fraction_level, const BitRefLevel cut_bit, const BitRefLevel trim_bit, const BitRefLevel bit, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Tag th=NULL, const bool update_meshsets=false, const bool debug=false)
 Merge edges. More...
 
MoFEMErrorCode setTagData (Tag th, const BitRefLevel bit=BitRefLevel())
 set coords to tag More...
 
MoFEMErrorCode setCoords (Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
 set coords from tag More...
 
const Range & getVolume () const
 
const Range & getSurface () const
 
const Range & getFront () const
 
const Range & getCutEdges () const
 
const Range & getNewCutVolumes () const
 
const Range & getNewCutSurfaces () const
 
const Range & getNewCutVertices () const
 
const Range & projectZeroDistanceEnts () const
 
const Range & getTrimEdges () const
 
const Range & getNewTrimVolumes () const
 
const Range & getNewTrimSurfaces () const
 
const Range & getNewTrimVertices () const
 
const Range & getMergedVolumes () const
 
const Range & getMergedSurfaces () const
 
const Range & getCutSurfaceVolumes () const
 
const Range & getCutFrontVolumes () const
 
void setTrimFixedEdges (const bool b)
 
MoFEMErrorCode saveCutEdges (const std::string prefix="")
 
MoFEMErrorCode saveTrimEdges ()
 
boost::shared_ptr< OrientedBoxTreeTool > & getTreeSurfPtr ()
 
MoFEMErrorCode clearMap ()
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (const MOFEMuuid &uuid, bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE , bool VERIFY = false>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, IFACE *&iface) const
 Get interface by uuid and return reference to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 
virtual MoFEMErrorCode getLibVersion (Version &version) const
 Get library version. More...
 
virtual const MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version) const
 Get database major version. More...
 
virtual MoFEMErrorCode getInterfaceVersion (Version &version) const
 Get database major version. More...
 
template<>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, UnknownInterface *&iface) const
 

Public Attributes

MoFEM::CorecOre
 
int lineSearchSteps
 
int nbMaxMergingCycles
 
double projectEntitiesQualityTrashold
 

Private Attributes

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

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

226  {
227  CoreInterface &m_field = cOre;
228  moab::Interface &moab = m_field.get_moab();
230  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
231  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
234 }
virtual moab::Interface & get_moab()=0
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define CHKERR
Inline error check.
Definition: definitions.h:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

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

◆ 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 }
virtual moab::Interface & get_moab()=0
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:125
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:108
#define CHKERR
Inline error check.
Definition: definitions.h:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ createFrontLevelSets()

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

Calculate distance from mesh nodes to surface front.

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

Definition at line 568 of file CutMeshInterface.cpp.

570  {
571  CoreInterface &m_field = cOre;
572  moab::Interface &moab = m_field.get_moab();
574 
575  auto create_tag = [&](const std::string name, const int dim) {
576  Tag th;
577  rval = moab.tag_get_handle(name.c_str(), th);
578  if (rval == MB_SUCCESS)
579  return th;
580  std::vector<double> def_val(dim, 0);
581  CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
582  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
583  return th;
584  };
585 
586  Range vol_vertices;
587  CHKERR moab.get_connectivity(vol, vol_vertices, true);
588 
589  std::vector<double> min_distances_from_front(vol_vertices.size(), -1);
590  std::vector<double> points_on_edges(3 * vol_vertices.size(), 0);
591  std::vector<EntityHandle> closest_edges(vol_vertices.size(), 0);
592 
593  std::vector<double> coords(3 * vol_vertices.size());
594  if (th)
595  CHKERR moab.tag_get_data(th, vol_vertices, &*coords.begin());
596  else
597  CHKERR moab.get_coords(vol_vertices, &*coords.begin());
598 
599  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
600  &*coords.begin(), vol_vertices.size(), fRont,
601  &*min_distances_from_front.begin(), &*points_on_edges.begin(),
602  &*closest_edges.begin());
603 
604  if (!points_on_edges.empty()) {
605  for (int i = 0; i != min_distances_from_front.size(); ++i) {
606  Range faces;
607  CHKERR moab.get_adjacencies(&closest_edges[0], 1, 2, false, faces);
608  auto point_in = getVectorAdaptor(&coords[3 * i], 3);
609  auto point_out = getVectorAdaptor(&points_on_edges[3 * i], 3);
610  point_out -= point_in;
611  }
612  }
613 
614  auto th_dist_front_vec = create_tag("DIST_FRONT_VECTOR", 3);
615  CHKERR moab.tag_set_data(th_dist_front_vec, vol_vertices,
616  &*points_on_edges.begin());
617 
619 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
const int dim
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:42
#define CHKERR
Inline error check.
Definition: definitions.h:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26

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

480  {
481  CoreInterface &m_field = cOre;
482  moab::Interface &moab = m_field.get_moab();
484  auto tools_interface = m_field.getInterface<Tools>();
485 
486  auto create_tag = [&](const std::string name, const int dim) {
487  Tag th;
488  rval = moab.tag_get_handle(name.c_str(), th);
489  if (rval == MB_SUCCESS)
490  return th;
491  std::vector<double> def_val(dim, 0);
492  CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
493  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
494 
495  return th;
496  };
497 
498  auto set_vol = [&](const Range &vol_verts, std::vector<double> &coords,
499  std::vector<double> &dist_surface_vec,
500  std::vector<double> &dist_surface_normal_vec,
501  std::vector<double> &dist_surface_digned_dist_vec) {
503 
504  coords.resize(3 * vol_verts.size());
505  dist_surface_vec.resize(3 * vol_verts.size());
506  dist_surface_normal_vec.resize(3 * vol_verts.size());
507  dist_surface_digned_dist_vec.resize(vol_verts.size());
508 
509  CHKERR moab.get_coords(vol_verts, &*coords.begin());
510  std::srand(0);
511 
512  for (auto v : vol_verts) {
513 
514  const int index = vol_verts.index(v);
515  auto point_in = getVectorAdaptor(&coords[3 * index], 3);
516  VectorDouble3 point_out(3);
517  EntityHandle facets_out;
518  CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
519  &point_out[0], facets_out);
520 
521  VectorDouble3 n(3);
522  CHKERR tools_interface->getTriNormal(facets_out, &*n.begin());
523  n /= norm_2(n);
524 
525  VectorDouble3 delta = point_out - point_in;
526  const double dot = inner_prod(delta, n);
527 
528  if (std::abs(dot) < std::numeric_limits<double>::epsilon()) {
529  if (std::rand() % 2 == 0)
530  delta += n * std::numeric_limits<double>::epsilon();
531  else
532  delta -= n * std::numeric_limits<double>::epsilon();
533  }
534 
535  auto dist_vec = getVectorAdaptor(&dist_surface_vec[3 * index], 3);
536  noalias(dist_vec) = delta;
537 
538  auto dist_normal_vec =
539  getVectorAdaptor(&dist_surface_normal_vec[3 * index], 3);
540 
541  dist_surface_digned_dist_vec[index] = dot;
542  noalias(dist_normal_vec) = dot * n;
543  }
544 
546  };
547 
548  std::vector<double> coords;
549  std::vector<double> dist_surface_vec;
550  std::vector<double> dist_surface_normal_vec;
551  std::vector<double> dist_surface_signed_dist_vec;
552  Range vol_verts;
553  CHKERR moab.get_connectivity(vOlume, vol_verts, true);
554 
555  CHKERR set_vol(vol_verts, coords, dist_surface_vec, dist_surface_normal_vec,
556  dist_surface_signed_dist_vec);
557 
558  CHKERR moab.tag_set_data(create_tag("DIST_SURFACE_VECTOR", 3), vol_verts,
559  &*dist_surface_vec.begin());
560  CHKERR moab.tag_set_data(create_tag("DIST_SURFACE_NORMAL_VECTOR", 3),
561  vol_verts, &*dist_surface_normal_vec.begin());
562  CHKERR moab.tag_set_data(create_tag("DIST_SURFACE_NORMAL_SIGNED", 1),
563  vol_verts, &*dist_surface_signed_dist_vec.begin());
564 
566 }
virtual moab::Interface & get_moab()=0
static constexpr double delta
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
const int dim
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:42
#define CHKERR
Inline error check.
Definition: definitions.h:601
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ cutAndTrim()

MoFEMErrorCode MoFEM::CutMeshInterface::cutAndTrim ( int &  first_bit,
Tag  th,
const double  tol_geometry,
const double  tol_cut_close,
const double  tol_trim_close,
Range *  fixed_edges = NULL,
Range *  corner_nodes = NULL,
const bool  update_meshsets = false,
const bool  debug = false 
)

Cut and trim.

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

Definition at line 322 of file CutMeshInterface.cpp.

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

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

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

◆ cutOnly()

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

Cut mesh only.

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

Definition at line 237 of file CutMeshInterface.cpp.

240  {
241  CoreInterface &m_field = cOre;
242  moab::Interface &moab = m_field.get_moab();
244 
245  // cut mesh
247  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_geometry,
248  tol_cut_close, QUIET, debug);
251 
252  CHKERR cOre.getInterface<BitRefManager>()->updateRange(constrainSurface,
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 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
MoFEMErrorCode cutEdgesInMiddle(const BitRefLevel bit, Range &cut_vols, Range &cut_surf, Range &cut_verts, const bool debug=false)
cut edges
MoFEMErrorCode findEdgesToCut(Range vol, 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:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ cutTrimAndMerge()

MoFEMErrorCode MoFEM::CutMeshInterface::cutTrimAndMerge ( int &  first_bit,
const int  fraction_level,
Tag  th,
const double  tol_geometry,
const double  tol_cut_close,
const double  tol_trim_close,
Range &  fixed_edges,
Range &  corner_nodes,
const bool  update_meshsets = false,
const bool  debug = false 
)

Cut, trim and merge.

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

Definition at line 392 of file CutMeshInterface.cpp.

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

◆ findEdgesToCut()

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

find edges to cut

Parameters
volis tetrahedrons search to cut
verbverbosity level
debugdebugging
Returns
MoFEMErrorCode

Definition at line 868 of file CutMeshInterface.cpp.

869  {
870  CoreInterface &m_field = cOre;
871  moab::Interface &moab = m_field.get_moab();
873 
874  edgesToCut.clear();
875  cutEdges.clear();
876 
877  Tag th_signed_dist;
878  CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_SIGNED", th_signed_dist);
879 
880  auto get_tag_edge_dist = [&](auto th, auto conn) {
881  std::array<double, 2> r;
882  CHKERR moab.tag_get_data(th, conn, 2, r.data());
883  return r;
884  };
885 
886  auto get_conn = [&](const auto e) {
887  int num_nodes;
888  const EntityHandle *conn;
889  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
890  return conn;
891  };
892 
893  auto get_adj = [&moab](const Range r, const int dim) {
894  Range a;
895  if (dim)
896  CHKERR moab.get_adjacencies(r, dim, false, a, moab::Interface::UNION);
897  else
898  CHKERR moab.get_connectivity(r, a, true);
899  return a;
900  };
901 
902  auto get_range = [](std::vector<EntityHandle> v) {
903  Range r;
904  r.insert_list(v.begin(), v.end());
905  return r;
906  };
907 
908  auto vol_edges = get_adj(vol, 1);
909 
910  aveLength = 0;
911  maxLength = 0;
912  int nb_ave_length = 0;
913  for (auto e : vol_edges) {
914 
915  auto conn = get_conn(e);
916  auto dist = get_tag_edge_dist(th_signed_dist, conn);
917  const auto dist_max = std::max(dist[0], dist[1]);
918  const auto dist_min = std::min(dist[0], dist[1]);
919  const auto dot = dist[0] * dist[1];
920 
921  VectorDouble6 coords(6);
922  CHKERR moab.get_coords(conn, 2, &coords[0]);
923  VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
924  VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
925  VectorDouble3 ray = n1 - n0;
926  const double ray_length = norm_2(ray);
927  ray /= ray_length;
928 
929  if (
930 
931  (dot < 0 && dist_max > std::numeric_limits<double>::epsilon()) ||
932  (std::abs(dist_min) < std::numeric_limits<double>::epsilon() &&
933  dist_max > std::numeric_limits<double>::epsilon())
934 
935  ) {
936 
937  // Edges is on two sides of the surface
938  const double s = dist[0] / (dist[0] - dist[1]);
939  const double dist = s * ray_length;
940 
941  auto add_edge = [&](auto dist) {
942  edgesToCut[e].dIst = dist;
943  edgesToCut[e].lEngth = ray_length;
944  edgesToCut[e].unitRayDir = ray;
945  edgesToCut[e].rayPoint = n0;
946  cutEdges.insert(e);
947 
948  aveLength += norm_2(ray);
949  maxLength = fmax(maxLength, norm_2(ray));
950  ++nb_ave_length;
951  };
952 
953  add_edge(dist);
954  }
955  }
956  aveLength /= nb_ave_length;
957 
958  if (debug)
959  CHKERR SaveData(m_field.get_moab())("cut_edges.vtk", cutEdges);
960 
961  if (debug)
962  CHKERR SaveData(m_field.get_moab())("cut_edges_vol.vtk",
963  get_adj(cutEdges, 3));
964 
966 }
virtual moab::Interface & get_moab()=0
double maxLength
Maximal edge length.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
double aveLength
Average edge length.
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:94
const int dim
static const bool debug
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:42
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:108
#define CHKERR
Inline error check.
Definition: definitions.h:601
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
map< EntityHandle, TreeData > edgesToCut

◆ findEdgesToTrim()

MoFEMErrorCode MoFEM::CutMeshInterface::findEdgesToTrim ( Range *  fixed_edges,
Range *  corner_nodes,
Tag  th = NULL,
const double  tol_geometry = 1e-4,
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 1482 of file CutMeshInterface.cpp.

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

◆ findLevelSetVolumes() [1/2]

MoFEMErrorCode MoFEM::CutMeshInterface::findLevelSetVolumes ( Tag  th,
Range &  vol_edges,
const bool  remove_adj_prims_edges,
int  verb = QUIET,
const bool  debug = false,
const std::string  edges_file_name = string() 
)

Find level set volumes.

Parameters
th
vol_edges
remove_adj_prims_edges
verb
debug
edges_file_name
Returns
MoFEMErrorCode

Definition at line 621 of file CutMeshInterface.cpp.

623  {
624  CoreInterface &m_field = cOre;
625  moab::Interface &moab = m_field.get_moab();
627 
628  auto get_tag_dim = [&](auto th) {
629  int dim;
630  CHKERR moab.tag_get_length(th, dim);
631  return dim;
632  };
633  auto dim = get_tag_dim(th);
634 
635  auto get_tag_data = [&](auto th, auto conn) {
636  const void *ptr;
637  CHKERR moab.tag_get_by_ptr(th, &conn, 1, &ptr);
638  return getVectorAdaptor(
639  const_cast<double *>(static_cast<const double *>(ptr)), 3);
640  };
641 
642  auto get_edge_ray = [&](auto conn) {
643  VectorDouble6 coords(6);
644  CHKERR moab.get_coords(conn, 2, &coords[0]);
645  VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
646  VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
647  VectorDouble3 ray = n1 - n0;
648  return ray;
649  };
650 
651  Range edges;
652  CHKERR moab.get_adjacencies(vOlume, 1, true, edges, moab::Interface::UNION);
653 
654  auto remove_prisms_edges = [&](Range &edges) {
656  Range prisms;
657  CHKERR moab.get_adjacencies(edges, 3, false, prisms,
658  moab::Interface::UNION);
659  prisms = prisms.subset_by_type(MBPRISM);
660  Range prisms_verts;
661  CHKERR moab.get_connectivity(prisms, prisms_verts, true);
662  Range prism_edges;
663  CHKERR moab.get_adjacencies(prisms_verts, 1, false, prism_edges,
664  moab::Interface::UNION);
665  edges = subtract(edges, prism_edges);
667  };
668  if (remove_adj_prims_edges)
669  CHKERR remove_prisms_edges(edges);
670 
671  std::vector<EntityHandle> cut_edges;
672  cut_edges.reserve(edges.size());
673 
674  auto get_cut_edges_vec = [&](auto th, auto &cut_edges, auto e, auto &&conn) {
676  auto ray = get_edge_ray(conn);
677  const double length = norm_2(ray);
678  ray /= length;
679  auto signed_norm = [&](const auto &v) { return inner_prod(ray, v); };
680  const auto dist0 = get_tag_data(th, conn[0]);
681  const auto dist1 = get_tag_data(th, conn[1]);
682  // const double min_dist = std::min(norm_2(dist0), norm_2(dist1));
683  // if (min_dist < 2 * length) {
684  auto opposite = inner_prod(dist0, dist1);
685  if (opposite <= std::numeric_limits<double>::epsilon()) {
686  const double sign_dist0 = signed_norm(dist0);
687  const double sign_dist1 = signed_norm(dist1);
688  if (sign_dist0 > -std::numeric_limits<double>::epsilon() &&
689  sign_dist1 < std::numeric_limits<double>::epsilon())
690  cut_edges.push_back(e);
691  }
692  // }
694  };
695 
696  auto get_cut_edges_signed_dist = [&](auto th, auto &cut_edges, auto e,
697  auto &&conn) {
699  auto get_tag_signed_dist = [&](auto conn) {
700  double dist;
701  CHKERR moab.tag_get_data(th, &conn, 1, &dist);
702  return dist;
703  };
704  const auto dist0 = get_tag_signed_dist(conn[0]);
705  const auto dist1 = get_tag_signed_dist(conn[1]);
706  const double opposite = dist0 * dist1;
707  if (opposite < 0)
708  cut_edges.push_back(e);
710  };
711 
712  auto get_conn = [&](auto e) {
713  int num_nodes;
714  const EntityHandle *conn;
715  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
716  return conn;
717  };
718 
719  if (dim == 3)
720  for (auto e : edges)
721  CHKERR get_cut_edges_vec(th, cut_edges, e, get_conn(e));
722  else
723  for (auto e : edges)
724  CHKERR get_cut_edges_signed_dist(th, cut_edges, e, get_conn(e));
725 
726  CHKERR moab.get_adjacencies(&*cut_edges.begin(), cut_edges.size(), 3, false,
727  vol_edges, moab::Interface::UNION);
728  vol_edges = intersect(vol_edges, vOlume);
729 
730  auto convert_to_ranger = [](auto &v) {
731  Range r;
732  r.insert_list(v.begin(), v.end());
733  return r;
734  };
735 
736  if (debug && !edges_file_name.empty())
737  CHKERR SaveData(m_field.get_moab())(edges_file_name,
738  convert_to_ranger(cut_edges));
739 
741 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:94
const int dim
static const bool debug
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:42
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:108
#define CHKERR
Inline error check.
Definition: definitions.h:601
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ findLevelSetVolumes() [2/2]

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

Find level set volumes.

Parameters
update_front
verb
debug
Returns
MoFEMErrorCode

Definition at line 743 of file CutMeshInterface.cpp.

744  {
745  CoreInterface &m_field = cOre;
746  moab::Interface &moab = m_field.get_moab();
748 
749  CHKERR createFrontLevelSets(vOlume, nullptr, verb, debug);
750  Tag th_dist_front_vec;
751  CHKERR moab.tag_get_handle("DIST_FRONT_VECTOR", th_dist_front_vec);
752  CHKERR findLevelSetVolumes(th_dist_front_vec, cutFrontVolumes, true, verb,
753  debug, "cutFrontEdges.vtk");
754 
756 
757  Tag th_dist_surface_vec;
758  CHKERR moab.tag_get_handle("DIST_SURFACE_VECTOR", th_dist_surface_vec);
759  cutSurfaceVolumes.clear();
760  CHKERR findLevelSetVolumes(th_dist_surface_vec, cutSurfaceVolumes, true, verb,
761  debug, "cutSurfaceEdges.vtk");
762  Tag th_dist_signed_normal;
763  CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_SIGNED",
764  th_dist_signed_normal);
765  Range signed_vol;
766  CHKERR findLevelSetVolumes(th_dist_signed_normal, signed_vol, true, verb,
767  debug, "cutSurfaceEdges.vtk");
768  cutSurfaceVolumes = intersect(cutSurfaceVolumes, signed_vol);
769 
770  if (debug)
771  CHKERR SaveData(m_field.get_moab())("level_sets.vtk", vOlume);
772  if (debug)
773  CHKERR SaveData(m_field.get_moab())("cutSurfaceVolumes.vtk",
775  if (debug)
776  CHKERR SaveData(m_field.get_moab())("cutFrontVolumes.vtk", cutFrontVolumes);
777 
779 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
MoFEMErrorCode createFrontLevelSets(Range vol, Tag th=nullptr, int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to surface front.
MoFEMErrorCode findLevelSetVolumes(Tag th, Range &vol_edges, const bool remove_adj_prims_edges, int verb=QUIET, const bool debug=false, const std::string edges_file_name=string())
Find level set volumes.

◆ getCutEdges()

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

Definition at line 484 of file CutMeshInterface.hpp.

484 { return cutEdges; }

◆ getCutFrontVolumes()

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

Definition at line 501 of file CutMeshInterface.hpp.

◆ getCutSurfaceVolumes()

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

Definition at line 500 of file CutMeshInterface.hpp.

◆ getFront()

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

Definition at line 482 of file CutMeshInterface.hpp.

482 { return fRont; }

◆ getMergedSurfaces()

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

Definition at line 498 of file CutMeshInterface.hpp.

◆ getMergedVolumes()

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

Definition at line 497 of file CutMeshInterface.hpp.

497 { return mergedVolumes; }

◆ getNewCutSurfaces()

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

Definition at line 486 of file CutMeshInterface.hpp.

◆ getNewCutVertices()

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

Definition at line 487 of file CutMeshInterface.hpp.

◆ getNewCutVolumes()

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

Definition at line 485 of file CutMeshInterface.hpp.

485 { return cutNewVolumes; }

◆ getNewTrimSurfaces()

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

Definition at line 494 of file CutMeshInterface.hpp.

◆ getNewTrimVertices()

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

Definition at line 495 of file CutMeshInterface.hpp.

◆ getNewTrimVolumes()

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

Definition at line 493 of file CutMeshInterface.hpp.

◆ getOptions()

MoFEMErrorCode MoFEM::CutMeshInterface::getOptions ( )

Get options from command line.

Returns
error code

Definition at line 51 of file CutMeshInterface.hpp.

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

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

481 { return sUrface; }

◆ getTreeSurfPtr()

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

Definition at line 509 of file CutMeshInterface.hpp.

509  {
510  return treeSurfPtr;
511  }
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr

◆ getTrimEdges()

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

Definition at line 492 of file CutMeshInterface.hpp.

492 { return trimEdges; }

◆ getVolume()

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

Definition at line 480 of file CutMeshInterface.hpp.

480 { 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 461 of file CutMeshInterface.cpp.

461  {
462  CoreInterface &m_field = cOre;
463  moab::Interface &moab = m_field.get_moab();
465  Skinner skin(&moab);
466  Range tets_skin;
467  CHKERR skin.find_skin(0, vOlume, false, tets_skin);
468  Range tets_skin_edges;
469  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
470  moab::Interface::UNION);
471  Range surface_skin;
472  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
473  fRont = subtract(surface_skin, tets_skin_edges);
474  if (debug)
475  CHKERR SaveData(moab)("front.vtk", fRont);
477 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

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

2209  {
2210  CoreInterface &m_field = cOre;
2211  moab::Interface &moab = m_field.get_moab();
2213 
2214  /**
2215  * \brief Merge nodes
2216  */
2217  struct MergeNodes {
2218  CoreInterface &mField;
2219  const bool onlyIfImproveQuality;
2220  Tag tH;
2221  bool updateMehsets;
2222 
2223  MergeNodes(CoreInterface &m_field, const bool only_if_improve_quality,
2224  Tag th, bool update_mehsets)
2225  : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
2226  tH(th), updateMehsets(update_mehsets) {
2227  mField.getInterface(nodeMergerPtr);
2228  }
2229  NodeMergerInterface *nodeMergerPtr;
2230  MoFEMErrorCode mergeNodes(int line_search, EntityHandle father,
2231  EntityHandle mother, Range &proc_tets,
2232  bool add_child = true) {
2233  moab::Interface &moab = mField.get_moab();
2235  const EntityHandle conn[] = {father, mother};
2236  Range vert_tets;
2237  CHKERR moab.get_adjacencies(conn, 2, 3, false, vert_tets,
2238  moab::Interface::UNION);
2239  vert_tets = intersect(vert_tets, proc_tets);
2240  Range out_tets;
2241  CHKERR nodeMergerPtr->mergeNodes(father, mother, out_tets, &vert_tets,
2242  onlyIfImproveQuality, 0, line_search,
2243  tH);
2244 
2245  if (add_child && nodeMergerPtr->getSuccessMerge()) {
2246 
2247  Range::iterator lo, hi = proc_tets.begin();
2248  for (auto pt = vert_tets.pair_begin(); pt != vert_tets.pair_end();
2249  ++pt) {
2250  lo = proc_tets.lower_bound(hi, proc_tets.end(), pt->first);
2251  if (lo != proc_tets.end()) {
2252  hi = proc_tets.upper_bound(lo, proc_tets.end(), pt->second);
2253  proc_tets.erase(lo, hi);
2254  } else
2255  break;
2256  }
2257  proc_tets.merge(out_tets);
2258 
2259  auto &parent_child_map = nodeMergerPtr->getParentChildMap();
2260 
2261  struct ChangeChild {
2262  EntityHandle child;
2263  ChangeChild(const EntityHandle child) : child(child) {}
2264  void operator()(NodeMergerInterface::ParentChild &p) {
2265  p.cHild = child;
2266  }
2267  };
2268 
2269  std::vector<decltype(parentsChildMap.get<0>().begin())> it_vec;
2270  it_vec.reserve(parentsChildMap.size());
2271 
2272  for (auto &p : parent_child_map) {
2273 
2274  it_vec.clear();
2275  for (auto it = parentsChildMap.get<0>().equal_range(p.pArent);
2276  it.first != it.second; ++it.first)
2277  it_vec.emplace_back(it.first);
2278 
2279  for (auto it = parentsChildMap.get<1>().equal_range(p.pArent);
2280  it.first != it.second; ++it.first)
2281  it_vec.emplace_back(parentsChildMap.project<0>(it.first));
2282 
2283  for (auto &it : it_vec)
2284  parentsChildMap.modify(it, ChangeChild(p.cHild));
2285  }
2286 
2287  parentsChildMap.insert(parent_child_map.begin(),
2288  parent_child_map.end());
2289  }
2291  }
2292 
2293  MoFEMErrorCode updateRangeByChilds(Range &new_surf, Range &edges_to_merge,
2294  Range &not_merged_edges,
2295  bool add_child) {
2296  moab::Interface &moab = mField.get_moab();
2298  if (add_child) {
2299 
2300  std::vector<EntityHandle> parents_ents_vec(parentsChildMap.size());
2301  for (auto &it : parentsChildMap)
2302  parents_ents_vec.emplace_back(it.pArent);
2303  Range parent_ents;
2304  parent_ents.insert_list(parents_ents_vec.begin(),
2305  parents_ents_vec.end());
2306 
2307  Range surf_parent_ents = intersect(new_surf, parent_ents);
2308  new_surf = subtract(new_surf, surf_parent_ents);
2309  Range child_surf_ents;
2310  CHKERR updateRangeByChilds(parentsChildMap, surf_parent_ents,
2311  child_surf_ents);
2312  new_surf.merge(child_surf_ents);
2313 
2314  Range edges_to_merge_parent_ents =
2315  intersect(edges_to_merge, parent_ents);
2316  edges_to_merge = subtract(edges_to_merge, edges_to_merge_parent_ents);
2317  Range merged_child_edge_ents;
2318  CHKERR updateRangeByChilds(parentsChildMap, edges_to_merge_parent_ents,
2319  merged_child_edge_ents);
2320 
2321  Range not_merged_edges_child_ents =
2322  intersect(not_merged_edges, parent_ents);
2323  not_merged_edges =
2324  subtract(not_merged_edges, not_merged_edges_child_ents);
2325  Range not_merged_child_edge_ents;
2326  CHKERR updateRangeByChilds(parentsChildMap, not_merged_edges_child_ents,
2327  not_merged_child_edge_ents);
2328 
2329  merged_child_edge_ents =
2330  subtract(merged_child_edge_ents, not_merged_child_edge_ents);
2331  edges_to_merge.merge(merged_child_edge_ents);
2332  not_merged_edges.merge(not_merged_child_edge_ents);
2333 
2334  if (updateMehsets) {
2336  (*mField.getInterface<MeshsetsManager>()), cubit_it)) {
2337  EntityHandle cubit_meshset = cubit_it->meshset;
2338  Range meshset_ents;
2339  CHKERR moab.get_entities_by_handle(cubit_meshset, meshset_ents,
2340  true);
2341  Range child_ents;
2342  CHKERR updateRangeByChilds(parentsChildMap, meshset_ents,
2343  child_ents);
2344  CHKERR moab.add_entities(cubit_meshset, child_ents);
2345  }
2346  }
2347  }
2348 
2350  };
2351 
2352  private:
2353  NodeMergerInterface::ParentChildMap parentsChildMap;
2354  std::vector<EntityHandle> childsVec;
2355 
2356  inline MoFEMErrorCode updateRangeByChilds(
2357  const NodeMergerInterface::ParentChildMap &parent_child_map,
2358  const Range &parents, Range &childs) {
2360  childsVec.clear();
2361  childsVec.reserve(parents.size());
2362  for (auto pit = parents.pair_begin(); pit != parents.pair_end(); pit++) {
2363  auto it = parent_child_map.lower_bound(pit->first);
2364  if (it != parent_child_map.end()) {
2365  for (auto hi_it = parent_child_map.upper_bound(pit->second);
2366  it != hi_it; ++it)
2367  childsVec.emplace_back(it->cHild);
2368  }
2369  }
2370  childs.insert_list(childsVec.begin(), childsVec.end());
2372  }
2373  };
2374 
2375  /**
2376  * \brief Calculate edge length
2377  */
2378  struct LengthMap {
2379  Tag tH;
2380  CoreInterface &mField;
2381  moab::Interface &moab;
2382  const double maxLength;
2383  LengthMap(CoreInterface &m_field, Tag th, double max_length)
2384  : tH(th), mField(m_field), moab(m_field.get_moab()),
2385  maxLength(max_length) {
2386  gammaL = 1.;
2387  gammaQ = 3.;
2388  acceptedThrasholdMergeQuality = 0.5;
2389  }
2390 
2391  double
2392  gammaL; ///< Controls importance of length when ranking edges for merge
2393  double
2394  gammaQ; ///< Controls importance of quality when ranking edges for merge
2395  double acceptedThrasholdMergeQuality; ///< Do not merge quality if quality
2396  ///< above accepted thrashold
2397 
2398  MoFEMErrorCode operator()(const Range &tets, const Range &edges,
2399  LengthMapData_multi_index &length_map,
2400  double &ave) const {
2401  int num_nodes;
2402  const EntityHandle *conn;
2403  std::array<double, 12> coords;
2405  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
2406  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
2407  VectorDouble3 delta(3);
2408 
2409  struct NodeQuality {
2410  EntityHandle node;
2411  double quality;
2412  NodeQuality(const EntityHandle node) : node(node), quality(1) {}
2413  };
2414 
2415  typedef multi_index_container<
2416  NodeQuality, indexed_by<
2417 
2418  sequenced<>,
2419 
2420  hashed_non_unique<tag<Ent_mi_tag>,
2421  member<NodeQuality, EntityHandle,
2422  &NodeQuality::node>>
2423 
2424  >>
2425  NodeQuality_sequence;
2426 
2427  NodeQuality_sequence node_quality_sequence;
2428 
2429  Range edges_nodes;
2430  CHKERR moab.get_connectivity(tets, edges_nodes, false);
2431  Range edges_tets;
2432  CHKERR moab.get_adjacencies(edges, 3, false, edges_tets,
2433  moab::Interface::UNION);
2434  edges_tets = intersect(edges_tets, tets);
2435 
2436  for (auto node : edges_nodes)
2437  node_quality_sequence.emplace_back(node);
2438 
2439  for (auto tet : edges_tets) {
2440 
2441  CHKERR moab.get_connectivity(tet, conn, num_nodes, true);
2442  if (tH)
2443  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2444  else
2445  CHKERR moab.get_coords(conn, num_nodes, coords.data());
2446 
2447  const double q = Tools::volumeLengthQuality(coords.data());
2448 
2449  for (auto n : {0, 1, 2, 3}) {
2450  auto it = node_quality_sequence.get<1>().find(conn[n]);
2451  if (it != node_quality_sequence.get<1>().end())
2452  const_cast<double &>(it->quality) = std::min(q, it->quality);
2453  }
2454  }
2455 
2456  for (auto edge : edges) {
2457 
2458  CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
2459 
2460  if (tH)
2461  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2462  else
2463  CHKERR moab.get_coords(conn, num_nodes, coords.data());
2464 
2465  double q = 1;
2466  for (auto n : {0, 1}) {
2467  auto it = node_quality_sequence.get<1>().find(conn[n]);
2468  if (it != node_quality_sequence.get<1>().end())
2469  q = std::min(q, it->quality);
2470  }
2471 
2472  if (q < acceptedThrasholdMergeQuality) {
2473  noalias(delta) = (s0 - s1) / maxLength;
2474  double dot = inner_prod(delta, delta);
2475  double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
2476  length_map.insert(LengthMapData(val, q, edge));
2477  }
2478  }
2479 
2480  ave = 0;
2481  for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
2482  length_map.get<0>().begin();
2483  mit != length_map.get<0>().end(); mit++) {
2484  ave += mit->qUality;
2485  }
2486  ave /= length_map.size();
2488  }
2489  };
2490 
2491  /**
2492  * \brief Topological relations
2493  */
2494  struct Toplogy {
2495 
2496  CoreInterface &mField;
2497  Tag tH;
2498  const double tOL;
2499  Toplogy(CoreInterface &m_field, Tag th, const double tol)
2500  : mField(m_field), tH(th), tOL(tol) {}
2501 
2502  enum TYPE {
2503  FREE = 0,
2504  SKIN = 1 << 0,
2505  SURFACE = 1 << 1,
2506  SURFACE_SKIN = 1 << 2,
2507  FRONT_ENDS = 1 << 3,
2508  FIX_EDGES = 1 << 4,
2509  FIX_CORNERS = 1 << 5
2510  };
2511 
2512  typedef map<int, Range> SetsMap;
2513 
2514  MoFEMErrorCode classifyVerts(const Range &surface, const Range &tets,
2515  const Range &fixed_edges,
2516  const Range &corner_nodes,
2517  const Range &constrain_surface,
2518  SetsMap &sets_map) const {
2519  moab::Interface &moab(mField.get_moab());
2520  Skinner skin(&moab);
2522 
2523  sets_map[FIX_CORNERS].merge(corner_nodes);
2524  Range fixed_verts;
2525  CHKERR moab.get_connectivity(fixed_edges, fixed_verts, true);
2526  sets_map[FIX_EDGES].swap(fixed_verts);
2527 
2528  Range tets_skin;
2529  CHKERR skin.find_skin(0, tets, false, tets_skin);
2530  Range tets_skin_edges;
2531  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2532  moab::Interface::UNION);
2533  tets_skin.merge(constrain_surface);
2534 
2535  // surface skin
2536  Range surface_skin;
2537  CHKERR skin.find_skin(0, surface, false, surface_skin);
2538  Range front_in_the_body;
2539  front_in_the_body = subtract(surface_skin, tets_skin_edges);
2540  Range front_ends;
2541  CHKERR skin.find_skin(0, front_in_the_body, false, front_ends);
2542  sets_map[FRONT_ENDS].swap(front_ends);
2543 
2544  Range surface_skin_verts;
2545  CHKERR moab.get_connectivity(surface_skin, surface_skin_verts, true);
2546  sets_map[SURFACE_SKIN].swap(surface_skin_verts);
2547 
2548  // surface
2549  Range surface_verts;
2550  CHKERR moab.get_connectivity(surface, surface_verts, true);
2551  sets_map[SURFACE].swap(surface_verts);
2552 
2553  // skin
2554  Range tets_skin_verts;
2555  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
2556  sets_map[SKIN].swap(tets_skin_verts);
2557 
2558  Range tets_verts;
2559  CHKERR moab.get_connectivity(tets, tets_verts, true);
2560  sets_map[FREE].swap(tets_verts);
2561 
2563  }
2564 
2565  MoFEMErrorCode getProcTets(const Range &tets, const Range &edges_to_merge,
2566  Range &proc_tets) const {
2567  moab::Interface &moab(mField.get_moab());
2569  Range edges_to_merge_verts;
2570  CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts, true);
2571  Range edges_to_merge_verts_tets;
2572  CHKERR moab.get_adjacencies(edges_to_merge_verts, 3, false,
2573  edges_to_merge_verts_tets,
2574  moab::Interface::UNION);
2575  edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
2576  proc_tets.swap(edges_to_merge_verts_tets);
2578  }
2579 
2580  MoFEMErrorCode removeBadEdges(const Range &surface, const Range &tets,
2581  const Range &fixed_edges,
2582  const Range &corner_nodes,
2583  Range &edges_to_merge,
2584  Range &not_merged_edges) {
2585  moab::Interface &moab(mField.get_moab());
2587 
2588  // find skin
2589  Skinner skin(&moab);
2590  Range tets_skin;
2591  CHKERR skin.find_skin(0, tets, false, tets_skin);
2592  Range surface_skin;
2593  CHKERR skin.find_skin(0, surface, false, surface_skin);
2594 
2595  // end nodes
2596  Range tets_skin_edges;
2597  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2598  moab::Interface::UNION);
2599  Range surface_front;
2600  surface_front = subtract(surface_skin, tets_skin_edges);
2601  Range surface_front_nodes;
2602  CHKERR moab.get_connectivity(surface_front, surface_front_nodes, true);
2603  Range ends_nodes;
2604  CHKERR skin.find_skin(0, surface_front, false, ends_nodes);
2605 
2606  // remove bad merges
2607 
2608  // get surface and body skin verts
2609  Range surface_edges;
2610  CHKERR moab.get_adjacencies(surface, 1, false, surface_edges,
2611  moab::Interface::UNION);
2612  // get nodes on the surface
2613  Range surface_edges_verts;
2614  CHKERR moab.get_connectivity(surface_edges, surface_edges_verts, true);
2615  // get vertices on the body skin
2616  Range tets_skin_edges_verts;
2617  CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts,
2618  true);
2619 
2620  Range edges_to_remove;
2621 
2622  // remove edges self connected to body skin
2623  {
2624  Range ents_nodes_and_edges;
2625  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2626  ents_nodes_and_edges.merge(tets_skin_edges);
2627  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2628  0, false);
2629  }
2630  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2631  not_merged_edges.merge(edges_to_remove);
2632 
2633  // remove edges self connected to surface
2634  {
2635  Range ents_nodes_and_edges;
2636  ents_nodes_and_edges.merge(surface_edges_verts);
2637  ents_nodes_and_edges.merge(surface_edges);
2638  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2639  ents_nodes_and_edges.merge(tets_skin_edges);
2640  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2641  0, false);
2642  }
2643  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2644  not_merged_edges.merge(edges_to_remove);
2645 
2646  // remove edges adjacent corner_nodes execpt those on fixed edges
2647  Range fixed_edges_nodes;
2648  CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes, true);
2649  {
2650  Range ents_nodes_and_edges;
2651  ents_nodes_and_edges.merge(fixed_edges_nodes);
2652  ents_nodes_and_edges.merge(ends_nodes);
2653  ents_nodes_and_edges.merge(corner_nodes);
2654  ents_nodes_and_edges.merge(fixed_edges);
2655  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2656  0, false);
2657  }
2658  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2659  not_merged_edges.merge(edges_to_remove);
2660 
2661  // remove edges self connected to surface
2662  CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove, 0, false);
2663  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2664  not_merged_edges.merge(edges_to_remove);
2665 
2666  // remove edges self contented on surface skin
2667  {
2668  Range ents_nodes_and_edges;
2669  ents_nodes_and_edges.merge(surface_skin);
2670  ents_nodes_and_edges.merge(fixed_edges_nodes);
2671  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2672  0, false);
2673  }
2674  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2675  not_merged_edges.merge(edges_to_remove);
2676 
2677  // remove edges connecting crack front and fixed edges, except those short
2678  {
2679  Range ents_nodes_and_edges;
2680  ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
2681  ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2682  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2683  0, false);
2684  }
2685  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2686  not_merged_edges.merge(edges_to_remove);
2687 
2688  // remove crack front nodes connected to the surface, except those short
2689  {
2690  Range ents_nodes_and_edges;
2691  ents_nodes_and_edges.merge(surface_front_nodes);
2692  ents_nodes_and_edges.merge(surface_front);
2693  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2694  ents_nodes_and_edges.merge(tets_skin_edges);
2695  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2696  tOL, false);
2697  }
2698  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2699  not_merged_edges.merge(edges_to_remove);
2700 
2702  }
2703 
2704  private:
2705  MoFEMErrorCode removeSelfConectingEdges(const Range &ents,
2706  Range &edges_to_remove,
2707  const bool length,
2708  bool debug) const {
2709  moab::Interface &moab(mField.get_moab());
2711  // get nodes
2712  Range ents_nodes = ents.subset_by_type(MBVERTEX);
2713  if (ents_nodes.empty()) {
2714  CHKERR moab.get_connectivity(ents, ents_nodes, true);
2715  }
2716  // edges adj. to nodes
2717  Range ents_nodes_edges;
2718  CHKERR moab.get_adjacencies(ents_nodes, 1, false, ents_nodes_edges,
2719  moab::Interface::UNION);
2720  // nodes of adj. edges
2721  Range ents_nodes_edges_nodes;
2722  CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2723  true);
2724  // hanging nodes
2725  ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2726  Range ents_nodes_edges_nodes_edges;
2727  CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1, false,
2728  ents_nodes_edges_nodes_edges,
2729  moab::Interface::UNION);
2730  // remove edges adj. to hanging edges
2731  ents_nodes_edges =
2732  subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2733  ents_nodes_edges =
2734  subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2735  if (length > 0) {
2736  Range::iterator eit = ents_nodes_edges.begin();
2737  for (; eit != ents_nodes_edges.end();) {
2738 
2739  int num_nodes;
2740  const EntityHandle *conn;
2741  CHKERR moab.get_connectivity(*eit, conn, num_nodes, true);
2742  double coords[6];
2743  if (tH)
2744  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
2745  else
2746  CHKERR moab.get_coords(conn, num_nodes, coords);
2747 
2748  auto get_edge_length = [coords]() {
2750  &coords[0], &coords[1], &coords[2]);
2753  t_delta(i) = t_coords(i);
2754  ++t_coords;
2755  t_delta(i) -= t_coords(i);
2756  return sqrt(t_delta(i) * t_delta(i));
2757  };
2758 
2759  if (get_edge_length() < tOL) {
2760  eit = ents_nodes_edges.erase(eit);
2761  } else {
2762  eit++;
2763  }
2764  }
2765  }
2766  edges_to_remove.swap(ents_nodes_edges);
2767  if (debug) {
2768  CHKERR SaveData(moab)("edges_to_remove.vtk", edges_to_remove);
2769  }
2771  }
2772  };
2773 
2774  Range not_merged_edges;
2775  const double tol = 1e-3;
2776  CHKERR Toplogy(m_field, th, tol * aveLength)
2777  .removeBadEdges(surface, tets, fixed_edges, corner_nodes, edges_to_merge,
2778  not_merged_edges);
2779  Toplogy::SetsMap sets_map;
2780  CHKERR Toplogy(m_field, th, tol * aveLength)
2781  .classifyVerts(surface, tets, fixed_edges, corner_nodes, constrainSurface,
2782  sets_map);
2783  if (debug) {
2784  for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2785  sit != sets_map.rend(); sit++) {
2786  std::string name = "classification_verts_" +
2787  boost::lexical_cast<std::string>(sit->first) + ".vtk";
2788  if (!sit->second.empty())
2789  CHKERR SaveData(moab)(name, sit->second);
2790  }
2791  }
2792  Range proc_tets;
2793  CHKERR Toplogy(m_field, th, tol * aveLength)
2794  .getProcTets(tets, edges_to_merge, proc_tets);
2795  out_tets = subtract(tets, proc_tets);
2796 
2797  if (bit_ptr) {
2798  Range all_out_ents = out_tets;
2799  for (int dd = 2; dd >= 0; dd--) {
2800  CHKERR moab.get_adjacencies(out_tets, dd, false, all_out_ents,
2801  moab::Interface::UNION);
2802  }
2803  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(all_out_ents,
2804  *bit_ptr);
2805  }
2806 
2807  int nb_nodes_merged = 0;
2808  LengthMapData_multi_index length_map;
2809  new_surf = surface;
2810 
2811  auto save_merge_step = [&](const int pp, const Range collapsed_edges) {
2813  Range adj_faces;
2814  CHKERR moab.get_adjacencies(proc_tets, 2, false, adj_faces,
2815  moab::Interface::UNION);
2816  std::string name;
2817  name = "node_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2818  CHKERR SaveData(moab)(
2819  name, unite(intersect(new_surf, adj_faces), collapsed_edges));
2820  name =
2821  "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2822  CHKERR SaveData(moab)(
2823  name, unite(intersect(new_surf, adj_faces), edges_to_merge));
2825  };
2826 
2827  if (debug)
2828  CHKERR save_merge_step(0, Range());
2829 
2830  double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2831  for (int pp = 0; pp != nbMaxMergingCycles; ++pp) {
2832 
2833  int nb_nodes_merged_p = nb_nodes_merged;
2834  length_map.clear();
2835  min_pp = min_p;
2836  min_p = min;
2837  CHKERR LengthMap(m_field, th, aveLength)(proc_tets, edges_to_merge,
2838  length_map, ave);
2839  min = length_map.get<2>().begin()->qUality;
2840  if (pp == 0) {
2841  ave0 = ave;
2842  }
2843 
2844  int nn = 0;
2845  Range collapsed_edges;
2846  MergeNodes merge_nodes(m_field, true, th, update_meshsets);
2847 
2848  for (auto mit = length_map.get<0>().begin();
2849  mit != length_map.get<0>().end(); mit++, nn++) {
2850 
2851  if (!mit->skip) {
2852 
2853  auto get_father_and_mother =
2854  [&](EntityHandle &father, EntityHandle &mother, int &line_search) {
2856  int num_nodes;
2857  const EntityHandle *conn;
2858  CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes, true);
2859  std::array<int, 2> conn_type = {0, 0};
2860  for (int nn = 0; nn != 2; nn++) {
2861  for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2862  sit != sets_map.rend(); sit++) {
2863  if (sit->second.find(conn[nn]) != sit->second.end()) {
2864  conn_type[nn] |= sit->first;
2865  }
2866  }
2867  }
2868  int type_father, type_mother;
2869  if (conn_type[0] > conn_type[1]) {
2870  father = conn[0];
2871  mother = conn[1];
2872  type_father = conn_type[0];
2873  type_mother = conn_type[1];
2874  } else {
2875  father = conn[1];
2876  mother = conn[0];
2877  type_father = conn_type[1];
2878  type_mother = conn_type[0];
2879  }
2880  if (type_father == type_mother) {
2881  line_search = lineSearchSteps;
2882  }
2884  };
2885 
2886  int line_search = 0;
2887  EntityHandle father, mother;
2888  CHKERR get_father_and_mother(father, mother, line_search);
2889  CHKERR merge_nodes.mergeNodes(line_search, father, mother, proc_tets);
2890  if (m_field.getInterface<NodeMergerInterface>()->getSuccessMerge()) {
2891  const EntityHandle father_and_mother[] = {father, mother};
2892  Range adj_tets;
2893  CHKERR moab.get_adjacencies(father_and_mother, 1, 3, false, adj_tets);
2894  Range adj_tets_nodes;
2895  CHKERR moab.get_connectivity(adj_tets, adj_tets_nodes, true);
2896  Range adj_edges;
2897  CHKERR moab.get_adjacencies(adj_tets_nodes, 1, false, adj_edges,
2898  moab::Interface::UNION);
2899  for (auto ait : adj_edges) {
2900  auto miit = length_map.get<1>().find(ait);
2901  if (miit != length_map.get<1>().end())
2902  (const_cast<LengthMapData &>(*miit)).skip = true;
2903  }
2904  nb_nodes_merged++;
2905  collapsed_edges.insert(mit->eDge);
2906  }
2907 
2908  if (nn > static_cast<int>(length_map.size() / fraction_level))
2909  break;
2910  if (mit->qUality > ave)
2911  break;
2912  }
2913  }
2914 
2915  CHKERR merge_nodes.updateRangeByChilds(new_surf, edges_to_merge,
2916  not_merged_edges, true);
2917 
2918  PetscPrintf(m_field.get_comm(),
2919  "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e\n", pp,
2920  nb_nodes_merged, ave, min);
2921 
2922  if (debug)
2923  CHKERR save_merge_step(pp + 1, collapsed_edges);
2924 
2925  if (nb_nodes_merged == nb_nodes_merged_p)
2926  break;
2927  if (min > 1e-2 && min == min_pp)
2928  break;
2929  if (min > ave0)
2930  break;
2931 
2932  Range adj_edges;
2933  CHKERR moab.get_adjacencies(proc_tets, 1, false, adj_edges,
2934  moab::Interface::UNION);
2935  edges_to_merge = intersect(edges_to_merge, adj_edges);
2936  CHKERR Toplogy(m_field, th, tol * aveLength)
2937  .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
2938  edges_to_merge, not_merged_edges);
2939  }
2940 
2941  if (bit_ptr)
2942  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(proc_tets,
2943  *bit_ptr);
2944 
2945  out_tets.merge(proc_tets);
2946  Range adj_faces;
2947  CHKERR moab.get_adjacencies(out_tets, 2, false, adj_faces,
2948  moab::Interface::UNION);
2949  new_surf = intersect(new_surf, adj_faces);
2950 
2952 }
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:34
virtual moab::Interface & get_moab()=0
double maxLength
Maximal edge length.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
static constexpr double delta
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
double aveLength
Average edge length.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
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:108
#define CHKERR
Inline error check.
Definition: definitions.h:601
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26

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

2958  {
2959  CoreInterface &m_field = cOre;
2961  Range tets_level;
2962  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2963  trim_bit, BitRefLevel().set(), MBTET, tets_level);
2964 
2965  Range edges_to_merge;
2966  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
2967  trim_bit, cut_bit | trim_bit, edges_to_merge);
2968 
2969  // get all entities not in database
2970  Range all_ents_not_in_database_before;
2971  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2972  all_ents_not_in_database_before);
2973 
2974  edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
2975  if (debug)
2976  CHKERR SaveData(m_field.get_moab())("edges_to_merge.vtk", edges_to_merge);
2977 
2978  Range out_new_tets, out_new_surf;
2979  CHKERR mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
2980  corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
2981  th, update_meshsets, &bit, debug);
2982 
2983  // get all entities not in database after merge
2984  Range all_ents_not_in_database_after;
2985  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2986  all_ents_not_in_database_after);
2987 
2988  // delete hanging entities
2989  all_ents_not_in_database_after =
2990  subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
2991  Range meshsets;
2992  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
2993  false);
2994  for (auto m : meshsets)
2995  CHKERR m_field.get_moab().remove_entities(m,
2996  all_ents_not_in_database_after);
2997 
2998  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
2999 
3000  mergedVolumes.swap(out_new_tets);
3001  mergedSurfaces.swap(out_new_surf);
3003 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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.
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:601
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ mergeSurface()

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

merge surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 131 of file CutMeshInterface.cpp.

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

◆ mergeVolumes()

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

merge volume entities

Parameters
volumeentities which going to be added
Returns
error code

Definition at line 137 of file CutMeshInterface.cpp.

137  {
139  vOlume.merge(volume);
141 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513

◆ moveMidNodesOnCutEdges()

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

projecting of mid edge nodes on new mesh on surface

Returns
error code

Definition at line 1447 of file CutMeshInterface.cpp.

1447  {
1449 
1450  CoreInterface &m_field = cOre;
1451  moab::Interface &moab = m_field.get_moab();
1453 
1454  // Range out_side_vertices;
1455  for (auto m : verticesOnCutEdges) {
1456  double dist = m.second.dIst;
1457  VectorDouble3 new_coors = m.second.rayPoint + dist * m.second.unitRayDir;
1458  if (th)
1459  CHKERR moab.tag_set_data(th, &m.first, 1, &new_coors[0]);
1460  else
1461  CHKERR moab.set_coords(&m.first, 1, &new_coors[0]);
1462  }
1463 
1465 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
#define CHKERR
Inline error check.
Definition: definitions.h:601
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
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:412

◆ moveMidNodesOnTrimmedEdges()

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

move trimmed edges mid nodes

Returns
error code

Definition at line 1467 of file CutMeshInterface.cpp.

1467  {
1468  CoreInterface &m_field = cOre;
1469  moab::Interface &moab = m_field.get_moab();
1471  for (auto &v : verticesOnTrimEdges) {
1472  double dist = v.second.dIst;
1473  VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1474  if (th)
1475  CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1476  else
1477  CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1478  }
1480 }
virtual moab::Interface & get_moab()=0
map< EntityHandle, TreeData > verticesOnTrimEdges
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define CHKERR
Inline error check.
Definition: definitions.h:601
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ projectZeroDistanceEnts() [1/2]

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

Find entities on cut surface which can be projected.

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

Definition at line 968 of file CutMeshInterface.cpp.

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

◆ projectZeroDistanceEnts() [2/2]

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

Definition at line 488 of file CutMeshInterface.hpp.

488  {
489  return zeroDistanceEnts;
490  }

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

785  {
786  CoreInterface &m_field = cOre;
787  moab::Interface &moab = m_field.get_moab();
788  MeshRefinement *refiner;
789  BitRefManager *bit_ref_manager;
791  CHKERR m_field.getInterface(refiner);
792  CHKERR m_field.getInterface(bit_ref_manager);
793 
794  auto add_bit = [&](const int bit) {
796  CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
797  verb);
798  Range adj_ents;
799  for (auto d : {2, 1, 0})
800  CHKERR moab.get_adjacencies(vOlume, d, true, adj_ents,
801  moab::Interface::UNION);
802  CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
803  verb);
805  };
806  CHKERR add_bit(init_bit_level);
807 
808  auto update_range = [&](Range *r_ptr) {
810  if (r_ptr) {
811  Range childs;
812  CHKERR bit_ref_manager->updateRange(*r_ptr, childs);
813  r_ptr->merge(childs);
814  }
816  };
817 
818  auto refine = [&](const BitRefLevel bit, const Range tets) {
820  Range verts;
821  CHKERR moab.get_connectivity(tets, verts, true);
822  Range ref_edges;
823  CHKERR moab.get_adjacencies(verts, 1, true, ref_edges,
824  moab::Interface::UNION);
825 
826  CHKERR refiner->add_vertices_in_the_middle_of_edges(ref_edges, bit);
827  CHKERR refiner->refine_TET(vOlume, bit, false, verb);
828 
829  CHKERR update_range(fixed_edges);
830  CHKERR update_range(&vOlume);
831  CHKERR m_field.getInterface<MeshsetsManager>()
832  ->updateAllMeshsetsByEntitiesChildren(bit);
833 
834  Range bit_tets;
835  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
836  bit, BitRefLevel().set(), MBTET, bit_tets);
837  vOlume = intersect(vOlume, bit_tets);
838 
839  Range bit_edges;
840  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
841  bit, BitRefLevel().set(), MBEDGE, bit_edges);
842  if (fixed_edges)
843  *fixed_edges = intersect(*fixed_edges, bit_edges);
844 
846  };
847 
848  for (int ll = init_bit_level; ll != init_bit_level + surf_levels; ++ll) {
850  CHKERR refine(BitRefLevel().set(ll + 1),
852  }
853 
854  for (int ll = init_bit_level + surf_levels;
855  ll != init_bit_level + surf_levels + front_levels; ++ll) {
857  CHKERR refine(BitRefLevel().set(ll + 1), cutFrontVolumes);
858  }
859 
861 
862  if (debug)
863  CHKERR SaveData(m_field.get_moab())("refinedVolume.vtk", vOlume);
864 
866 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:601
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
MoFEMErrorCode findLevelSetVolumes(Tag th, Range &vol_edges, const bool remove_adj_prims_edges, int verb=QUIET, const bool debug=false, const std::string edges_file_name=string())
Find level set volumes.

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

2139  {
2140  CoreInterface &m_field = cOre;
2141  moab::Interface &moab = m_field.get_moab();
2142  PrismInterface *interface;
2144  CHKERR m_field.getInterface(interface);
2145  // Remove tris on surface front
2146  {
2147  Range front_tris;
2148  EntityHandle meshset;
2149  CHKERR moab.create_meshset(MESHSET_SET, meshset);
2150  CHKERR moab.add_entities(meshset, ents);
2151  CHKERR interface->findFacesWithThreeNodesOnInternalSurfaceSkin(
2152  meshset, split_bit, true, front_tris);
2153  CHKERR moab.delete_entities(&meshset, 1);
2154  ents = subtract(ents, front_tris);
2155  }
2156  // Remove entities on skin
2157  Range tets;
2158  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2159  split_bit, BitRefLevel().set(), MBTET, tets);
2160  // Remove entities on skin
2161  Skinner skin(&moab);
2162  Range tets_skin;
2163  CHKERR skin.find_skin(0, tets, false, tets_skin);
2164  ents = subtract(ents, tets_skin);
2165 
2167 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define CHKERR
Inline error check.
Definition: definitions.h:601
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ saveCutEdges()

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

Definition at line 3038 of file CutMeshInterface.cpp.

3038  {
3039  CoreInterface &m_field = cOre;
3040  moab::Interface &moab = m_field.get_moab();
3042  CHKERR SaveData(moab)(prefix + "out_vol.vtk", vOlume);
3043  CHKERR SaveData(moab)(prefix + "out_surface.vtk", sUrface);
3044  CHKERR SaveData(moab)(prefix + "out_cut_edges.vtk", cutEdges);
3045  CHKERR SaveData(moab)(prefix + "out_cut_new_volumes.vtk", cutNewVolumes);
3046  CHKERR SaveData(moab)(prefix + "out_cut_new_surfaces.vtk", cutNewSurfaces);
3047  CHKERR SaveData(moab)(prefix + "out_cut_zero_distance_ents.vtk",
3050 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define CHKERR
Inline error check.
Definition: definitions.h:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ saveTrimEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveTrimEdges ( )

Definition at line 3052 of file CutMeshInterface.cpp.

3052  {
3053  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
3055  CHKERR SaveData(moab)("out_trim_surface.vtk", sUrface);
3056  CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
3057  CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
3058  CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
3060 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ setConstrainSurface()

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

Set the constrain surface object.

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

Parameters
surf
Returns
MoFEMErrorCode

Definition at line 125 of file CutMeshInterface.cpp.

125  {
127  constrainSurface = surf;
129 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513

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

3022  {
3023  CoreInterface &m_field = cOre;
3024  moab::Interface &moab = m_field.get_moab();
3026  Range nodes;
3027  if (bit.none())
3028  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3029  else
3030  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3031  bit, mask, MBVERTEX, nodes);
3032  std::vector<double> coords(3 * nodes.size());
3033  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
3034  CHKERR moab.set_coords(nodes, &coords[0]);
3036 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define CHKERR
Inline error check.
Definition: definitions.h:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

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

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

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

3005  {
3006  CoreInterface &m_field = cOre;
3007  moab::Interface &moab = m_field.get_moab();
3009  Range nodes;
3010  if (bit.none())
3011  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3012  else
3013  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3014  bit, BitRefLevel().set(), MBVERTEX, nodes);
3015  std::vector<double> coords(3 * nodes.size());
3016  CHKERR moab.get_coords(nodes, &coords[0]);
3017  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
3019 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define CHKERR
Inline error check.
Definition: definitions.h:601
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ setTrimFixedEdges()

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

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

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

147  {
148  CoreInterface &m_field = cOre;
149  moab::Interface &moab = m_field.get_moab();
151 
152  // Get cutting surface skin
153  Skinner skin(&moab);
154  Range surface_skin;
155  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
156 
157  CHKERR snapSurfaceToEdges(surface_skin, fixed_edges, rel_tol, abs_tol, th,
158  debug);
159 
161 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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)
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

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

167  {
168  CoreInterface &m_field = cOre;
169  moab::Interface &moab = m_field.get_moab();
172 
173  map<EntityHandle, double> map_verts_length;
174 
175  for (auto f : surface_edges) {
176  int num_nodes;
177  const EntityHandle *conn_skin;
178  CHKERR moab.get_connectivity(f, conn_skin, num_nodes, true);
179  VectorDouble6 coords_skin(6);
180  if (th)
181  CHKERR moab.tag_get_data(th, conn_skin, num_nodes, &coords_skin[0]);
182  else
183  CHKERR moab.get_coords(conn_skin, num_nodes, &coords_skin[0]);
185  &coords_skin[0], &coords_skin[1], &coords_skin[2]);
187  &coords_skin[3], &coords_skin[4], &coords_skin[5]);
188  FTensor::Tensor1<double, 3> t_skin_delta;
189  t_skin_delta(i) = t_n1(i) - t_n0(i);
190  const double skin_edge_length = sqrt(t_skin_delta(i) * t_skin_delta(i));
191  for (int nn = 0; nn != 2; ++nn) {
192  auto it = map_verts_length.find(conn_skin[nn]);
193  if (it != map_verts_length.end())
194  it->second = std::min(it->second, skin_edge_length);
195  else
196  map_verts_length[conn_skin[nn]] = skin_edge_length;
197  }
198  }
199 
200  for (auto m : map_verts_length) {
201 
203  if (th)
204  CHKERR moab.tag_get_data(th, &m.first, 1, &t_n(0));
205  else
206  CHKERR moab.get_coords(&m.first, 1, &t_n(0));
207 
208  double min_dist = rel_tol * m.second;
209  FTensor::Tensor1<double, 3> t_min_coords;
210  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
211  &t_n(0), 1, fixed_edges, &min_dist, &t_min_coords(0));
212 
213  if (min_dist < rel_tol * m.second || min_dist < abs_tol) {
214  if (debug)
215  cerr << "Snap " << min_dist << endl;
216  if (th)
217  CHKERR moab.tag_set_data(th, &m.first, 1, &t_min_coords(0));
218  else
219  CHKERR moab.set_coords(&m.first, 1, &t_min_coords(0));
220  }
221  }
222 
224 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:94
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26

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

2171  {
2172  CoreInterface &m_field = cOre;
2173  moab::Interface &moab = m_field.get_moab();
2174  PrismInterface *interface;
2176  CHKERR m_field.getInterface(interface);
2177  EntityHandle meshset_volume;
2178  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
2179  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2180  split_bit, BitRefLevel().set(), MBTET, meshset_volume);
2181  EntityHandle meshset_surface;
2182  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
2183  CHKERR moab.add_entities(meshset_surface, ents);
2184  CHKERR interface->getSides(meshset_surface, split_bit, true);
2185  CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
2186  true);
2187  CHKERR moab.delete_entities(&meshset_volume, 1);
2188  CHKERR moab.delete_entities(&meshset_surface, 1);
2189  if (th) {
2190  Range prisms;
2191  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2192  bit, BitRefLevel().set(), MBPRISM, prisms);
2193  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
2194  int num_nodes;
2195  const EntityHandle *conn;
2196  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
2197  MatrixDouble data(3, 3);
2198  CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
2199  CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
2200  }
2201  }
2203 }
virtual moab::Interface & get_moab()=0
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define CHKERR
Inline error check.
Definition: definitions.h:601
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

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

1949  {
1950  CoreInterface &m_field = cOre;
1951  moab::Interface &moab = m_field.get_moab();
1952  MeshRefinement *refiner;
1953  const RefEntity_multiIndex *ref_ents_ptr;
1955 
1956  CHKERR m_field.getInterface(refiner);
1957  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1958  CHKERR refiner->add_vertices_in_the_middle_of_edges(trimEdges, bit);
1959  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
1960  debug ? &trimEdges : NULL);
1961 
1962  trimNewVolumes.clear();
1963  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1964  bit, BitRefLevel().set(), MBTET, trimNewVolumes);
1965 
1966  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
1967  mit != edgesToTrim.end(); mit++) {
1968  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1969  boost::make_tuple(MBVERTEX, mit->first));
1970  if (vit ==
1971  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end())
1972  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1973  "No vertex on trim edges, that make no sense");
1974 
1975  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1976  if ((ref_ent->getBitRefLevel() & bit).any()) {
1977  EntityHandle vert = ref_ent->getRefEnt();
1978  trimNewVertices.insert(vert);
1979  verticesOnTrimEdges[vert] = mit->second;
1980  }
1981  }
1982 
1983  // Get faces which are trimmed
1984  trimNewSurfaces.clear();
1985  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1986  bit, bit, MBTRI, trimNewSurfaces);
1987 
1988  Range trim_new_surfaces_nodes;
1989  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
1990  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
1991  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
1992  Range faces_not_on_surface;
1993  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
1994  faces_not_on_surface, moab::Interface::UNION);
1995  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
1996 
1997  // Get surfaces which are not trimmed and add them to surface
1998  Range all_surfaces_on_bit_level;
1999  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2000  bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
2001  all_surfaces_on_bit_level =
2002  intersect(all_surfaces_on_bit_level, cutNewSurfaces);
2003 
2004  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
2005 
2007 }
map< EntityHandle, TreeData > edgesToTrim
virtual moab::Interface & get_moab()=0
map< EntityHandle, TreeData > verticesOnTrimEdges
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:601
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
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
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

◆ trimOnly()

MoFEMErrorCode MoFEM::CutMeshInterface::trimOnly ( const BitRefLevel  trim_bit,
Tag  th,
const double  tol_geometry,
const double  tol_trim_close,
Range *  fixed_edges = NULL,
Range *  corner_nodes = NULL,
const bool  update_meshsets = false,
const bool  debug = false 
)

Trim mesh only.

Parameters
trim_bit
th
tol_cut_close
fixed_edges
corner_nodes
update_meshsets
debug
Returns
MoFEMErrorCode

Definition at line 274 of file CutMeshInterface.cpp.

280  {
281  CoreInterface &m_field = cOre;
282  moab::Interface &moab = m_field.get_moab();
284 
285  // trim mesh
286  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_geometry,
287  tol_trim_close, debug);
288  CHKERR trimEdgesInTheMiddle(trim_bit, debug);
289 
290  CHKERR cOre.getInterface<BitRefManager>()->updateRange(constrainSurface,
292  if (fixed_edges)
293  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
294  *fixed_edges);
295 
296  if (corner_nodes)
297  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
298  *corner_nodes);
299 
300  if (update_meshsets)
301  CHKERR m_field.getInterface<MeshsetsManager>()
302  ->updateAllMeshsetsByEntitiesChildren(trim_bit);
303 
304  // move nodes
306 
307  // remove faces
308  CHKERR trimSurface(fixed_edges, corner_nodes, debug);
309 
310  if (debug) {
312  Range bit2_edges;
313  CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
314  trim_bit, BitRefLevel().set(), MBEDGE, bit2_edges);
315  CHKERR SaveData(moab)("trim_fixed_edges.vtk",
316  intersect(*fixed_edges, bit2_edges));
317  }
318 
320 }
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
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:601
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol_geometry=1e-4, const double tol=1e-4, const bool debug=false)
Find edges to trimEdges.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
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 2009 of file CutMeshInterface.cpp.

2011  {
2012  CoreInterface &m_field = cOre;
2013  moab::Interface &moab = m_field.get_moab();
2015 
2016  Skinner skin(&moab);
2017  Range trim_tets_skin;
2018  CHKERR skin.find_skin(0, trimNewVolumes, false, trim_tets_skin);
2019  Range trim_tets_skin_edges;
2020  CHKERR moab.get_adjacencies(trim_tets_skin, 1, false, trim_tets_skin_edges,
2021  moab::Interface::UNION);
2022 
2023  Range barrier_vertices(trimNewVertices);
2024 
2025  if (fixed_edges && trimFixedEdges) {
2026 
2027  // get all vertices on fixed edges and surface
2028  Range trim_surface_edges;
2029  CHKERR moab.get_adjacencies(trimNewSurfaces, 1, false, trim_surface_edges,
2030  moab::Interface::UNION);
2031  Range fixed_edges_on_trim_surface;
2032  fixed_edges_on_trim_surface = intersect(trim_surface_edges, *fixed_edges);
2033  Range fixed_edges_on_trim_surface_verts;
2034  CHKERR moab.get_connectivity(fixed_edges_on_trim_surface,
2035  fixed_edges_on_trim_surface_verts, false);
2036 
2037  // get faces adjacent to barrier_vertices
2038  Range barrier_vertices_faces;
2039  CHKERR moab.get_adjacencies(barrier_vertices, 2, false,
2040  barrier_vertices_faces, moab::Interface::UNION);
2041  barrier_vertices_faces = intersect(barrier_vertices_faces, trimNewSurfaces);
2042 
2043  // get vertices on fixed edges
2044  Range fixed_edges_vertices;
2045  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_vertices, false);
2046  fixed_edges_vertices = intersect(barrier_vertices, fixed_edges_vertices);
2047  fixed_edges_vertices =
2048  subtract(fixed_edges_vertices, fixed_edges_on_trim_surface_verts);
2049  if (corner_nodes)
2050  fixed_edges_vertices.merge(intersect(barrier_vertices, *corner_nodes));
2051 
2052  // get faces adjacent to vertices on fixed edges
2053  Range fixed_edges_faces;
2054  CHKERR moab.get_adjacencies(fixed_edges_vertices, 2, false,
2055  fixed_edges_faces, moab::Interface::UNION);
2056  fixed_edges_faces = intersect(fixed_edges_faces, barrier_vertices_faces);
2057 
2058  if (debug && !fixed_edges_faces.empty())
2059  CHKERR SaveData(m_field.get_moab())("fixed_edges_faces.vtk",
2060  fixed_edges_faces);
2061 
2062  // get nodes on faces
2063  Range fixed_edges_faces_vertices;
2064  CHKERR moab.get_connectivity(fixed_edges_faces, fixed_edges_faces_vertices,
2065  false);
2066  barrier_vertices.merge(fixed_edges_faces_vertices);
2067  }
2068 
2069  auto remove_faces_on_skin = [&]() {
2071  Range skin_faces = intersect(trimNewSurfaces, trim_tets_skin);
2072  if (!skin_faces.empty()) {
2073  Range add_to_barrier;
2074  CHKERR moab.get_connectivity(skin_faces, add_to_barrier, false);
2075  barrier_vertices.merge(add_to_barrier);
2076  for (auto f : skin_faces)
2077  trimNewSurfaces.erase(f);
2078  }
2080  };
2081 
2082  auto get_trim_free_edges = [&]() {
2083  // get current surface skin
2084  Range trim_surf_skin;
2085  CHKERR skin.find_skin(0, trimNewSurfaces, false, trim_surf_skin);
2086  trim_surf_skin = subtract(trim_surf_skin, trim_tets_skin_edges);
2087 
2088  Range trim_surf_skin_verts;
2089  CHKERR moab.get_connectivity(trim_surf_skin, trim_surf_skin_verts, true);
2090  for (auto e : barrier_vertices)
2091  trim_surf_skin_verts.erase(e);
2092 
2093  Range free_edges;
2094  CHKERR moab.get_adjacencies(trim_surf_skin_verts, 1, false, free_edges,
2095  moab::Interface::UNION);
2096  free_edges = intersect(free_edges, trim_surf_skin);
2097 
2098  return free_edges;
2099  };
2100 
2101  CHKERR remove_faces_on_skin();
2102 
2103  if (debug && !barrier_vertices.empty())
2104  CHKERR SaveData(m_field.get_moab())("barrier_vertices.vtk",
2105  barrier_vertices);
2106 
2107  int nn = 0;
2108  Range out_edges;
2109  while (!(out_edges = get_trim_free_edges()).empty()) {
2110 
2111  if (debug && !trimNewSurfaces.empty())
2112  CHKERR SaveData(m_field.get_moab())(
2113  "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
2114  trimNewSurfaces);
2115 
2116  if (debug && !out_edges.empty())
2117  CHKERR SaveData(m_field.get_moab())(
2118  "trimNewSurfacesOutsideVerts_" +
2119  boost::lexical_cast<std::string>(nn) + ".vtk",
2120  out_edges);
2121 
2122  Range outside_faces;
2123  CHKERR moab.get_adjacencies(out_edges, 2, false, outside_faces,
2124  moab::Interface::UNION);
2125  trimNewSurfaces = subtract(trimNewSurfaces, outside_faces);
2126  ++nn;
2127  }
2128 
2129  if (debug && !trimNewSurfaces.empty())
2130  CHKERR SaveData(m_field.get_moab())(
2131  "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
2132  trimNewSurfaces);
2133 
2135 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:601
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412

Member Data Documentation

◆ aveLength

double MoFEM::CutMeshInterface::aveLength
private

Average edge length.

Definition at line 591 of file CutMeshInterface.hpp.

◆ constrainSurface

Range MoFEM::CutMeshInterface::constrainSurface
private

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

◆ cutFrontVolumes

Range MoFEM::CutMeshInterface::cutFrontVolumes
private

Definition at line 595 of file CutMeshInterface.hpp.

◆ cutNewSurfaces

Range MoFEM::CutMeshInterface::cutNewSurfaces
private

Definition at line 566 of file CutMeshInterface.hpp.

◆ cutNewVertices

Range MoFEM::CutMeshInterface::cutNewVertices
private

Definition at line 569 of file CutMeshInterface.hpp.

◆ cutNewVolumes

Range MoFEM::CutMeshInterface::cutNewVolumes
private

Definition at line 565 of file CutMeshInterface.hpp.

◆ cutSurfaceVolumes

Range MoFEM::CutMeshInterface::cutSurfaceVolumes
private

Definition at line 594 of file CutMeshInterface.hpp.

◆ edgesToCut

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

Definition at line 586 of file CutMeshInterface.hpp.

◆ edgesToTrim

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

Definition at line 588 of file CutMeshInterface.hpp.

◆ fRont

Range MoFEM::CutMeshInterface::fRont
private

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

◆ mergedSurfaces

Range MoFEM::CutMeshInterface::mergedSurfaces
private

Definition at line 577 of file CutMeshInterface.hpp.

◆ mergedVolumes

Range MoFEM::CutMeshInterface::mergedVolumes
private

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

◆ sUrface

Range MoFEM::CutMeshInterface::sUrface
private

Definition at line 556 of file CutMeshInterface.hpp.

◆ treeSurfPtr

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

Definition at line 561 of file CutMeshInterface.hpp.

◆ trimEdges

Range MoFEM::CutMeshInterface::trimEdges
private

Definition at line 574 of file CutMeshInterface.hpp.

◆ trimFixedEdges

bool MoFEM::CutMeshInterface::trimFixedEdges
private

Definition at line 559 of file CutMeshInterface.hpp.

◆ trimNewSurfaces

Range MoFEM::CutMeshInterface::trimNewSurfaces
private

Definition at line 573 of file CutMeshInterface.hpp.

◆ trimNewVertices

Range MoFEM::CutMeshInterface::trimNewVertices
private

Definition at line 572 of file CutMeshInterface.hpp.

◆ trimNewVolumes

Range MoFEM::CutMeshInterface::trimNewVolumes
private

Definition at line 571 of file CutMeshInterface.hpp.

◆ verticesOnCutEdges

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

Definition at line 587 of file CutMeshInterface.hpp.

◆ verticesOnTrimEdges

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

Definition at line 589 of file CutMeshInterface.hpp.

◆ vOlume

Range MoFEM::CutMeshInterface::vOlume
private

Definition at line 557 of file CutMeshInterface.hpp.

◆ zeroDistanceEnts

Range MoFEM::CutMeshInterface::zeroDistanceEnts
private

Definition at line 567 of file CutMeshInterface.hpp.

◆ zeroDistanceVerts

Range MoFEM::CutMeshInterface::zeroDistanceVerts
private

Definition at line 568 of file CutMeshInterface.hpp.


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