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_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
 Trim mesh only. More...
 
MoFEMErrorCode cutAndTrim (int &first_bit, Tag th, const double tol_geometry, const double tol_cut_close, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
 Cut and trim. More...
 
MoFEMErrorCode cutTrimAndMerge (int &first_bit, const int fraction_level, Tag th, const double tol_geometry, const double tol_cut_close, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)
 Cut, trim and merge. More...
 
MoFEMErrorCode makeFront (const bool debug=false)
 Create front from the surface. More...
 
MoFEMErrorCode createSurfaceLevelSets (int verb=QUIET, const bool debug=false)
 Calculate distance from mesh nodes to cut surface. More...
 
MoFEMErrorCode createFrontLevelSets (Range vol, Tag th=nullptr, int verb=QUIET, const bool debug=false)
 Calculate distance from mesh nodes to surface front. More...
 
MoFEMErrorCode findLevelSetVolumes (Tag th, Range &vol_edges, const bool remove_adj_prims_edges, int verb=QUIET, const bool debug=false, const std::string edges_file_name=string())
 Find level set volumes. More...
 
MoFEMErrorCode findLevelSetVolumes (int verb=QUIET, const bool debug=false)
 Find level set volumes. More...
 
MoFEMErrorCode refineMesh (const int init_bit_level, const int surf_levels, const int front_levels, Range *fixed_edges=nullptr, int verb=QUIET, const bool debug=false)
 Refine and set level sets. More...
 
MoFEMErrorCode findEdgesToCut (Range vol, int verb=QUIET, const bool debug=false)
 find edges to cut More...
 
MoFEMErrorCode projectZeroDistanceEnts (Range *fixed_edges, Range *corner_nodes, const double geometry_tol=0, const double close_tol=0, const int verb=QUIET, const bool debug=false)
 Find entities on cut surface which can be projected. More...
 
MoFEMErrorCode cutEdgesInMiddle (const BitRefLevel bit, Range &cut_vols, Range &cut_surf, Range &cut_verts, const bool debug=false)
 cut edges More...
 
MoFEMErrorCode moveMidNodesOnCutEdges (Tag th=NULL)
 projecting of mid edge nodes on new mesh on surface More...
 
MoFEMErrorCode findEdgesToTrim (Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, const bool debug=false)
 Find edges to trimEdges. More...
 
MoFEMErrorCode trimEdgesInTheMiddle (const BitRefLevel bit, const bool debug=false)
 trim edges More...
 
MoFEMErrorCode moveMidNodesOnTrimmedEdges (Tag th=NULL)
 move trimmed edges mid nodes More...
 
MoFEMErrorCode trimSurface (Range *fixed_edge, Range *corner_nodes, const double tol=1e-4, Tag th=NULL, const bool debug=false)
 Trim surface from faces beyond front. More...
 
MoFEMErrorCode removePathologicalFrontTris (const BitRefLevel split_bit, Range &ents)
 Remove pathological elements on surface internal front. More...
 
MoFEMErrorCode splitSides (const BitRefLevel split_bit, const BitRefLevel bit, const Range &ents, Tag th=NULL)
 split sides More...
 
MoFEMErrorCode mergeBadEdges (const int fraction_level, const Range &tets, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Range &merged_nodes, Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets=false, const BitRefLevel *bit_ptr=NULL, const bool debug=false)
 Merge edges. More...
 
MoFEMErrorCode mergeBadEdges (const int fraction_level, const BitRefLevel cut_bit, const BitRefLevel trim_bit, const BitRefLevel bit, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Tag th=NULL, const bool update_meshsets=false, const bool debug=false)
 Merge edges. More...
 
MoFEMErrorCode setTagData (Tag th, const BitRefLevel bit=BitRefLevel())
 set coords to tag More...
 
MoFEMErrorCode setCoords (Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
 set coords from tag More...
 
const 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
 
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
 
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 549 of file CutMeshInterface.hpp.

Constructor & Destructor Documentation

◆ CutMeshInterface()

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

Definition at line 39 of file CutMeshInterface.cpp.

40  : cOre(const_cast<Core &>(core)) {
41  lineSearchSteps = 10;
42  nbMaxMergingCycles = 200;
44 }

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

231  {
232  CoreInterface &m_field = cOre;
233  moab::Interface &moab = m_field.get_moab();
235  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
236  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
239 }
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:483
#define CHKERR
Inline error check.
Definition: definitions.h:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413

◆ clearMap()

MoFEMErrorCode MoFEM::CutMeshInterface::clearMap ( )

Definition at line 46 of file CutMeshInterface.cpp.

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

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

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

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

618  {
619  CoreInterface &m_field = cOre;
620  moab::Interface &moab = m_field.get_moab();
622 
623  auto create_tag = [&](const std::string name, const int dim) {
624  Tag th;
625  rval = moab.tag_get_handle(name.c_str(), th);
626  if (rval == MB_SUCCESS)
627  return th;
628  std::vector<double> def_val(dim, 0);
629  CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
630  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
631  return th;
632  };
633 
634  Range vol_vertices;
635  CHKERR moab.get_connectivity(vol, vol_vertices, true);
636 
637  std::vector<double> min_distances_from_front(vol_vertices.size(), -1);
638  std::vector<double> points_on_edges(3 * vol_vertices.size(), 0);
639  std::vector<EntityHandle> closest_edges(vol_vertices.size(), 0);
640 
641  std::vector<double> coords(3 * vol_vertices.size());
642  if (th)
643  CHKERR moab.tag_get_data(th, vol_vertices, &*coords.begin());
644  else
645  CHKERR moab.get_coords(vol_vertices, &*coords.begin());
646 
647  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
648  &*coords.begin(), vol_vertices.size(), fRont,
649  &*min_distances_from_front.begin(), &*points_on_edges.begin(),
650  &*closest_edges.begin());
651 
652  if (!points_on_edges.empty()) {
653  for (int i = 0; i != min_distances_from_front.size(); ++i) {
654  Range faces;
655  CHKERR moab.get_adjacencies(&closest_edges[0], 1, 2, false, faces);
656  auto point_in = getVectorAdaptor(&coords[3 * i], 3);
657  auto point_out = getVectorAdaptor(&points_on_edges[3 * i], 3);
658  point_out -= point_in;
659  }
660  }
661 
662  auto th_dist_front_vec = create_tag("DIST_FRONT_VECTOR", 3);
663  CHKERR moab.tag_set_data(th_dist_front_vec, vol_vertices,
664  &*points_on_edges.begin());
665 
667 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
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:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
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 484 of file CutMeshInterface.cpp.

485  {
486  CoreInterface &m_field = cOre;
487  moab::Interface &moab = m_field.get_moab();
489  auto tools_interface = m_field.getInterface<Tools>();
490 
491  auto create_tag = [&](const std::string name, const int dim) {
492  Tag th;
493  rval = moab.tag_get_handle(name.c_str(), th);
494  if (rval == MB_SUCCESS)
495  return th;
496  std::vector<double> def_val(dim, 0);
497  CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
498  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
499 
500  return th;
501  };
502 
503  auto set_vol = [&](const Range &vol_verts, std::vector<double> &coords,
504  std::vector<double> &dist_surface_vec,
505  std::vector<double> &dist_surface_normal_vec,
506  std::vector<double> &dist_surface_signed_dist_vec) {
508 
509  coords.resize(3 * vol_verts.size());
510  dist_surface_vec.resize(3 * vol_verts.size());
511  dist_surface_normal_vec.resize(3 * vol_verts.size());
512  dist_surface_signed_dist_vec.resize(vol_verts.size());
513 
514  CHKERR moab.get_coords(vol_verts, &*coords.begin());
515  std::srand(0);
516 
517  for (auto v : vol_verts) {
518 
519  const int index = vol_verts.index(v);
520  auto point_in = getVectorAdaptor(&coords[3 * index], 3);
521  VectorDouble3 point_out(3);
522  EntityHandle facets_out;
523  CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
524  &point_out[0], facets_out);
525 
526  VectorDouble3 delta = point_out - point_in;
527  VectorDouble3 n_first(3);
528  CHKERR tools_interface->getTriNormal(facets_out, &*n_first.begin());
529 
530  // Check of only triangle in proximity of point out is one triangle, if is
531  // more than one, use normal of lager triangle to set orientation.
532  auto check_triangle_orientation = [&](auto n) {
533  int num_nodes;
534  const EntityHandle *conn;
535  CHKERR moab.get_connectivity(facets_out, conn, num_nodes, true);
536  MatrixDouble3by3 coords(3, 3);
537  CHKERR moab.get_coords(conn, 3, &coords(0, 0));
538  VectorDouble3 center(3);
539  center.clear();
540  for (auto ii : {0, 1, 2})
541  for (auto jj : {0, 1, 2})
542  center(jj) += coords(ii, jj) / 3;
543 
544  std::vector<EntityHandle> triangles_out;
545  std::vector<double> distance_out;
546  auto a_max = norm_2(n);
547  VectorDouble3 ray = -n / a_max;
548  VectorDouble3 pt = center - ray * sqrt(a_max);
549  const double ray_length = 2 * sqrt(a_max);
550 
551  CHKERR treeSurfPtr->ray_intersect_triangles(
552  distance_out, triangles_out, rootSetSurf,
553  std::numeric_limits<float>::epsilon(), &pt[0], &ray[0],
554  &ray_length);
555 
556  if (triangles_out.size() > 1) {
557  VectorDouble3 nt(3);
558  for (auto t : triangles_out) {
559  CHKERR tools_interface->getTriNormal(t, &*nt.begin());
560  double at = norm_2(nt);
561  if (at > a_max) {
562  n = nt;
563  a_max = at;
564  }
565  }
566  }
567 
568  return n;
569  };
570 
571  auto n = check_triangle_orientation(n_first);
572  n /= norm_2(n);
573  const double dot = inner_prod(delta, n);
574 
575 
576  if (std::abs(dot) < std::numeric_limits<double>::epsilon()) {
577  if (std::rand() % 2 == 0)
578  delta += n * std::numeric_limits<double>::epsilon();
579  else
580  delta -= n * std::numeric_limits<double>::epsilon();
581  }
582 
583  auto dist_vec = getVectorAdaptor(&dist_surface_vec[3 * index], 3);
584  noalias(dist_vec) = delta;
585 
586  auto dist_normal_vec =
587  getVectorAdaptor(&dist_surface_normal_vec[3 * index], 3);
588 
589  dist_surface_signed_dist_vec[index] = dot;
590  noalias(dist_normal_vec) = dot * n;
591  }
592 
594  };
595 
596  std::vector<double> coords;
597  std::vector<double> dist_surface_vec;
598  std::vector<double> dist_surface_normal_vec;
599  std::vector<double> dist_surface_signed_dist_vec;
600  Range vol_verts;
601  CHKERR moab.get_connectivity(vOlume, vol_verts, true);
602 
603  CHKERR set_vol(vol_verts, coords, dist_surface_vec, dist_surface_normal_vec,
604  dist_surface_signed_dist_vec);
605 
606  CHKERR moab.tag_set_data(create_tag("DIST_SURFACE_VECTOR", 3), vol_verts,
607  &*dist_surface_vec.begin());
608  CHKERR moab.tag_set_data(create_tag("DIST_SURFACE_NORMAL_VECTOR", 3),
609  vol_verts, &*dist_surface_normal_vec.begin());
610  CHKERR moab.tag_set_data(create_tag("DIST_SURFACE_NORMAL_SIGNED", 1),
611  vol_verts, &*dist_surface_signed_dist_vec.begin());
612 
614 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
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:483
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
const int dim
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:42
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:101
#define CHKERR
Inline error check.
Definition: definitions.h:602
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413

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

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

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

1320  {
1321  CoreInterface &m_field = cOre;
1322  moab::Interface &moab = m_field.get_moab();
1323  MeshRefinement *refiner;
1324  auto ref_ents_ptr = m_field.get_ref_ents();
1326 
1327  if (cutEdges.size() != edgesToCut.size())
1328  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1329 
1330  auto refine_mesh = [&]() {
1332  CHKERR m_field.getInterface(refiner);
1333  CHKERR refiner->add_vertices_in_the_middle_of_edges(cutEdges, bit);
1334  CHKERR refiner->refine_TET(vOlume, bit, false, QUIET,
1335  debug ? &cutEdges : NULL);
1337  };
1338 
1339  CHKERR refine_mesh();
1340 
1341  if (debug)
1342  if (cutEdges.size() != edgesToCut.size())
1343  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1344 
1345  cut_vols.clear();
1346  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1347  bit, BitRefLevel().set(), MBTET, cut_vols);
1348  cut_surf.clear();
1349  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1350  bit, bit, MBTRI, cut_surf);
1351 
1352  if (debug)
1353  CHKERR SaveData(moab)("cut_surf_from_bit.vtk", cut_surf);
1354 
1355  // Find new vertices on cut edges
1356  cut_verts.clear();
1357  cut_verts.merge(zeroDistanceVerts);
1358 
1359  for (auto &m : edgesToCut) {
1360  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1361  boost::make_tuple(MBVERTEX, m.first));
1362  if (vit ==
1363  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1364  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1365  "No vertex on cut edges, that make no sense");
1366  }
1367  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1368  if ((ref_ent->getBitRefLevel() & bit).any()) {
1369  EntityHandle vert = ref_ent->getRefEnt();
1370  cut_verts.insert(vert);
1371  verticesOnCutEdges[vert] = m.second;
1372  } else {
1373  SETERRQ1(
1374  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1375  "Vertex has wrong bit ref level %s",
1376  boost::lexical_cast<std::string>(ref_ent->getBitRefLevel()).c_str());
1377  }
1378  }
1379 
1380  // Add zero distance entities faces
1381  Range tets_skin;
1382  Skinner skin(&moab);
1383  CHKERR skin.find_skin(0, cut_vols, false, tets_skin);
1384  tets_skin.merge(constrainSurface);
1385 
1386  // At that point cut_surf has all newly created faces, now take all
1387  // nodes on those faces and subtract nodes on cut edges. Faces adjacent to
1388  // nodes which left are not part of surface.
1389  cut_surf.merge(zeroDistanceEnts.subset_by_type(MBTRI));
1390  Range diff_verts;
1391  CHKERR moab.get_connectivity(cut_surf, diff_verts, true);
1392  diff_verts = subtract(diff_verts, cut_verts);
1393 
1394  Range subtract_faces;
1395  CHKERR moab.get_adjacencies(diff_verts, 2, false, subtract_faces,
1396  moab::Interface::UNION);
1397  cut_surf = subtract(cut_surf, unite(subtract_faces, tets_skin));
1398 
1399  cut_verts.clear();
1400  CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1401 
1402  // Check non-mainfolds
1403  auto check_for_non_minfold = [&]() {
1405  Range surf_edges;
1406  CHKERR moab.get_adjacencies(cut_surf, 1, false, surf_edges,
1407  moab::Interface::UNION);
1408  for (auto e : surf_edges) {
1409 
1410  Range faces;
1411  CHKERR moab.get_adjacencies(&e, 1, 2, false, faces);
1412  faces = intersect(faces, cut_surf);
1413  if (faces.size() > 2) {
1414 
1415  bool resolved = false;
1416 
1417  // Check for haning node
1418  Range nodes;
1419  CHKERR moab.get_connectivity(faces, nodes, true);
1420  for (auto n : nodes) {
1421  Range adj_faces;
1422  CHKERR moab.get_adjacencies(&n, 1, 2, false, adj_faces);
1423  adj_faces = intersect(adj_faces, cut_surf);
1424  if (adj_faces.size() == 1) {
1425  cut_surf.erase(adj_faces[0]);
1426  resolved = true;
1427  }
1428  }
1429 
1430  // Check for two edges minfold
1431  Range adj_edges;
1432  CHKERR moab.get_adjacencies(faces, 1, false, adj_edges,
1433  moab::Interface::UNION);
1434  adj_edges = intersect(adj_edges, surf_edges);
1435  adj_edges.erase(e);
1436  for (auto other_e : adj_edges) {
1437  Range other_faces;
1438  CHKERR moab.get_adjacencies(&other_e, 1, 2, false, other_faces);
1439  other_faces = intersect(other_faces, cut_surf);
1440  if (other_faces.size() > 2) {
1441  other_faces = intersect(other_faces, faces);
1442  cut_surf = subtract(cut_surf, other_faces);
1443  resolved = true;
1444  }
1445  }
1446 
1447  if (!resolved && !debug)
1448  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1449  "Non-mainfold surface");
1450 
1451  cut_verts.clear();
1452  CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1453  }
1454  }
1456  };
1457 
1458  CHKERR check_for_non_minfold();
1459 
1460  if (debug)
1461  CHKERR SaveData(moab)("cut_surf.vtk", cut_surf);
1462 
1464 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
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:483
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:602
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:1879
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:413
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 242 of file CutMeshInterface.cpp.

245  {
246  CoreInterface &m_field = cOre;
247  moab::Interface &moab = m_field.get_moab();
249 
250  // cut mesh
252  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_geometry,
253  tol_cut_close, QUIET, debug);
256 
257  CHKERR cOre.getInterface<BitRefManager>()->updateRange(constrainSurface,
259  if (fixed_edges)
260  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
261  *fixed_edges);
262  if (corner_nodes)
263  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
264  *corner_nodes);
265  if (update_meshsets)
266  CHKERR m_field.getInterface<MeshsetsManager>()
267  ->updateAllMeshsetsByEntitiesChildren(cut_bit);
269 
270  if (debug) {
272  if (fixed_edges)
273  CHKERR SaveData(moab)("out_cut_new_fixed_edges.vtk", *fixed_edges);
274  }
275 
277 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879

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

398  {
399  CoreInterface &m_field = cOre;
401 
402  std::vector<BitRefLevel> bit_levels;
403 
404  auto get_back_bit_levels = [&]() {
405  bit_levels.push_back(BitRefLevel());
406  bit_levels.back().set(first_bit + bit_levels.size() - 1);
407  return bit_levels.back();
408  };
409 
410  if (debug) {
411  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
412  "ents_not_in_database.vtk", "VTK", "");
413  }
414 
415  CHKERR cutAndTrim(first_bit, th, tol_geometry, tol_cut_close, tol_trim_close,
416  &fixed_edges, &corner_nodes, update_meshsets, debug);
417  if (debug)
418  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
419  "cut_trim_ents_not_in_database.vtk", "VTK", "");
420 
421  BitRefLevel bit_level1 = BitRefLevel().set(first_bit - 1);
422  BitRefLevel bit_level2 = get_back_bit_levels();
423  BitRefLevel bit_level3 = get_back_bit_levels();
424 
425  CHKERR mergeBadEdges(fraction_level, bit_level1, bit_level2, bit_level3,
426  getNewTrimSurfaces(), fixed_edges, corner_nodes, th,
427  update_meshsets, debug);
428 
429  auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
430  Range tets_level; // test at level
431  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
432  bit, BitRefLevel().set(), MBTET, tets_level);
433  double min_q = 1;
434  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
435  if (min_q < 0 && debug) {
436  CHKERR m_field.getInterface<Tools>()->writeTetsWithQuality(
437  "negative_tets.vtk", "VTK", "", tets_level, th);
438  }
439  return min_q;
440  };
441 
442  MOFEM_LOG_C("WORLD", Sev::inform, "Min quality node merge %6.4g",
443  get_min_quality(bit_level3, th));
444 
445  CHKERR cOre.getInterface<BitRefManager>()->updateRange(constrainSurface,
447  CHKERR cOre.getInterface<BitRefManager>()->updateRange(fixed_edges,
448  fixed_edges);
449  CHKERR cOre.getInterface<BitRefManager>()->updateRange(corner_nodes,
450  corner_nodes);
451 
452  first_bit += bit_levels.size() - 1;
453 
454  if (debug) {
455  CHKERR cOre.getInterface<BitRefManager>()->writeBitLevelByType(
456  bit_level3, BitRefLevel().set(), MBTET, "out_tets_merged.vtk", "VTK",
457  "");
458  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
459  "cut_trim_merge_ents_not_in_database.vtk", "VTK", "");
460  CHKERR SaveData(m_field.get_moab())("merged_surface.vtk", mergedSurfaces);
461  }
462 
464 }
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:295
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:602
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
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 903 of file CutMeshInterface.cpp.

904  {
905  CoreInterface &m_field = cOre;
906  moab::Interface &moab = m_field.get_moab();
908 
909  edgesToCut.clear();
910  cutEdges.clear();
911 
912  Tag th_signed_dist;
913  CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_SIGNED", th_signed_dist);
914 
915  auto get_tag_edge_dist = [&](auto th, auto conn) {
916  std::array<double, 2> r;
917  CHKERR moab.tag_get_data(th, conn, 2, r.data());
918  return r;
919  };
920 
921  auto get_conn = [&](const auto e) {
922  int num_nodes;
923  const EntityHandle *conn;
924  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
925  return conn;
926  };
927 
928  auto get_adj = [&moab](const Range r, const int dim) {
929  Range a;
930  if (dim)
931  CHKERR moab.get_adjacencies(r, dim, false, a, moab::Interface::UNION);
932  else
933  CHKERR moab.get_connectivity(r, a, true);
934  return a;
935  };
936 
937  auto get_range = [](std::vector<EntityHandle> v) {
938  Range r;
939  r.insert_list(v.begin(), v.end());
940  return r;
941  };
942 
943  auto vol_edges = get_adj(vol, 1);
944 
945  aveLength = 0;
946  maxLength = 0;
947  int nb_ave_length = 0;
948  for (auto e : vol_edges) {
949 
950  auto conn = get_conn(e);
951  auto dist = get_tag_edge_dist(th_signed_dist, conn);
952  const auto dist_max = std::max(dist[0], dist[1]);
953  const auto dist_min = std::min(dist[0], dist[1]);
954  const auto dot = dist[0] * dist[1];
955 
956  VectorDouble6 coords(6);
957  CHKERR moab.get_coords(conn, 2, &coords[0]);
958  VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
959  VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
960  VectorDouble3 ray = n1 - n0;
961  const double ray_length = norm_2(ray);
962  ray /= ray_length;
963 
964  if (
965 
966  (dot < 0 && dist_max > std::numeric_limits<double>::epsilon()) ||
967  (std::abs(dist_min) < std::numeric_limits<double>::epsilon() &&
968  dist_max > std::numeric_limits<double>::epsilon())
969 
970  ) {
971 
972  // Edges is on two sides of the surface
973  const double s = dist[0] / (dist[0] - dist[1]);
974  const double dist = s * ray_length;
975 
976  auto add_edge = [&](auto dist) {
977  edgesToCut[e].dIst = dist;
978  edgesToCut[e].lEngth = ray_length;
979  edgesToCut[e].unitRayDir = ray;
980  edgesToCut[e].rayPoint = n0;
981  cutEdges.insert(e);
982 
983  aveLength += norm_2(ray);
984  maxLength = fmax(maxLength, norm_2(ray));
985  ++nb_ave_length;
986  };
987 
988  add_edge(dist);
989  }
990  }
991  aveLength /= nb_ave_length;
992 
993  if (debug)
994  CHKERR SaveData(m_field.get_moab())("cut_edges.vtk", cutEdges);
995 
996  if (debug)
997  CHKERR SaveData(m_field.get_moab())("cut_edges_vol.vtk",
998  get_adj(cutEdges, 3));
999 
1001 }
virtual moab::Interface & get_moab()=0
double maxLength
Maximal edge length.
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
map< EntityHandle, TreeData > edgesToCut

◆ findEdgesToTrim()

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

Find edges to trimEdges.

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

Parameters
verbverbosity level
Returns
error code

Definition at line 1499 of file CutMeshInterface.cpp.

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

671  {
672  CoreInterface &m_field = cOre;
673  moab::Interface &moab = m_field.get_moab();
675 
676  auto get_tag_dim = [&](auto th) {
677  int dim;
678  CHKERR moab.tag_get_length(th, dim);
679  return dim;
680  };
681  auto dim = get_tag_dim(th);
682 
683  auto get_tag_data = [&](auto th, auto conn) {
684  const void *ptr;
685  CHKERR moab.tag_get_by_ptr(th, &conn, 1, &ptr);
686  return getVectorAdaptor(
687  const_cast<double *>(static_cast<const double *>(ptr)), 3);
688  };
689 
690  auto get_edge_ray = [&](auto conn) {
691  VectorDouble6 coords(6);
692  CHKERR moab.get_coords(conn, 2, &coords[0]);
693  VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
694  VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
695  VectorDouble3 ray = n1 - n0;
696  return ray;
697  };
698 
699  Range edges;
700  CHKERR moab.get_adjacencies(vOlume, 1, true, edges, moab::Interface::UNION);
701 
702  auto remove_prisms_edges = [&](Range &edges) {
704  Range prisms;
705  CHKERR moab.get_adjacencies(edges, 3, false, prisms,
706  moab::Interface::UNION);
707  prisms = prisms.subset_by_type(MBPRISM);
708  Range prisms_verts;
709  CHKERR moab.get_connectivity(prisms, prisms_verts, true);
710  Range prism_edges;
711  CHKERR moab.get_adjacencies(prisms_verts, 1, false, prism_edges,
712  moab::Interface::UNION);
713  edges = subtract(edges, prism_edges);
715  };
716  if (remove_adj_prims_edges)
717  CHKERR remove_prisms_edges(edges);
718 
719  std::vector<EntityHandle> cut_edges;
720  cut_edges.reserve(edges.size());
721 
722  auto get_cut_edges_vec = [&](auto th, auto &cut_edges, auto e, auto &&conn) {
724 
725  auto ray = get_edge_ray(conn);
726  const double length = norm_2(ray);
727  ray /= length;
728  const auto dist0 = get_tag_data(th, conn[0]);
729  const auto dist1 = get_tag_data(th, conn[1]);
730  const double max_dist = std::max(norm_2(dist0), norm_2(dist1));
731  if (max_dist < length) {
732  cut_edges.push_back(e);
733  }
734 
736  };
737 
738  auto get_cut_edges_signed_dist = [&](auto th, auto &cut_edges, auto e,
739  auto &&conn) {
741  auto get_tag_signed_dist = [&](auto conn) {
742  double dist;
743  CHKERR moab.tag_get_data(th, &conn, 1, &dist);
744  return dist;
745  };
746  const auto dist0 = get_tag_signed_dist(conn[0]);
747  const auto dist1 = get_tag_signed_dist(conn[1]);
748  const double opposite = dist0 * dist1;
749  if (opposite <= 0)
750  cut_edges.push_back(e);
752  };
753 
754  auto get_conn = [&](auto e) {
755  int num_nodes;
756  const EntityHandle *conn;
757  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
758  return conn;
759  };
760 
761  if (dim == 3)
762  for (auto e : edges)
763  CHKERR get_cut_edges_vec(th, cut_edges, e, get_conn(e));
764  else
765  for (auto e : edges)
766  CHKERR get_cut_edges_signed_dist(th, cut_edges, e, get_conn(e));
767 
768  CHKERR moab.get_adjacencies(&*cut_edges.begin(), cut_edges.size(), 3, false,
769  vol_edges, moab::Interface::UNION);
770  vol_edges = intersect(vol_edges, vOlume);
771 
772  auto convert_to_ranger = [](auto &v) {
773  Range r;
774  r.insert_list(v.begin(), v.end());
775  return r;
776  };
777 
778  if (debug && !edges_file_name.empty())
779  CHKERR SaveData(m_field.get_moab())(edges_file_name,
780  convert_to_ranger(cut_edges));
781 
783 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413

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

786  {
787  CoreInterface &m_field = cOre;
788  moab::Interface &moab = m_field.get_moab();
790 
791  CHKERR createFrontLevelSets(vOlume, nullptr, verb, debug);
792  Tag th_dist_front_vec;
793  CHKERR moab.tag_get_handle("DIST_FRONT_VECTOR", th_dist_front_vec);
794  CHKERR findLevelSetVolumes(th_dist_front_vec, cutFrontVolumes, true, verb,
795  debug, "cutFrontEdges.vtk");
796 
798 
799  Tag th_dist_surface_vec;
800  CHKERR moab.tag_get_handle("DIST_SURFACE_VECTOR", th_dist_surface_vec);
801  cutSurfaceVolumes.clear();
802  CHKERR findLevelSetVolumes(th_dist_surface_vec, cutSurfaceVolumes, true, verb,
803  debug, "cutSurfaceEdges.vtk");
804 
805  if (debug)
806  CHKERR SaveData(m_field.get_moab())("level_sets.vtk", vOlume);
807  if (debug)
808  CHKERR SaveData(m_field.get_moab())("cutSurfaceVolumes.vtk",
810  if (debug)
811  CHKERR SaveData(m_field.get_moab())("cutFrontVolumes.vtk", cutFrontVolumes);
812 
814 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
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 483 of file CutMeshInterface.hpp.

483 { return cutEdges; }

◆ getCutFrontVolumes()

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

Definition at line 500 of file CutMeshInterface.hpp.

◆ getCutSurfaceVolumes()

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

Definition at line 499 of file CutMeshInterface.hpp.

◆ getFront()

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

Definition at line 481 of file CutMeshInterface.hpp.

481 { return fRont; }

◆ getMergedSurfaces()

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

Definition at line 497 of file CutMeshInterface.hpp.

◆ getMergedVolumes()

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

Definition at line 496 of file CutMeshInterface.hpp.

496 { return mergedVolumes; }

◆ getNewCutSurfaces()

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

Definition at line 485 of file CutMeshInterface.hpp.

◆ getNewCutVertices()

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

Definition at line 486 of file CutMeshInterface.hpp.

◆ getNewCutVolumes()

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

Definition at line 484 of file CutMeshInterface.hpp.

484 { return cutNewVolumes; }

◆ getNewTrimSurfaces()

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

Definition at line 493 of file CutMeshInterface.hpp.

◆ getNewTrimVertices()

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

Definition at line 494 of file CutMeshInterface.hpp.

◆ getNewTrimVolumes()

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

Definition at line 492 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:483
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:550
#define CHKERR
Inline error check.
Definition: definitions.h:602
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413

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

480 { return sUrface; }

◆ getTreeSurfPtr()

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

Definition at line 506 of file CutMeshInterface.hpp.

506  {
507  return treeSurfPtr;
508  }
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr

◆ getTrimEdges()

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

Definition at line 491 of file CutMeshInterface.hpp.

491 { return trimEdges; }

◆ getVolume()

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

Definition at line 479 of file CutMeshInterface.hpp.

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

466  {
467  CoreInterface &m_field = cOre;
468  moab::Interface &moab = m_field.get_moab();
470  Skinner skin(&moab);
471  Range tets_skin;
472  CHKERR skin.find_skin(0, vOlume, false, tets_skin);
473  Range tets_skin_edges;
474  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
475  moab::Interface::UNION);
476  Range surface_skin;
477  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
478  fRont = subtract(surface_skin, tets_skin_edges);
479  if (debug)
480  CHKERR SaveData(moab)("front.vtk", fRont);
482 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879

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

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

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

3083  {
3084  CoreInterface &m_field = cOre;
3086  Range tets_level;
3087  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3088  trim_bit, BitRefLevel().set(), MBTET, tets_level);
3089 
3090  Range edges_to_merge;
3091  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
3092  trim_bit, cut_bit | trim_bit, edges_to_merge);
3093 
3094  // get all entities not in database
3095  Range all_ents_not_in_database_before;
3096  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
3097  all_ents_not_in_database_before);
3098 
3099  edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
3100  if (debug)
3101  CHKERR SaveData(m_field.get_moab())("edges_to_merge.vtk", edges_to_merge);
3102 
3103  Range out_new_tets, out_new_surf;
3104  CHKERR mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
3105  corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
3106  th, update_meshsets, &bit, debug);
3107 
3108  // get all entities not in database after merge
3109  Range all_ents_not_in_database_after;
3110  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
3111  all_ents_not_in_database_after);
3112 
3113  // delete hanging entities
3114  all_ents_not_in_database_after =
3115  subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
3116  Range meshsets;
3117  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
3118  false);
3119  for (auto m : meshsets)
3120  CHKERR m_field.get_moab().remove_entities(m,
3121  all_ents_not_in_database_after);
3122 
3123  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
3124 
3125  mergedVolumes.swap(out_new_tets);
3126  mergedSurfaces.swap(out_new_surf);
3127 
3129 }
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50

◆ mergeSurface()

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

merge surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 136 of file CutMeshInterface.cpp.

136  {
138  sUrface.merge(surface);
140 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514

◆ mergeVolumes()

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

merge volume entities

Parameters
volumeentities which going to be added
Returns
error code

Definition at line 142 of file CutMeshInterface.cpp.

142  {
144  vOlume.merge(volume);
146 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514

◆ moveMidNodesOnCutEdges()

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

projecting of mid edge nodes on new mesh on surface

Returns
error code

Definition at line 1466 of file CutMeshInterface.cpp.

1466  {
1467  CoreInterface &m_field = cOre;
1468  moab::Interface &moab = m_field.get_moab();
1470 
1471  // Range out_side_vertices;
1472  for (auto m : verticesOnCutEdges) {
1473  double dist = m.second.dIst;
1474  VectorDouble3 new_coors = m.second.rayPoint + dist * m.second.unitRayDir;
1475  if (th)
1476  CHKERR moab.tag_set_data(th, &m.first, 1, &new_coors[0]);
1477  else
1478  CHKERR moab.set_coords(&m.first, 1, &new_coors[0]);
1479  }
1480 
1482 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
#define CHKERR
Inline error check.
Definition: definitions.h:602
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
map< EntityHandle, TreeData > verticesOnCutEdges

◆ moveMidNodesOnTrimmedEdges()

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

move trimmed edges mid nodes

Returns
error code

Definition at line 1484 of file CutMeshInterface.cpp.

1484  {
1485  CoreInterface &m_field = cOre;
1486  moab::Interface &moab = m_field.get_moab();
1488  for (auto &v : verticesOnTrimEdges) {
1489  double dist = v.second.dIst;
1490  VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1491  if (th)
1492  CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1493  else
1494  CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1495  }
1497 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
map< EntityHandle, TreeData > verticesOnTrimEdges
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
#define CHKERR
Inline error check.
Definition: definitions.h:602
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879

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

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

487  {
488  return zeroDistanceEnts;
489  }

◆ query_interface()

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

Implements MoFEM::UnknownInterface.

Definition at line 27 of file CutMeshInterface.cpp.

28  {
30  *iface = NULL;
31  if (uuid == IDD_MOFEMCutMesh) {
32  *iface = const_cast<CutMeshInterface *>(this);
34  }
35  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
37 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
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 816 of file CutMeshInterface.cpp.

820  {
821  CoreInterface &m_field = cOre;
822  moab::Interface &moab = m_field.get_moab();
823  MeshRefinement *refiner;
824  BitRefManager *bit_ref_manager;
826  CHKERR m_field.getInterface(refiner);
827  CHKERR m_field.getInterface(bit_ref_manager);
828 
829  auto add_bit = [&](const int bit) {
831  CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
832  verb);
833  Range adj_ents;
834  for (auto d : {2, 1, 0})
835  CHKERR moab.get_adjacencies(vOlume, d, true, adj_ents,
836  moab::Interface::UNION);
837  CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
838  verb);
840  };
841  CHKERR add_bit(init_bit_level);
842 
843  auto update_range = [&](Range *r_ptr) {
845  if (r_ptr) {
846  Range childs;
847  CHKERR bit_ref_manager->updateRange(*r_ptr, childs);
848  r_ptr->merge(childs);
849  }
851  };
852 
853  auto refine = [&](const BitRefLevel bit, const Range tets) {
855  Range verts;
856  CHKERR moab.get_connectivity(tets, verts, true);
857  Range ref_edges;
858  CHKERR moab.get_adjacencies(verts, 1, true, ref_edges,
859  moab::Interface::UNION);
860 
861  CHKERR refiner->add_vertices_in_the_middle_of_edges(ref_edges, bit);
862  CHKERR refiner->refine_TET(vOlume, bit, false, verb);
863 
864  CHKERR update_range(fixed_edges);
865  CHKERR update_range(&vOlume);
866  CHKERR m_field.getInterface<MeshsetsManager>()
867  ->updateAllMeshsetsByEntitiesChildren(bit);
868 
869  Range bit_tets;
870  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
871  bit, BitRefLevel().set(), MBTET, bit_tets);
872  vOlume = intersect(vOlume, bit_tets);
873 
874  Range bit_edges;
875  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
876  bit, BitRefLevel().set(), MBEDGE, bit_edges);
877  if (fixed_edges)
878  *fixed_edges = intersect(*fixed_edges, bit_edges);
879 
881  };
882 
883  for (int ll = init_bit_level; ll != init_bit_level + surf_levels; ++ll) {
885  CHKERR refine(BitRefLevel().set(ll + 1),
887  }
888 
889  for (int ll = init_bit_level + surf_levels;
890  ll != init_bit_level + surf_levels + front_levels; ++ll) {
892  CHKERR refine(BitRefLevel().set(ll + 1), cutFrontVolumes);
893  }
894 
896 
897  if (debug)
898  CHKERR SaveData(m_field.get_moab())("refinedVolume.vtk", vOlume);
899 
901 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
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:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
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 2254 of file CutMeshInterface.cpp.

2255  {
2256  CoreInterface &m_field = cOre;
2257  moab::Interface &moab = m_field.get_moab();
2258  PrismInterface *interface;
2260  CHKERR m_field.getInterface(interface);
2261  // Remove tris on surface front
2262  {
2263  Range front_tris;
2264  EntityHandle meshset;
2265  CHKERR moab.create_meshset(MESHSET_SET, meshset);
2266  CHKERR moab.add_entities(meshset, ents);
2267  CHKERR interface->findFacesWithThreeNodesOnInternalSurfaceSkin(
2268  meshset, split_bit, true, front_tris);
2269  CHKERR moab.delete_entities(&meshset, 1);
2270  ents = subtract(ents, front_tris);
2271  }
2272  // Remove entities on skin
2273  Range tets;
2274  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2275  split_bit, BitRefLevel().set(), MBTET, tets);
2276  // Remove entities on skin
2277  Skinner skin(&moab);
2278  Range tets_skin;
2279  CHKERR skin.find_skin(0, tets, false, tets_skin);
2280  ents = subtract(ents, tets_skin);
2281 
2283 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
#define CHKERR
Inline error check.
Definition: definitions.h:602
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:1879

◆ saveCutEdges()

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

Definition at line 3164 of file CutMeshInterface.cpp.

3164  {
3165  CoreInterface &m_field = cOre;
3166  moab::Interface &moab = m_field.get_moab();
3168  CHKERR SaveData(moab)(prefix + "out_vol.vtk", vOlume);
3169  CHKERR SaveData(moab)(prefix + "out_surface.vtk", sUrface);
3170  CHKERR SaveData(moab)(prefix + "out_cut_edges.vtk", cutEdges);
3171  CHKERR SaveData(moab)(prefix + "out_cut_new_volumes.vtk", cutNewVolumes);
3172  CHKERR SaveData(moab)(prefix + "out_cut_new_surfaces.vtk", cutNewSurfaces);
3173  CHKERR SaveData(moab)(prefix + "out_cut_zero_distance_ents.vtk",
3176 }
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:483
#define CHKERR
Inline error check.
Definition: definitions.h:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413

◆ saveTrimEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveTrimEdges ( )

Definition at line 3178 of file CutMeshInterface.cpp.

3178  {
3179  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
3181  CHKERR SaveData(moab)("out_trim_surface.vtk", sUrface);
3182  CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
3183  CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
3184  CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
3186 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413

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

130  {
132  constrainSurface = surf;
134 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514

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

3148  {
3149  CoreInterface &m_field = cOre;
3150  moab::Interface &moab = m_field.get_moab();
3152  Range nodes;
3153  if (bit.none())
3154  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3155  else
3156  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3157  bit, mask, MBVERTEX, nodes);
3158  std::vector<double> coords(3 * nodes.size());
3159  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
3160  CHKERR moab.set_coords(nodes, &coords[0]);
3162 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
#define CHKERR
Inline error check.
Definition: definitions.h:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879

◆ setFront()

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

Definition at line 52 of file CutMeshInterface.cpp.

52  {
54  fRont = front;
56 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514

◆ setSurface()

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

set surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 58 of file CutMeshInterface.cpp.

58  {
60  sUrface = surface;
62 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514

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

3131  {
3132  CoreInterface &m_field = cOre;
3133  moab::Interface &moab = m_field.get_moab();
3135  Range nodes;
3136  if (bit.none())
3137  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3138  else
3139  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3140  bit, BitRefLevel().set(), MBVERTEX, nodes);
3141  std::vector<double> coords(3 * nodes.size());
3142  CHKERR moab.get_coords(nodes, &coords[0]);
3143  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
3145 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
#define CHKERR
Inline error check.
Definition: definitions.h:602
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:1879

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

124  {
126  vOlume = volume;
128 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514

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

152  {
153  CoreInterface &m_field = cOre;
154  moab::Interface &moab = m_field.get_moab();
156 
157  // Get cutting surface skin
158  Skinner skin(&moab);
159  Range surface_skin;
160  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
161 
162  CHKERR snapSurfaceToEdges(surface_skin, fixed_edges, rel_tol, abs_tol, th,
163  debug);
164 
166 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879

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

172  {
173  CoreInterface &m_field = cOre;
174  moab::Interface &moab = m_field.get_moab();
177 
178  map<EntityHandle, double> map_verts_length;
179 
180  for (auto f : surface_edges) {
181  int num_nodes;
182  const EntityHandle *conn_skin;
183  CHKERR moab.get_connectivity(f, conn_skin, num_nodes, true);
184  VectorDouble6 coords_skin(6);
185  if (th)
186  CHKERR moab.tag_get_data(th, conn_skin, num_nodes, &coords_skin[0]);
187  else
188  CHKERR moab.get_coords(conn_skin, num_nodes, &coords_skin[0]);
190  &coords_skin[0], &coords_skin[1], &coords_skin[2]);
192  &coords_skin[3], &coords_skin[4], &coords_skin[5]);
193  FTensor::Tensor1<double, 3> t_skin_delta;
194  t_skin_delta(i) = t_n1(i) - t_n0(i);
195  const double skin_edge_length = sqrt(t_skin_delta(i) * t_skin_delta(i));
196  for (int nn = 0; nn != 2; ++nn) {
197  auto it = map_verts_length.find(conn_skin[nn]);
198  if (it != map_verts_length.end())
199  it->second = std::min(it->second, skin_edge_length);
200  else
201  map_verts_length[conn_skin[nn]] = skin_edge_length;
202  }
203  }
204 
205  for (auto m : map_verts_length) {
206 
208  if (th)
209  CHKERR moab.tag_get_data(th, &m.first, 1, &t_n(0));
210  else
211  CHKERR moab.get_coords(&m.first, 1, &t_n(0));
212 
213  double min_dist = rel_tol * m.second;
214  FTensor::Tensor1<double, 3> t_min_coords;
215  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
216  &t_n(0), 1, fixed_edges, &min_dist, &t_min_coords(0));
217 
218  if (min_dist < rel_tol * m.second || min_dist < abs_tol) {
219  if (debug)
220  MOFEM_LOG("WORLD", Sev::noisy) << "Snap " << min_dist;
221  if (th)
222  CHKERR moab.tag_set_data(th, &m.first, 1, &t_min_coords(0));
223  else
224  CHKERR moab.set_coords(&m.first, 1, &t_min_coords(0));
225  }
226  }
227 
229 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:292
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
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 2285 of file CutMeshInterface.cpp.

2287  {
2288  CoreInterface &m_field = cOre;
2289  moab::Interface &moab = m_field.get_moab();
2290  PrismInterface *interface;
2292  CHKERR m_field.getInterface(interface);
2293  EntityHandle meshset_volume;
2294  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
2295  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2296  split_bit, BitRefLevel().set(), MBTET, meshset_volume);
2297  EntityHandle meshset_surface;
2298  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
2299  CHKERR moab.add_entities(meshset_surface, ents);
2300  CHKERR interface->getSides(meshset_surface, split_bit, true);
2301  CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
2302  true);
2303  CHKERR moab.delete_entities(&meshset_volume, 1);
2304  CHKERR moab.delete_entities(&meshset_surface, 1);
2305  if (th) {
2306  Range prisms;
2307  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2308  bit, BitRefLevel().set(), MBPRISM, prisms);
2309  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
2310  int num_nodes;
2311  const EntityHandle *conn;
2312  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
2313  MatrixDouble data(3, 3);
2314  CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
2315  CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
2316  }
2317  }
2319 }
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
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:483
#define CHKERR
Inline error check.
Definition: definitions.h:602
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:1879

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

1988  {
1989  CoreInterface &m_field = cOre;
1990  moab::Interface &moab = m_field.get_moab();
1991  MeshRefinement *refiner;
1992  auto ref_ents_ptr = m_field.get_ref_ents();
1994 
1995  CHKERR m_field.getInterface(refiner);
1996  CHKERR refiner->add_vertices_in_the_middle_of_edges(trimEdges, bit);
1997  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
1998  debug ? &trimEdges : NULL);
1999 
2000  trimNewVolumes.clear();
2001  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2002  bit, BitRefLevel().set(), MBTET, trimNewVolumes);
2003 
2004  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
2005  mit != edgesToTrim.end(); mit++) {
2006  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
2007  boost::make_tuple(MBVERTEX, mit->first));
2008  if (vit ==
2009  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end())
2010  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2011  "No vertex on trim edges, that make no sense");
2012 
2013  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
2014  if ((ref_ent->getBitRefLevel() & bit).any()) {
2015  EntityHandle vert = ref_ent->getRefEnt();
2016  trimNewVertices.insert(vert);
2017  verticesOnTrimEdges[vert] = mit->second;
2018  }
2019  }
2020 
2021  // Get faces which are trimmed
2022  trimNewSurfaces.clear();
2023  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2024  bit, bit, MBTRI, trimNewSurfaces);
2025 
2026  Range trim_new_surfaces_nodes;
2027  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
2028  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
2029  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
2030  Range faces_not_on_surface;
2031  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
2032  faces_not_on_surface, moab::Interface::UNION);
2033  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
2034 
2035  // Get surfaces which are not trimmed and add them to surface
2036  Range all_surfaces_on_bit_level;
2037  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2038  bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
2039  all_surfaces_on_bit_level =
2040  intersect(all_surfaces_on_bit_level, cutNewSurfaces);
2041 
2042  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
2043 
2045 }
map< EntityHandle, TreeData > edgesToTrim
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
map< EntityHandle, TreeData > verticesOnTrimEdges
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:602
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:1879

◆ trimOnly()

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

Trim mesh only.

Parameters
trim_bit
th
tol_cut_close
fixed_edges
corner_nodes
update_meshsets
debug
Returns
MoFEMErrorCode

Definition at line 279 of file CutMeshInterface.cpp.

284  {
285  CoreInterface &m_field = cOre;
286  moab::Interface &moab = m_field.get_moab();
288 
289  // trim mesh
290  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim_close, debug);
291  CHKERR trimEdgesInTheMiddle(trim_bit, debug);
292 
293  CHKERR cOre.getInterface<BitRefManager>()->updateRange(constrainSurface,
295  if (fixed_edges)
296  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
297  *fixed_edges);
298 
299  if (corner_nodes)
300  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
301  *corner_nodes);
302 
303  if (update_meshsets)
304  CHKERR m_field.getInterface<MeshsetsManager>()
305  ->updateAllMeshsetsByEntitiesChildren(trim_bit);
306 
307  // move nodes
309 
310  // remove faces
311  CHKERR trimSurface(fixed_edges, corner_nodes, tol_trim_close, th, debug);
312 
313  if (debug) {
315  Range bit2_edges;
316  CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
317  trim_bit, BitRefLevel().set(), MBEDGE, bit2_edges);
318  CHKERR SaveData(moab)("trim_fixed_edges.vtk",
319  intersect(*fixed_edges, bit2_edges));
320  }
321 
323 }
MoFEMErrorCode trimEdgesInTheMiddle(const BitRefLevel bit, const bool debug=false)
trim edges
virtual moab::Interface & get_moab()=0
#define CutMeshFunctionBegin
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:602
MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, const bool debug=false)
Find edges to trimEdges.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
MoFEMErrorCode trimSurface(Range *fixed_edge, Range *corner_nodes, const double tol=1e-4, Tag th=NULL, const bool debug=false)
Trim surface from faces beyond front.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
MoFEMErrorCode saveTrimEdges()
MoFEMErrorCode moveMidNodesOnTrimmedEdges(Tag th=NULL)
move trimmed edges mid nodes

◆ trimSurface()

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

Trim surface from faces beyond front.

Parameters
fixed_edge
corner_nodes
debug
Returns
MoFEMErrorCode

Definition at line 2047 of file CutMeshInterface.cpp.

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

Member Data Documentation

◆ aveLength

double MoFEM::CutMeshInterface::aveLength
private

Average edge length.

Definition at line 586 of file CutMeshInterface.hpp.

◆ constrainSurface

Range MoFEM::CutMeshInterface::constrainSurface
private

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

◆ cutFrontVolumes

Range MoFEM::CutMeshInterface::cutFrontVolumes
private

Definition at line 590 of file CutMeshInterface.hpp.

◆ cutNewSurfaces

Range MoFEM::CutMeshInterface::cutNewSurfaces
private

Definition at line 561 of file CutMeshInterface.hpp.

◆ cutNewVertices

Range MoFEM::CutMeshInterface::cutNewVertices
private

Definition at line 564 of file CutMeshInterface.hpp.

◆ cutNewVolumes

Range MoFEM::CutMeshInterface::cutNewVolumes
private

Definition at line 560 of file CutMeshInterface.hpp.

◆ cutSurfaceVolumes

Range MoFEM::CutMeshInterface::cutSurfaceVolumes
private

Definition at line 589 of file CutMeshInterface.hpp.

◆ edgesToCut

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

Definition at line 581 of file CutMeshInterface.hpp.

◆ edgesToTrim

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

Definition at line 583 of file CutMeshInterface.hpp.

◆ fRont

Range MoFEM::CutMeshInterface::fRont
private

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

◆ mergedSurfaces

Range MoFEM::CutMeshInterface::mergedSurfaces
private

Definition at line 572 of file CutMeshInterface.hpp.

◆ mergedVolumes

Range MoFEM::CutMeshInterface::mergedVolumes
private

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

◆ sUrface

Range MoFEM::CutMeshInterface::sUrface
private

Definition at line 553 of file CutMeshInterface.hpp.

◆ treeSurfPtr

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

Definition at line 556 of file CutMeshInterface.hpp.

◆ trimEdges

Range MoFEM::CutMeshInterface::trimEdges
private

Definition at line 569 of file CutMeshInterface.hpp.

◆ trimNewSurfaces

Range MoFEM::CutMeshInterface::trimNewSurfaces
private

Definition at line 568 of file CutMeshInterface.hpp.

◆ trimNewVertices

Range MoFEM::CutMeshInterface::trimNewVertices
private

Definition at line 567 of file CutMeshInterface.hpp.

◆ trimNewVolumes

Range MoFEM::CutMeshInterface::trimNewVolumes
private

Definition at line 566 of file CutMeshInterface.hpp.

◆ verticesOnCutEdges

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

Definition at line 582 of file CutMeshInterface.hpp.

◆ verticesOnTrimEdges

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

Definition at line 584 of file CutMeshInterface.hpp.

◆ vOlume

Range MoFEM::CutMeshInterface::vOlume
private

Definition at line 554 of file CutMeshInterface.hpp.

◆ zeroDistanceEnts

Range MoFEM::CutMeshInterface::zeroDistanceEnts
private

Definition at line 562 of file CutMeshInterface.hpp.

◆ zeroDistanceVerts

Range MoFEM::CutMeshInterface::zeroDistanceVerts
private

Definition at line 563 of file CutMeshInterface.hpp.


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