v0.8.20
Classes | 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  TreeData
 

Public Member Functions

MoFEMErrorCode query_interface (const MOFEMuuid &uuid, UnknownInterface **iface) const
 
 CutMeshInterface (const MoFEM::Core &core)
 
 ~CutMeshInterface ()
 
MoFEMErrorCode getSubInterfaceOptions ()
 
MoFEMErrorCode getOptions ()
 Get options from command line. More...
 
MoFEMErrorCode setFront (const Range &surface)
 
MoFEMErrorCode setSurface (const Range &surface)
 set surface entities More...
 
MoFEMErrorCode copySurface (const Range &surface, Tag th=NULL, double *shift=NULL, double *origin=NULL, double *transform=NULL, const std::string save_mesh="")
 copy surface entities More...
 
MoFEMErrorCode setVolume (const Range &volume)
 set volume entities More...
 
MoFEMErrorCode mergeSurface (const Range &surface)
 merge surface entities More...
 
MoFEMErrorCode mergeVolumes (const Range &volume)
 merge volume entities More...
 
MoFEMErrorCode snapSurfaceSkinToEdges (const Range &fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
 
MoFEMErrorCode snapSurfaceToEdges (const Range &surface_edges, const Range &fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
 
MoFEMErrorCode buildTree ()
 build tree More...
 
MoFEMErrorCode cutAndTrim (int &first_bit, const int before_trim_levels, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
 
MoFEMErrorCode cutTrimAndMerge (int &first_bit, const int fraction_level, const int before_trim_levels, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)
 
MoFEMErrorCode refCutTrimAndMerge (int &first_bit, const int fraction_level, const int ref_before_cut_levels, const int befor_trim_levels, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)
 
MoFEMErrorCode findEdgesToCut (Range *fixed_edges, Range *corner_nodes, const double low_tol, int verb=QUIET, const bool debug=false)
 find edges to cut More...
 
MoFEMErrorCode refineBeforeCut (const BitRefLevel &bit, Range *fixed_edges, const bool update_meshsets, const bool debug=false)
 
MoFEMErrorCode refineBeforeTrim (const BitRefLevel &bit, Range *fixed_edges, const bool update_meshsets, const bool debug=false)
 
MoFEMErrorCode projectZeroDistanceEnts (Range *fixed_edges, Range *corner_nodes, const double low_tol=0, const int verb=QUIET, const bool debug=false)
 
MoFEMErrorCode cutEdgesInMiddle (const BitRefLevel bit, 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 double tol_close=0, const bool debug=false)
 Find edges to trimEdges. More...
 
MoFEMErrorCode trimEdgesInTheMiddle (const BitRefLevel bit, Tag th=NULL, const double tol=1e-4, const bool debug=false)
 trim edges More...
 
MoFEMErrorCode moveMidNodesOnTrimmedEdges (Tag th=NULL)
 move trimmed edges mid nodes 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 rebuildMeshWithTetGen (vector< string > &switches, const BitRefLevel &mesh_bit, const BitRefLevel &bit, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Tag th=NULL, const bool debug=false)
 
MoFEMErrorCode setTagData (Tag th, const BitRefLevel bit=BitRefLevel())
 set coords to tag More...
 
MoFEMErrorCode setCoords (Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
 set coords from tag More...
 
const Range & getVolume () const
 
const Range & getSurface () const
 
const Range & getFront () const
 
const Range & getCutEdges () const
 
const Range & getCutVolumes () const
 
const Range & getNewCutVolumes () const
 
const Range & getNewCutSurfaces () const
 
const Range & getNewCutVertices () const
 
const Range & projectZeroDistanceEnts () const
 
const Range & getTrimEdges () const
 
const Range & getNewTrimVolumes () const
 
const Range & getNewTrimSurfaces () const
 
const Range & getNewTrimVertices () const
 
const Range & getMergedVolumes () const
 
const Range & getMergedSurfaces () const
 
const Range & getTetgenSurfaces () 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 ()
 
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 cutVolumes
 
Range cutNewVolumes
 
Range cutNewSurfaces
 
Range zeroDistanceEnts
 
Range zeroDistanceVerts
 
Range cutNewVertices
 
Range trimNewVolumes
 
Range trimNewVertices
 
Range trimNewSurfaces
 
Range trimEdges
 
Range mergedVolumes
 
Range mergedSurfaces
 
Range tetgenSurfaces
 
map< EntityHandle, TreeDataedgesToCut
 
map< EntityHandle, TreeDataverticesOnCutEdges
 
map< EntityHandle, TreeDataedgesToTrim
 
map< EntityHandle, TreeDataverticesOnTrimEdges
 
map< EntityHandle, unsigned long > moabTetGenMap
 
map< unsigned long, EntityHandletetGenMoabMap
 
boost::ptr_vector< tetgenio > tetGenData
 
double aveLength
 Average edge length. More...
 
double maxLength
 Maximal edge length. More...
 

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.

Constructor & Destructor Documentation

◆ CutMeshInterface()

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

Definition at line 52 of file CutMeshInterface.cpp.

53  : cOre(const_cast<Core &>(core)) {
54  lineSearchSteps = 10;
55  nbMaxMergingCycles = 200;
57 }

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

246  {
247  CoreInterface &m_field = cOre;
248  moab::Interface &moab = m_field.get_moab();
250  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
251  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
254 }
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ clearMap()

MoFEMErrorCode MoFEM::CutMeshInterface::clearMap ( )

Definition at line 59 of file CutMeshInterface.cpp.

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

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

80  {
81  CoreInterface &m_field = cOre;
82  moab::Interface &moab = m_field.get_moab();
84  std::map<EntityHandle, EntityHandle> verts_map;
85  for (Range::const_iterator tit = surface.begin(); tit != surface.end();
86  tit++) {
87  int num_nodes;
88  const EntityHandle *conn;
89  CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
90  MatrixDouble coords(num_nodes, 3);
91  if (th) {
92  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords(0, 0));
93  } else {
94  CHKERR moab.get_coords(conn, num_nodes, &coords(0, 0));
95  }
96  EntityHandle new_verts[num_nodes];
97  for (int nn = 0; nn != num_nodes; nn++) {
98  if (verts_map.find(conn[nn]) != verts_map.end()) {
99  new_verts[nn] = verts_map[conn[nn]];
100  } else {
101  if (transform) {
102  ublas::matrix_row<MatrixDouble> mr(coords, nn);
103  if (origin) {
104  VectorAdaptor vec_origin(
105  3, ublas::shallow_array_adaptor<double>(3, origin));
106  mr = mr - vec_origin;
107  }
108  MatrixAdaptor mat_transform = MatrixAdaptor(
109  3, 3, ublas::shallow_array_adaptor<double>(9, transform));
110  mr = prod(mat_transform, mr);
111  if (origin) {
112  VectorAdaptor vec_origin(
113  3, ublas::shallow_array_adaptor<double>(3, origin));
114  mr = mr + vec_origin;
115  }
116  }
117  if (shift) {
118  ublas::matrix_row<MatrixDouble> mr(coords, nn);
119  VectorAdaptor vec_shift(
120  3, ublas::shallow_array_adaptor<double>(3, shift));
121  mr = mr + vec_shift;
122  }
123  CHKERR moab.create_vertex(&coords(nn, 0), new_verts[nn]);
124  verts_map[conn[nn]] = new_verts[nn];
125  }
126  }
127  EntityHandle ele;
128  CHKERR moab.create_element(MBTRI, new_verts, num_nodes, ele);
129  sUrface.insert(ele);
130  }
131  if (!save_mesh.empty())
132  CHKERR SaveData(m_field.get_moab())(save_mesh, sUrface);
134 }
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
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:141
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:99
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ cutAndTrim()

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

Definition at line 265 of file CutMeshInterface.cpp.

269  {
270  CoreInterface &m_field = cOre;
271  moab::Interface &moab = m_field.get_moab();
273 
274  std::vector<BitRefLevel> bit_levels;
275 
276  auto get_back_bit_levels = [&]() {
277  bit_levels.push_back(BitRefLevel());
278  bit_levels.back().set(first_bit + bit_levels.size() - 1);
279  return bit_levels.back();
280  };
281 
282  BitRefLevel bit_cut = get_back_bit_levels();
283 
284  // cut mesh
285  CHKERR findEdgesToCut(fixed_edges, corner_nodes, tol_cut, QUIET, debug);
286  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_cut_close,
287  QUIET, debug);
288  CHKERR cutEdgesInMiddle(bit_cut, debug);
289  if (fixed_edges) {
290  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
291  *fixed_edges);
292  }
293  if (corner_nodes) {
294  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
295  *corner_nodes);
296  }
297  if (update_meshsets) {
298  CHKERR UpdateMeshsets()(cOre, bit_cut);
299  }
301 
302  auto get_min_quality = [&m_field](const BitRefLevel bit, Tag th) {
303  Range tets_level; // test at level
304  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
305  bit, BitRefLevel().set(), MBTET, tets_level);
306  double min_q = 1;
307  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
308  return min_q;
309  };
310 
311  PetscPrintf(PETSC_COMM_WORLD, "Min quality cut %6.4g\n",
312  get_min_quality(bit_cut, th));
313 
314  if (debug)
316 
317  Range starting_volume = cutNewVolumes;
318  Range starting_volume_nodes;
319  CHKERR m_field.get_moab().get_connectivity(starting_volume,
320  starting_volume_nodes, true);
321  std::vector<double> staring_volume_coords(3 * starting_volume_nodes.size());
322  CHKERR m_field.get_moab().get_coords(starting_volume_nodes,
323  &*staring_volume_coords.begin());
324 
325  if (th) {
326  std::vector<double> staring_volume_th_coords(3 *
327  starting_volume_nodes.size());
328  CHKERR m_field.get_moab().tag_get_data(th, starting_volume_nodes,
329  &*staring_volume_th_coords.begin());
330  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
331  &*staring_volume_th_coords.begin());
332  }
333 
334  if (th)
335  CHKERR setTagData(th);
336 
337  for (int ll = 0; ll != before_trim_levels; ++ll) {
338  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim,
339  tol_trim_close, debug);
340  BitRefLevel bit = get_back_bit_levels();
341  CHKERR refineBeforeTrim(bit, fixed_edges, true);
342  if (th)
343  CHKERR setTagData(th);
344  }
345 
346  BitRefLevel bit_trim = get_back_bit_levels();
347  // trim mesh
348  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim,
349  tol_trim_close, debug);
350  CHKERR trimEdgesInTheMiddle(bit_trim, th, tol_trim_close, debug);
351  if (fixed_edges) {
352  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
353  *fixed_edges);
354  }
355  if (corner_nodes) {
356  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
357  *corner_nodes);
358  }
359  if (update_meshsets) {
360  CHKERR UpdateMeshsets()(cOre, bit_trim);
361  }
363 
364  PetscPrintf(PETSC_COMM_WORLD, "Min quality trim %3.2g\n",
365  get_min_quality(bit_trim, th));
366 
367  if (debug) {
369  Range bit2_edges;
370  CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
371  bit_trim, BitRefLevel().set(), MBEDGE, bit2_edges);
372  CHKERR SaveData(moab)("trim_fixed_edges.vtk",
373  intersect(*fixed_edges, bit2_edges));
374  }
375 
376  first_bit += bit_levels.size() - 1;
377 
378  if (th)
379  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
380  &*staring_volume_coords.begin());
381 
383 }
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, const double tol_close=0, const bool debug=false)
Find edges to trimEdges.
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="")
MoFEMErrorCode findEdgesToCut(Range *fixed_edges, Range *corner_nodes, const double low_tol, int verb=QUIET, const bool debug=false)
find edges to cut
static const bool debug
MoFEMErrorCode moveMidNodesOnCutEdges(Tag th=NULL)
projecting of mid edge nodes on new mesh on surface
MoFEMErrorCode cutEdgesInMiddle(const BitRefLevel bit, const bool debug=false)
cut edges
#define CHKERR
Inline error check.
Definition: definitions.h:594
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
MoFEMErrorCode saveTrimEdges()
MoFEMErrorCode trimEdgesInTheMiddle(const BitRefLevel bit, Tag th=NULL, const double tol=1e-4, const bool debug=false)
trim edges
MoFEMErrorCode refineBeforeTrim(const BitRefLevel &bit, Range *fixed_edges, const bool update_meshsets, const bool debug=false)
MoFEMErrorCode moveMidNodesOnTrimmedEdges(Tag th=NULL)
move trimmed edges mid nodes

◆ cutEdgesInMiddle()

MoFEMErrorCode MoFEM::CutMeshInterface::cutEdgesInMiddle ( const BitRefLevel  bit,
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 1272 of file CutMeshInterface.cpp.

1273  {
1274  CoreInterface &m_field = cOre;
1275  moab::Interface &moab = m_field.get_moab();
1276  MeshRefinement *refiner;
1277  const RefEntity_multiIndex *ref_ents_ptr;
1279  if (cutEdges.size() != edgesToCut.size()) {
1280  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1281  }
1282  CHKERR m_field.getInterface(refiner);
1283  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1284  CHKERR refiner->add_verices_in_the_middel_of_edges(cutEdges, bit);
1285  CHKERR refiner->refine_TET(vOlume, bit, false, QUIET,
1286  debug ? &cutEdges : NULL);
1287  cutNewVolumes.clear();
1288  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1289  bit, BitRefLevel().set(), MBTET, cutNewVolumes);
1290  cutNewSurfaces.clear();
1291  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1292  bit, BitRefLevel().set(), MBTRI, cutNewSurfaces);
1293  // Find new vertices on cut edges
1294  cutNewVertices.clear();
1295  CHKERR moab.get_connectivity(zeroDistanceEnts, cutNewVertices, true);
1297  for (auto &m : edgesToCut) {
1298  RefEntity_multiIndex::index<
1299  Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vit =
1300  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1301  boost::make_tuple(MBVERTEX, m.first));
1302  if (vit ==
1303  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1304  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1305  "No vertex on cut edges, that make no sense");
1306  }
1307  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1308  if ((ref_ent->getBitRefLevel() & bit).any()) {
1309  EntityHandle vert = ref_ent->getRefEnt();
1310  cutNewVertices.insert(vert);
1311  verticesOnCutEdges[vert] = m.second;
1312  } else {
1313  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1314  "Vertex has wrong bit ref level");
1315  }
1316  }
1317 
1318  // Add zero distance entities faces
1319  Range tets_skin;
1320  Skinner skin(&moab);
1321  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
1322  cutNewSurfaces.merge(zeroDistanceEnts.subset_by_type(MBTRI));
1323 
1324  // At that point cutNewSurfaces has all newly created faces, now take all
1325  // nodes on those faces and subtract nodes on catted edges. Faces adjacent to
1326  // nodes which left are not part of surface.
1327  Range diff_verts;
1328  CHKERR moab.get_connectivity(unite(cutNewSurfaces, zeroDistanceEnts),
1329  diff_verts, true);
1330  diff_verts = subtract(diff_verts, cutNewVertices);
1331  Range subtract_faces;
1332  CHKERR moab.get_adjacencies(diff_verts, 2, false, subtract_faces,
1333  moab::Interface::UNION);
1334  cutNewSurfaces = subtract(cutNewSurfaces, unite(subtract_faces, tets_skin));
1335  cutNewVertices.clear();
1336  CHKERR moab.get_connectivity(cutNewSurfaces, cutNewVertices, true);
1337 
1339 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:594
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
map< EntityHandle, TreeData > edgesToCut

◆ cutTrimAndMerge()

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

Definition at line 385 of file CutMeshInterface.cpp.

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

◆ findEdgesToCut()

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

find edges to cut

Parameters
verbverbosity level
Returns
error code

Definition at line 519 of file CutMeshInterface.cpp.

522  {
523  CoreInterface &m_field = cOre;
524  moab::Interface &moab = m_field.get_moab();
526 
527  edgesToCut.clear();
528  cutEdges.clear();
529 
530  zeroDistanceVerts.clear();
531  zeroDistanceEnts.clear();
532  verticesOnCutEdges.clear();
533 
534  double ray_length;
535  double ray_point[3], unit_ray_dir[3];
536  VectorAdaptor vec_unit_ray_dir(
537  3, ublas::shallow_array_adaptor<double>(3, unit_ray_dir));
538  VectorAdaptor vec_ray_point(
539  3, ublas::shallow_array_adaptor<double>(3, ray_point));
540 
541  Tag th_dist;
542  rval = moab.tag_get_handle("DIST", th_dist);
543  if (rval == MB_SUCCESS)
544  CHKERR moab.tag_delete(th_dist);
545  else
546  rval = MB_SUCCESS;
547 
548  Tag th_dist_normal;
549  rval = moab.tag_get_handle("DIST_NORMAL", th_dist_normal);
550  if (rval == MB_SUCCESS)
551  CHKERR moab.tag_delete(th_dist_normal);
552  else
553  rval = MB_SUCCESS;
554 
555  double def_val[] = {0, 0, 0};
556  CHKERR moab.tag_get_handle("DIST", 1, MB_TYPE_DOUBLE, th_dist,
557  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
558  CHKERR moab.tag_get_handle("DIST_NORMAL", 3, MB_TYPE_DOUBLE, th_dist_normal,
559  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
560 
561  // Calculate not signed and not signed distances from nodes to surface
562  Range vol_vertices;
563  CHKERR moab.get_connectivity(vOlume, vol_vertices, true);
564  for (auto v : vol_vertices) {
565  VectorDouble3 point_in(3);
566  CHKERR moab.get_coords(&v, 1, &point_in[0]);
567  VectorDouble3 point_out(3);
568  EntityHandle facets_out;
569  CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
570  &point_out[0], facets_out);
571  VectorDouble3 n(3);
572  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
573  VectorDouble3 delta = point_out - point_in;
574  double dist = norm_2(delta);
575  n /= norm_2(n);
576  VectorDouble3 dist_normal = inner_prod(delta, n) * n;
577  // not signed distance
578  CHKERR moab.tag_set_data(th_dist, &v, 1, &dist);
579  // signed distance
580  CHKERR moab.tag_set_data(th_dist_normal, &v, 1, &dist_normal[0]);
581  }
582 
583  auto get_normal_dist = [](const double *normal) {
584  FTensor::Tensor1<const double *, 3> t_n(normal, &normal[1], &normal[2]);
586  return sqrt(t_n(i) * t_n(i));
587  };
588 
589  auto get_edge_crossed = [&moab, get_normal_dist,
590  th_dist_normal](EntityHandle v0, EntityHandle v1) {
591  VectorDouble3 dist_normal0(3);
592  CHKERR moab.tag_get_data(th_dist_normal, &v0, 1, &dist_normal0[0]);
593  VectorDouble3 dist_normal1(3);
594  CHKERR moab.tag_get_data(th_dist_normal, &v1, 1, &dist_normal1[0]);
595  return (inner_prod(dist_normal0, dist_normal1) < 0);
596  };
597 
598  auto get_normal_dist_from_conn = [&moab, get_normal_dist,
599  th_dist_normal](EntityHandle v) {
600  double dist_normal[3];
601  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, dist_normal);
602  return get_normal_dist(dist_normal);
603  };
604 
605  auto get_ray_for_edge = [&](const EntityHandle ent, VectorAdaptor &ray_point,
606  VectorAdaptor &unit_ray_dir, double &ray_length) {
608  int num_nodes;
609  const EntityHandle *conn;
610  CHKERR moab.get_connectivity(ent, conn, num_nodes, true);
611  double coords[6];
612  CHKERR moab.get_coords(conn, num_nodes, coords);
613  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
614  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
615  noalias(ray_point) = s0;
616  noalias(unit_ray_dir) = s1 - s0;
617  ray_length = norm_2(unit_ray_dir);
618  if (ray_length > 0)
619  unit_ray_dir /= ray_length;
621  };
622 
623  auto set_edge_cut_point_from_level_set = [&](const EntityHandle e,
624  const auto dist_normal) {
626  CHKERR get_ray_for_edge(e, vec_ray_point, vec_unit_ray_dir, ray_length);
627  const double s = dist_normal[0] / (dist_normal[1] + dist_normal[0]);
628  edgesToCut[e].dIst = s * ray_length;
629  edgesToCut[e].lEngth = ray_length;
630  edgesToCut[e].unitRayDir = vec_unit_ray_dir;
631  edgesToCut[e].rayPoint = vec_ray_point;
632  cutEdges.insert(e);
634  };
635 
636  auto not_project_node = [this, &moab](const EntityHandle v) {
638  VectorDouble3 s0(3);
639  CHKERR moab.get_coords(&v, 1, &s0[0]);
640  verticesOnCutEdges[v].dIst = 0;
641  verticesOnCutEdges[v].lEngth = 0;
642  verticesOnCutEdges[v].unitRayDir.resize(3, false);
643  verticesOnCutEdges[v].unitRayDir.clear();
644  verticesOnCutEdges[v].rayPoint = s0;
646  };
647 
648  auto get_ave_edge_length = [&](const EntityHandle ent,
649  const Range &vol_edges) {
650  Range adj_edges;
651  if (moab.type_from_handle(ent) == MBVERTEX) {
652  CHKERR moab.get_adjacencies(&ent, 1, 1, false, adj_edges);
653  } else {
654  adj_edges.insert(ent);
655  }
656  adj_edges = intersect(adj_edges, vol_edges);
657  double ave_l = 0;
658  for (auto e : adj_edges) {
659  int num_nodes;
660  const EntityHandle *conn;
661  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
662  VectorDouble6 coords(6);
663  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
665  &coords[0], &coords[1], &coords[2]);
667  &coords[3], &coords[4], &coords[5]);
669  t_n0(i) -= t_n1(i);
670  ave_l += sqrt(t_n0(i) * t_n0(i));
671  }
672  return ave_l / adj_edges.size();
673  };
674 
675  Range vol_edges;
676  CHKERR moab.get_adjacencies(vOlume, 1, true, vol_edges,
677  moab::Interface::UNION);
678 
679  aveLength = 0;
680  maxLength = 0;
681  int nb_ave_length = 0;
682  for (auto e : vol_edges) {
683  int num_nodes;
684  const EntityHandle *conn;
685  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
686  double dist[num_nodes];
687  CHKERR moab.tag_get_data(th_dist, conn, num_nodes, dist);
688  CHKERR get_ray_for_edge(e, vec_ray_point, vec_unit_ray_dir, ray_length);
689  const double tol = ray_length * low_tol;
690 
691  // Edges is on two sides of the surface
692  if (get_edge_crossed(conn[0], conn[1])) {
693  std::vector<double> distances_out;
694  std::vector<EntityHandle> facets_out;
695  // sen the ray in the direction of the edge
696  CHKERR treeSurfPtr->ray_intersect_triangles(distances_out, facets_out,
697  rootSetSurf, tol, ray_point,
698  unit_ray_dir, &ray_length);
699  // if ray punching through the surface
700  if (!distances_out.empty()) {
701  // minimal distance
702  const auto dist_ptr =
703  std::min_element(distances_out.begin(), distances_out.end());
704  const double dist = *dist_ptr;
705  // if distance is shorter than edge length
706  if (dist <= ray_length) {
707  aveLength += ray_length;
708  maxLength = std::max(maxLength, ray_length);
709  nb_ave_length++;
710  edgesToCut[e].dIst = dist;
711  edgesToCut[e].lEngth = ray_length;
712  edgesToCut[e].unitRayDir = vec_unit_ray_dir;
713  edgesToCut[e].rayPoint = vec_ray_point;
714  cutEdges.insert(e);
715  }
716  }
717  }
718 
719  // check if cutting point is close to the end of the edges
720  if (fabs(dist[0]) < tol && fabs(dist[1]) < tol) {
721  aveLength += ray_length;
722  maxLength = fmax(maxLength, ray_length);
723  nb_ave_length++;
724  for (int n : {0, 1})
725  CHKERR not_project_node(conn[n]);
726  zeroDistanceEnts.insert(e);
727  }
728  }
729  aveLength /= nb_ave_length;
730 
731  for (auto v : vol_vertices) {
732  double dist;
733  CHKERR moab.tag_get_data(th_dist, &v, 1, &dist);
734  const double tol = get_ave_edge_length(v, vol_edges) * low_tol;
735  if (fabs(dist) < tol) {
736  CHKERR not_project_node(v);
737  zeroDistanceVerts.insert(v);
738  }
739  }
740 
741  cutVolumes.clear();
742  Range cut_edges_verts;
743  CHKERR moab.get_connectivity(unite(cutEdges, zeroDistanceEnts),
744  cut_edges_verts, true);
745  CHKERR moab.get_adjacencies(unite(cut_edges_verts, zeroDistanceVerts), 3,
746  false, cutVolumes, moab::Interface::UNION);
747  cutVolumes = intersect(cutVolumes, vOlume);
748 
749  // get edges on the cut volumes
750  Range edges;
751  CHKERR moab.get_adjacencies(cutVolumes, 1, false, edges,
752  moab::Interface::UNION);
753  edges = subtract(edges, cutEdges);
754 
755  // add to cut set edges which are cutted by extension of cutting surface
756  for (auto e : edges) {
757  int num_nodes;
758  const EntityHandle *conn;
759  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
760  const double tol = get_ave_edge_length(e, vol_edges) * low_tol;
761  double dist_normal[2];
762  dist_normal[0] = get_normal_dist_from_conn(conn[0]);
763  dist_normal[1] = get_normal_dist_from_conn(conn[1]);
764  if (get_edge_crossed(conn[0], conn[1]))
765  CHKERR set_edge_cut_point_from_level_set(e, dist_normal);
766  else if (fabs(dist_normal[0]) < tol && fabs(dist_normal[1]) < tol) {
767  for (int n : {0, 1})
768  CHKERR not_project_node(conn[n]);
769  zeroDistanceEnts.insert(e);
770  }
771  }
772 
773  for (auto f : zeroDistanceEnts) {
774  int num_nodes;
775  const EntityHandle *conn;
776  CHKERR moab.get_connectivity(f, conn, num_nodes, true);
777  Range adj_edges;
778  CHKERR moab.get_adjacencies(conn, num_nodes, 1, false, adj_edges,
779  moab::Interface::UNION);
780  for (auto e : adj_edges) {
781  cutEdges.erase(e);
782  edgesToCut.erase(e);
783  }
784  }
785 
786  for (auto v : zeroDistanceVerts) {
787  Range adj_edges;
788  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges,
789  moab::Interface::UNION);
790  for (auto e : adj_edges) {
791  cutEdges.erase(e);
792  edgesToCut.erase(e);
793  }
794  }
795 
796  if (debug)
797  CHKERR SaveData(m_field.get_moab())("cut_edges.vtk", cutEdges);
798 
799  if (debug)
800  CHKERR SaveData(m_field.get_moab())(
801  "cut_zero_dist.vtk", unite(zeroDistanceVerts, zeroDistanceEnts));
802 
804 }
double maxLength
Maximal edge length.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
double aveLength
Average edge length.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:86
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
double tol
static const bool debug
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:99
#define CHKERR
Inline error check.
Definition: definitions.h:594
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
map< EntityHandle, TreeData > edgesToCut

◆ findEdgesToTrim()

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

1382  {
1383  CoreInterface &m_field = cOre;
1384  moab::Interface &moab = m_field.get_moab();
1386 
1387  // takes body skin
1388  Skinner skin(&moab);
1389  Range tets_skin;
1390  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
1391  // vertices on the skin
1392  Range tets_skin_verts;
1393  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1394  // edges on the skin
1395  Range tets_skin_edges;
1396  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1397  moab::Interface::UNION);
1398  // get edges on new surface
1399  Range edges;
1400  CHKERR moab.get_adjacencies(cutNewSurfaces, 1, false, edges,
1401  moab::Interface::UNION);
1402 
1403  Range cut_surface_edges_on_fixed_edges;
1404  if (fixed_edges)
1405  cut_surface_edges_on_fixed_edges = intersect(edges, *fixed_edges);
1406 
1407  Range cut_surface_edges_on_fixed_edges_verts;
1408  CHKERR moab.get_connectivity(cut_surface_edges_on_fixed_edges,
1409  cut_surface_edges_on_fixed_edges_verts, true);
1410 
1411  Range fixed_edges_verts;
1412  if (fixed_edges)
1413  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts, true);
1414 
1415  Range surface_skin;
1416  if (fRont.empty())
1417  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
1418  else
1419  surface_skin = fRont;
1420 
1421  auto get_point_coords = [&th, &moab](EntityHandle v) {
1422  VectorDouble3 point(3);
1423  if (th)
1424  CHKERR moab.tag_get_data(th, &v, 1, &point[0]);
1425  else
1426  CHKERR moab.get_coords(&v, 1, &point[0]);
1427  return point;
1428  };
1429 
1430  struct VertMap {
1431  const double d;
1432  const EntityHandle v;
1433  const EntityHandle e;
1434  VertMap(const double d, const EntityHandle v, const EntityHandle e)
1435  : d(d), v(v), e(e) {}
1436  };
1437 
1438  typedef multi_index_container<
1439  VertMap,
1440  indexed_by<
1441  ordered_non_unique<member<VertMap, const double, &VertMap::d>>,
1442  ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::v>>
1443 
1444  >>
1445  VertMapMultiIndex;
1446 
1447  VertMapMultiIndex verts_map;
1448 
1449  auto add_vert = [&](const EntityHandle v, const EntityHandle e,
1450  const double dist) {
1451  verts_map.insert(VertMap(dist, v, e));
1452  };
1453 
1454  // clear data ranges
1455  trimEdges.clear();
1456  edgesToTrim.clear();
1457  verticesOnTrimEdges.clear();
1458  trimNewVertices.clear();
1459 
1460  if (debug)
1461  CHKERR SaveData(moab)("edges_potentially_to_trim.vtk", edges);
1462 
1463  if (debug)
1464  CHKERR SaveData(moab)("surface_skin_to_trim.vtk", surface_skin);
1465 
1466  auto cut_this_edge = [&](const EntityHandle e, const double length, auto &ray,
1467  auto &ray_point) {
1468  trimEdges.insert(e);
1469  edgesToTrim[e].dIst = 1;
1470  edgesToTrim[e].lEngth = length;
1471  edgesToTrim[e].unitRayDir.resize(3, false);
1472  edgesToTrim[e].rayPoint.resize(3, false);
1473  for (int ii = 0; ii != 3; ++ii)
1474  edgesToTrim[e].unitRayDir[ii] = ray(ii);
1475  for (int ii = 0; ii != 3; ++ii)
1476  edgesToTrim[e].rayPoint[ii] = ray_point(ii);
1477  };
1478 
1480  int num_nodes;
1481 
1482  MatrixDouble edge_face_normal(3, surface_skin.size());
1483  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_edge_face_normal =
1484  getFTensor1FromMat<3>(edge_face_normal);
1485  for (auto s : surface_skin) {
1486  Range adj_face;
1487  moab.get_adjacencies(&s, 1, 2, false, adj_face, moab::Interface::UNION);
1488  adj_face = intersect(adj_face, sUrface);
1489  if (adj_face.size() == 1)
1490  Util::normal(&moab, adj_face[0], t_edge_face_normal(0),
1491  t_edge_face_normal(1), t_edge_face_normal(2));
1492  else
1493  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1494  "Should be only one face in the range");
1495 
1496  t_edge_face_normal(i) /=
1497  sqrt(t_edge_face_normal(i) * t_edge_face_normal(i));
1498  ++t_edge_face_normal;
1499  }
1500 
1501  for (auto e : edges) {
1502 
1503  // Get edge connectivity and coords
1504  const EntityHandle *conn_edge;
1505  CHKERR moab.get_connectivity(e, conn_edge, num_nodes, true);
1506  double coords_edge[3 * num_nodes];
1507  CHKERR moab.get_coords(conn_edge, num_nodes, coords_edge);
1508 
1509  FTensor::Tensor1<double *, 3> t_e0(&coords_edge[0], &coords_edge[1],
1510  &coords_edge[2]);
1511  FTensor::Tensor1<double *, 3> t_e1(&coords_edge[3], &coords_edge[4],
1512  &coords_edge[5]);
1513 
1514  FTensor::Tensor1<double, 3> t_edge_delta;
1515  t_edge_delta(i) = t_e1(i) - t_e0(i);
1516  const double edge_length2 = t_edge_delta(i) * t_edge_delta(i);
1517  const double edge_length = sqrt(edge_length2);
1518  if (edge_length == 0)
1519  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Zero edge length");
1520 
1521  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_edge_face_normal =
1522  getFTensor1FromMat<3>(edge_face_normal);
1523 
1524  auto t_project = [&](auto &t_vec) {
1526  t_prj(i) =
1527  t_vec(i) - (t_edge_face_normal(i) * t_vec(i)) * t_edge_face_normal(i);
1528  return t_prj;
1529  };
1530 
1531  auto get_edge_coors = [&](const auto e) {
1532  const EntityHandle *conn;
1533  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1534  VectorDouble6 coords(6);
1535  if (th)
1536  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1537  else
1538  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1539  return coords;
1540  };
1541 
1542  // iterate entities on inner surface font to find that edge cross
1543  for (auto s : surface_skin) {
1544 
1545  auto coords_front = get_edge_coors(s);
1546 
1547  FTensor::Tensor1<double *, 3> t_f0(&coords_front[0], &coords_front[1],
1548  &coords_front[2]);
1549  FTensor::Tensor1<double *, 3> t_f1(&coords_front[3], &coords_front[4],
1550  &coords_front[5]);
1551 
1552  auto t_e0_prj = t_project(t_e0);
1553  auto t_e1_prj = t_project(t_e1);
1554  auto t_f0_prj = t_project(t_f0);
1555  auto t_f1_prj = t_project(t_f1);
1556 
1557  double t_edge, t_front;
1558  auto res = Tools::minDistanceFromSegments(&t_e0_prj(0), &t_e1_prj(0),
1559  &t_f0_prj(0), &t_f1_prj(0),
1560  &t_edge, &t_front);
1561 
1562  if (res != Tools::NO_SOLUTION) {
1563 
1564  const double overlap_tol = 1e-2;
1565  if (t_edge >= 0 && t_edge <= 1 && t_front >= -overlap_tol &&
1566  t_front <= 1 + overlap_tol) {
1567 
1568  FTensor::Tensor1<double, 3> t_front_delta_prj;
1569  t_front_delta_prj(i) = t_f1_prj(i) - t_f0_prj(i);
1570  FTensor::Tensor1<double, 3> t_edge_delta_prj;
1571  t_edge_delta_prj(i) = t_e1_prj(i) - t_e0_prj(i);
1572  FTensor::Tensor1<double, 3> t_edge_point_prj, t_front_point_prj;
1573  t_edge_point_prj(i) = t_e0_prj(i) + t_edge * t_edge_delta_prj(i);
1574  t_front_point_prj(i) = t_f0_prj(i) + t_front * t_front_delta_prj(i);
1575  FTensor::Tensor1<double, 3> t_ray_prj;
1576  t_ray_prj(i) = t_front_point_prj(i) - t_edge_point_prj(i);
1577  const double dist = sqrt(t_ray_prj(i) * t_ray_prj(i));
1578 
1579  if ((dist / edge_length) < 0.1) {
1580 
1581  auto check_to_add_edge = [&](const EntityHandle e,
1582  const double dist) {
1583  auto eit = edgesToTrim.find(e);
1584  if (eit != edgesToTrim.end())
1585  if (eit->second.dIst < dist)
1586  return false;
1587  return true;
1588  };
1589 
1591  if (t_edge < 0.5) {
1592  t_ray(i) = t_edge * t_edge_delta(i);
1593  const double ray_length = sqrt(t_ray(i) * t_ray(i));
1594  if (check_to_add_edge(e, ray_length)) {
1595  add_vert(conn_edge[0], e, fabs(t_edge));
1596  add_vert(conn_edge[1], e, fabs(t_edge - 1));
1597  cut_this_edge(e, edge_length, t_ray, t_e0);
1598  }
1599  } else {
1600  FTensor::Tensor1<double, 3> t_edge_point;
1601  t_edge_point(i) = t_e0(i) + t_edge * t_edge_delta(i);
1602  t_ray(i) = t_edge_point(i) - t_e1(i);
1603  const double ray_length = sqrt(t_ray(i) * t_ray(i));
1604  if (check_to_add_edge(e, ray_length)) {
1605  add_vert(conn_edge[0], e, fabs(t_edge));
1606  add_vert(conn_edge[1], e, fabs(t_edge - 1));
1607  cut_this_edge(e, edge_length, t_ray, t_e1);
1608  }
1609  }
1610  }
1611  }
1612  }
1613 
1614  ++t_edge_face_normal;
1615  }
1616  }
1617 
1618  if (debug)
1619  CHKERR SaveData(moab)("edges_selected_to_trim.vtk", trimEdges);
1620 
1621  Range corners;
1622  if (corner_nodes)
1623  corners = *corner_nodes;
1624  Range fix_edges;
1625  if (fixed_edges)
1626  fix_edges = *fixed_edges;
1627 
1628  // Iterate over all vertives close to surface front and check if those can be
1629  // moved
1630  auto lo_m = verts_map.get<0>().lower_bound(0);
1631  auto hi_m = verts_map.get<0>().lower_bound(tol);
1632  for (; lo_m != hi_m; ++lo_m) {
1633 
1634  const double dist = lo_m->d;
1635  EntityHandle v = lo_m->v;
1636  if (verticesOnTrimEdges.find(v) == verticesOnTrimEdges.end()) {
1637 
1638  VectorDouble3 ray_point = get_point_coords(v);
1639  EntityHandle e = lo_m->e;
1640 
1641  VectorDouble3 new_pos(3);
1642  new_pos.clear();
1643 
1644  auto get_quality_change = [&](const Range &adj_tets) {
1645  double q0 = 1;
1646  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0, th);
1647 
1648  double min_q = 1;
1649  for (auto t : adj_tets) {
1650  int num_nodes;
1651  const EntityHandle *conn;
1652  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1653  VectorDouble12 coords(12);
1654  if (th)
1655  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1656  else
1657  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1658 
1659  for (int n = 0; n != 4; ++n) {
1660  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1661  if (conn[n] == v) {
1662  noalias(n_coords) = new_pos;
1663  } else {
1664  auto m = verticesOnTrimEdges.find(conn[n]);
1665  if (m != verticesOnTrimEdges.end())
1666  noalias(n_coords) =
1667  m->second.rayPoint + m->second.dIst * m->second.unitRayDir;
1668  }
1669  }
1670 
1671  const double q = Tools::volumeLengthQuality(&coords[0]);
1672  if (!std::isnormal(q)) {
1673  min_q = -2;
1674  break;
1675  }
1676  min_q = std::min(min_q, q);
1677  }
1678  return min_q / q0;
1679  };
1680 
1681  auto add_trim_vert = [&]() {
1682  auto &r = edgesToTrim.at(e);
1683  noalias(new_pos) = r.rayPoint + r.dIst * r.unitRayDir;
1684  VectorDouble3 unit_ray_dir = ray_point - new_pos;
1685 
1686  Range adj_tets;
1687  CHKERR moab.get_adjacencies(&v, 1, 3, false, adj_tets);
1688  adj_tets = intersect(adj_tets, cutNewVolumes);
1689 
1690  double q = get_quality_change(adj_tets);
1692  VectorDouble3 ray_dir = new_pos - ray_point;
1693  double dist = norm_2(ray_dir);
1694  verticesOnTrimEdges[v].dIst = 1;
1695  verticesOnTrimEdges[v].lEngth = dist;
1696  verticesOnTrimEdges[v].unitRayDir = ray_dir;
1697  verticesOnTrimEdges[v].rayPoint = ray_point;
1698  trimNewVertices.insert(v);
1699  }
1700  };
1701 
1702  if (trimNewVertices.find(v) == trimNewVertices.end()) {
1703  if (edgesToTrim.find(e) != edgesToTrim.end()) {
1704  if (dist < tol_close) {
1705 
1706  int rank_vertex = 0;
1707  if (corners.find(v) != corners.end())
1708  ++rank_vertex;
1709  if (fixed_edges_verts.find(v) != fixed_edges_verts.end())
1710  ++rank_vertex;
1711  if (tets_skin_verts.find(v) != tets_skin_verts.end())
1712  ++rank_vertex;
1713 
1714  int rank_edge = 0;
1715  if (fix_edges.find(e) != fix_edges.end())
1716  ++rank_edge;
1717  if (tets_skin_edges.find(e) != tets_skin_edges.end())
1718  ++rank_edge;
1719 
1720  if (rank_vertex > rank_edge) {
1721  auto &m = edgesToTrim.at(e);
1722  verticesOnTrimEdges[v] = m;
1723  verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
1724  verticesOnTrimEdges[v].lEngth = 0;
1725  verticesOnTrimEdges[v].dIst = 0;
1726  trimNewVertices.insert(v);
1727  } else if (rank_vertex == rank_edge)
1728  add_trim_vert();
1729  else
1730  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1731  "Imposible case");
1732 
1733  } else if (corners.find(v) != corners.end()) {
1734  auto &m = edgesToTrim.at(e);
1735  verticesOnTrimEdges[v] = m;
1736  verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
1737  verticesOnTrimEdges[v].lEngth = 0;
1738  verticesOnTrimEdges[v].dIst = 0;
1739  trimNewVertices.insert(v);
1740  } else if (fixed_edges_verts.find(v) != fixed_edges_verts.end()) {
1741  if (fix_edges.find(e) != fix_edges.end()) {
1742  add_trim_vert();
1743  }
1744  } else if (tets_skin_verts.find(v) != tets_skin_verts.end()) {
1745  if (tets_skin_edges.find(e) != tets_skin_edges.end()) {
1746  add_trim_vert();
1747  }
1748  } else {
1749  add_trim_vert();
1750  }
1751  }
1752  }
1753  }
1754  }
1755 
1756  for (auto m : verticesOnTrimEdges) {
1757  EntityHandle v = m.first;
1758  // Range adj_vertex_edges;
1759  // CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_vertex_edges,
1760  // moab::Interface::UNION);
1761  // for(auto e : adj_vertex_edges) {
1762  // edgesToTrim.erase(e);
1763  // trimEdges.erase(e);
1764  // }
1765  auto lo = verts_map.get<1>().lower_bound(v);
1766  auto hi = verts_map.get<1>().upper_bound(v);
1767  for (; lo != hi; ++lo) {
1768  edgesToTrim.erase(lo->e);
1769  trimEdges.erase(lo->e);
1770  }
1771  }
1772 
1773  if (!trimNewVertices.empty()) {
1774  Range adj_trim_vertex_edges;
1775  CHKERR moab.get_adjacencies(trimNewVertices, 1, false,
1776  adj_trim_vertex_edges, moab::Interface::UNION);
1777  adj_trim_vertex_edges = intersect(adj_trim_vertex_edges, trimEdges);
1778  Range adj_trim_vertex_edges_verts;
1779  CHKERR moab.get_connectivity(adj_trim_vertex_edges,
1780  adj_trim_vertex_edges_verts, true);
1781  adj_trim_vertex_edges_verts =
1782  subtract(adj_trim_vertex_edges_verts, trimNewVertices);
1783  Range edges_connected_to_close_verts;
1784  CHKERR moab.get_adjacencies(adj_trim_vertex_edges_verts, 1, false,
1785  edges_connected_to_close_verts,
1786  moab::Interface::UNION);
1787  edges_connected_to_close_verts =
1788  subtract(adj_trim_vertex_edges, edges_connected_to_close_verts);
1789 
1790  if (debug)
1791  if (!edges_connected_to_close_verts.empty())
1792  CHKERR SaveData(moab)("trim_edges_connected_to_close_verts.vtk",
1793  edges_connected_to_close_verts);
1794 
1795  for (auto e : edges_connected_to_close_verts) {
1796  edgesToTrim.erase(e);
1797  trimEdges.erase(e);
1798  }
1799  }
1800 
1801  if (debug)
1802  if (!trimNewVertices.empty())
1803  CHKERR SaveData(moab)("trim_close_vertices.vtk", trimNewVertices);
1804 
1805  if (debug)
1806  if (!trimEdges.empty())
1807  CHKERR SaveData(moab)("trim_edges.vtk", trimEdges);
1808 
1810 } // namespace MoFEM
map< EntityHandle, TreeData > edgesToTrim
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:32
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:284
map< EntityHandle, TreeData > verticesOnTrimEdges
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
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
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:86
auto getVectorAdaptor
Get Vector adaptor.
Definition: Types.hpp:121
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Types.hpp:88
double tol
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:594
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ getCutEdges()

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

Definition at line 324 of file CutMeshInterface.hpp.

324 { return cutEdges; }

◆ getCutVolumes()

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

Definition at line 325 of file CutMeshInterface.hpp.

325 { return cutVolumes; }

◆ getFront()

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

Definition at line 321 of file CutMeshInterface.hpp.

321 { return fRont; }

◆ getMergedSurfaces()

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

Definition at line 339 of file CutMeshInterface.hpp.

◆ getMergedVolumes()

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

Definition at line 338 of file CutMeshInterface.hpp.

338 { return mergedVolumes; }

◆ getNewCutSurfaces()

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

Definition at line 327 of file CutMeshInterface.hpp.

◆ getNewCutVertices()

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

Definition at line 328 of file CutMeshInterface.hpp.

◆ getNewCutVolumes()

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

Definition at line 326 of file CutMeshInterface.hpp.

326 { return cutNewVolumes; }

◆ getNewTrimSurfaces()

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

Definition at line 335 of file CutMeshInterface.hpp.

◆ getNewTrimVertices()

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

Definition at line 336 of file CutMeshInterface.hpp.

◆ getNewTrimVolumes()

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

Definition at line 334 of file CutMeshInterface.hpp.

◆ getOptions()

MoFEMErrorCode MoFEM::CutMeshInterface::getOptions ( )

Get options from command line.

Returns
error code

Definition at line 51 of file CutMeshInterface.hpp.

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

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

320 { return sUrface; }

◆ getTetgenSurfaces()

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

Definition at line 341 of file CutMeshInterface.hpp.

◆ getTreeSurfPtr()

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

Definition at line 347 of file CutMeshInterface.hpp.

347  {
348  return treeSurfPtr;
349  }
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr

◆ getTrimEdges()

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

Definition at line 333 of file CutMeshInterface.hpp.

333 { return trimEdges; }

◆ getVolume()

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

Definition at line 319 of file CutMeshInterface.hpp.

319 { return vOlume; }

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

2004  {
2005  CoreInterface &m_field = cOre;
2006  moab::Interface &moab = m_field.get_moab();
2008 
2009  /**
2010  * \brief Merge nodes
2011  */
2012  struct MergeNodes {
2013  CoreInterface &mField;
2014  const bool onlyIfImproveQuality;
2015  const int lineSearch;
2016  Tag tH;
2017  bool updateMehsets;
2018 
2019  MergeNodes(CoreInterface &m_field, const bool only_if_improve_quality,
2020  const int line_search, Tag th, bool update_mehsets)
2021  : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
2022  lineSearch(line_search), tH(th), updateMehsets(update_mehsets) {
2023  mField.getInterface(nodeMergerPtr);
2024  }
2025  NodeMergerInterface *nodeMergerPtr;
2026  MoFEMErrorCode operator()(EntityHandle father, EntityHandle mother,
2027  Range &proc_tets, Range &new_surf,
2028  Range &edges_to_merge, Range &not_merged_edges,
2029  bool add_child = true) const {
2030  moab::Interface &moab = mField.get_moab();
2032  const EntityHandle conn[] = {father, mother};
2033  Range vert_tets;
2034  CHKERR moab.get_adjacencies(conn, 2, 3, false, vert_tets,
2035  moab::Interface::UNION);
2036  vert_tets = intersect(vert_tets, proc_tets);
2037  Range out_tets;
2038  CHKERR nodeMergerPtr->mergeNodes(father, mother, out_tets, &vert_tets,
2039  onlyIfImproveQuality, 0, lineSearch, tH);
2040  out_tets.merge(subtract(proc_tets, vert_tets));
2041  proc_tets.swap(out_tets);
2042 
2043  if (add_child && nodeMergerPtr->getSucessMerge()) {
2044 
2045  auto &parent_child_map = nodeMergerPtr->getParentChildMap();
2046 
2047  Range child_ents;
2048  for (auto &it : parent_child_map)
2049  child_ents.insert(it.pArent);
2050 
2051  Range new_surf_child_ents = intersect(new_surf, child_ents);
2052  new_surf = subtract(new_surf, new_surf_child_ents);
2053  Range child_surf_ents;
2054  CHKERR updateRangeByChilds(parent_child_map, new_surf_child_ents,
2055  child_surf_ents);
2056  new_surf.merge(child_surf_ents);
2057 
2058  Range edges_to_merge_child_ents = intersect(edges_to_merge, child_ents);
2059  edges_to_merge = subtract(edges_to_merge, edges_to_merge_child_ents);
2060  Range merged_child_edge_ents;
2061  CHKERR updateRangeByChilds(parent_child_map, edges_to_merge_child_ents,
2062  merged_child_edge_ents);
2063 
2064  Range not_merged_edges_child_ents =
2065  intersect(not_merged_edges, child_ents);
2066  not_merged_edges =
2067  subtract(not_merged_edges, not_merged_edges_child_ents);
2068  Range not_merged_child_edge_ents;
2069  CHKERR updateRangeByChilds(parent_child_map,
2070  not_merged_edges_child_ents,
2071  not_merged_child_edge_ents);
2072 
2073  merged_child_edge_ents =
2074  subtract(merged_child_edge_ents, not_merged_child_edge_ents);
2075  edges_to_merge.merge(merged_child_edge_ents);
2076  not_merged_edges.merge(not_merged_child_edge_ents);
2077 
2078  if (updateMehsets) {
2080  (*mField.getInterface<MeshsetsManager>()), cubit_it)) {
2081  EntityHandle cubit_meshset = cubit_it->meshset;
2082  Range parent_ents;
2083  CHKERR moab.get_entities_by_handle(cubit_meshset, parent_ents,
2084  true);
2085  Range child_ents;
2086  CHKERR updateRangeByChilds(parent_child_map, parent_ents,
2087  child_ents);
2088  CHKERR moab.add_entities(cubit_meshset, child_ents);
2089  }
2090  }
2091  }
2093  }
2094 
2095  private:
2096  MoFEMErrorCode updateRangeByChilds(
2097  const NodeMergerInterface::ParentChildMap &parent_child_map,
2098  const Range &parents, Range &childs) const {
2100  NodeMergerInterface::ParentChildMap::nth_index<0>::type::iterator it;
2101  for (Range::const_iterator eit = parents.begin(); eit != parents.end();
2102  eit++) {
2103  it = parent_child_map.get<0>().find(*eit);
2104  if (it == parent_child_map.get<0>().end())
2105  continue;
2106  childs.insert(it->cHild);
2107  }
2109  }
2110  };
2111 
2112  /**
2113  * \brief Calculate edge length
2114  */
2115  struct LengthMap {
2116  Tag tH;
2117  CoreInterface &mField;
2118  moab::Interface &moab;
2119  const double maxLength;
2120  LengthMap(CoreInterface &m_field, Tag th, double max_length)
2121  : tH(th), mField(m_field), moab(m_field.get_moab()),
2122  maxLength(max_length) {
2123  gammaL = 1.;
2124  gammaQ = 3.;
2125  acceptedThrasholdMergeQuality = 0.5;
2126  }
2127 
2128  double
2129  gammaL; ///< Controls importance of length when ranking edges for merge
2130  double
2131  gammaQ; ///< Controls importance of quality when ranking edges for merge
2132  double acceptedThrasholdMergeQuality; ///< Do not merge quality if quality
2133  ///< above accepted thrashold
2134 
2135  MoFEMErrorCode operator()(const Range &tets, const Range &edges,
2136  LengthMapData_multi_index &length_map,
2137  double &ave) const {
2138  int num_nodes;
2139  const EntityHandle *conn;
2140  double coords[6];
2142  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
2143  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
2144  VectorDouble3 delta(3);
2145  for (auto edge : edges) {
2146  CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
2147  Range adj_tet;
2148  CHKERR moab.get_adjacencies(conn, num_nodes, 3, false, adj_tet);
2149  adj_tet = intersect(adj_tet, tets);
2150  if (tH) {
2151  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
2152  } else {
2153  CHKERR moab.get_coords(conn, num_nodes, coords);
2154  }
2155  double q = 1;
2156  CHKERR mField.getInterface<Tools>()->minTetsQuality(adj_tet, q, tH);
2157  if (q < acceptedThrasholdMergeQuality) {
2158  noalias(delta) = (s0 - s1) / maxLength;
2159  double dot = inner_prod(delta, delta);
2160  double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
2161  length_map.insert(LengthMapData(val, q, edge));
2162  }
2163  }
2164  ave = 0;
2165  for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
2166  length_map.get<0>().begin();
2167  mit != length_map.get<0>().end(); mit++) {
2168  ave += mit->qUality;
2169  }
2170  ave /= length_map.size();
2172  }
2173  };
2174 
2175  /**
2176  * \brief Topological relations
2177  */
2178  struct Toplogy {
2179 
2180  CoreInterface &mField;
2181  Tag tH;
2182  const double tOL;
2183  Toplogy(CoreInterface &m_field, Tag th, const double tol)
2184  : mField(m_field), tH(th), tOL(tol) {}
2185 
2186  enum TYPE {
2187  FREE = 0,
2188  SKIN = 1 << 0,
2189  SURFACE = 1 << 1,
2190  SURFACE_SKIN = 1 << 2,
2191  FRONT_ENDS = 1 << 3,
2192  FIX_EDGES = 1 << 4,
2193  FIX_CORNERS = 1 << 5
2194  };
2195 
2196  typedef map<int, Range> SetsMap;
2197 
2198  MoFEMErrorCode classifyVerts(const Range &surface, const Range &tets,
2199  const Range &fixed_edges,
2200  const Range &corner_nodes,
2201  SetsMap &sets_map) const {
2202  moab::Interface &moab(mField.get_moab());
2203  Skinner skin(&moab);
2205 
2206  sets_map[FIX_CORNERS].merge(corner_nodes);
2207  Range fixed_verts;
2208  CHKERR moab.get_connectivity(fixed_edges, fixed_verts, true);
2209  sets_map[FIX_EDGES].swap(fixed_verts);
2210 
2211  Range tets_skin;
2212  CHKERR skin.find_skin(0, tets, false, tets_skin);
2213  Range tets_skin_edges;
2214  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2215  moab::Interface::UNION);
2216 
2217  // surface skin
2218  Range surface_skin;
2219  CHKERR skin.find_skin(0, surface, false, surface_skin);
2220  Range front_in_the_body;
2221  front_in_the_body = subtract(surface_skin, tets_skin_edges);
2222  Range front_ends;
2223  CHKERR skin.find_skin(0, front_in_the_body, false, front_ends);
2224  sets_map[FRONT_ENDS].swap(front_ends);
2225 
2226  Range surface_skin_verts;
2227  CHKERR moab.get_connectivity(surface_skin, surface_skin_verts, true);
2228  sets_map[SURFACE_SKIN].swap(surface_skin_verts);
2229 
2230  // surface
2231  Range surface_verts;
2232  CHKERR moab.get_connectivity(surface, surface_verts, true);
2233  sets_map[SURFACE].swap(surface_verts);
2234 
2235  // skin
2236  Range tets_skin_verts;
2237  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
2238  sets_map[SKIN].swap(tets_skin_verts);
2239 
2240  Range tets_verts;
2241  CHKERR moab.get_connectivity(tets, tets_verts, true);
2242  sets_map[FREE].swap(tets_verts);
2243 
2245  }
2246 
2247  MoFEMErrorCode getProcTets(const Range &tets, const Range &edges_to_merge,
2248  Range &proc_tets) const {
2249  moab::Interface &moab(mField.get_moab());
2251  Range edges_to_merge_verts;
2252  CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts, true);
2253  Range edges_to_merge_verts_tets;
2254  CHKERR moab.get_adjacencies(edges_to_merge_verts, 3, false,
2255  edges_to_merge_verts_tets,
2256  moab::Interface::UNION);
2257  edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
2258  proc_tets.swap(edges_to_merge_verts_tets);
2260  }
2261 
2262  MoFEMErrorCode edgesToMerge(const Range &surface, const Range &tets,
2263  Range &edges_to_merge) const {
2264  moab::Interface &moab(mField.get_moab());
2266  Range surface_verts;
2267  CHKERR moab.get_connectivity(surface, surface_verts, true);
2268  Range surface_verts_edges;
2269  CHKERR moab.get_adjacencies(surface_verts, 1, false, surface_verts_edges,
2270  moab::Interface::UNION);
2271  edges_to_merge.merge(surface_verts_edges);
2272  Range tets_edges;
2273  CHKERR moab.get_adjacencies(tets, 1, false, tets_edges,
2274  moab::Interface::UNION);
2275  edges_to_merge = intersect(edges_to_merge, tets_edges);
2277  }
2278 
2279  MoFEMErrorCode removeBadEdges(const Range &surface, const Range &tets,
2280  const Range &fixed_edges,
2281  const Range &corner_nodes,
2282  Range &edges_to_merge,
2283  Range &not_merged_edges) {
2284  moab::Interface &moab(mField.get_moab());
2286 
2287  // find skin
2288  Skinner skin(&moab);
2289  Range tets_skin;
2290  CHKERR skin.find_skin(0, tets, false, tets_skin);
2291  Range surface_skin;
2292  CHKERR skin.find_skin(0, surface, false, surface_skin);
2293 
2294  // end nodes
2295  Range tets_skin_edges;
2296  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2297  moab::Interface::UNION);
2298  Range surface_front;
2299  surface_front = subtract(surface_skin, tets_skin_edges);
2300  Range surface_front_nodes;
2301  CHKERR moab.get_connectivity(surface_front, surface_front_nodes, true);
2302  Range ends_nodes;
2303  CHKERR skin.find_skin(0, surface_front, false, ends_nodes);
2304 
2305  // remove bad merges
2306 
2307  // get surface and body skin verts
2308  Range surface_edges;
2309  CHKERR moab.get_adjacencies(surface, 1, false, surface_edges,
2310  moab::Interface::UNION);
2311  // get nodes on the surface
2312  Range surface_edges_verts;
2313  CHKERR moab.get_connectivity(surface_edges, surface_edges_verts, true);
2314  // get vertices on the body skin
2315  Range tets_skin_edges_verts;
2316  CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts,
2317  true);
2318 
2319  Range edges_to_remove;
2320 
2321  // remove edges self connected to body skin
2322  {
2323  Range ents_nodes_and_edges;
2324  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2325  ents_nodes_and_edges.merge(tets_skin_edges);
2326  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2327  0, false);
2328  }
2329  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2330  not_merged_edges.merge(edges_to_remove);
2331 
2332  // remove edges self connected to surface
2333  {
2334  Range ents_nodes_and_edges;
2335  ents_nodes_and_edges.merge(surface_edges_verts);
2336  ents_nodes_and_edges.merge(surface_edges);
2337  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2338  ents_nodes_and_edges.merge(tets_skin_edges);
2339  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2340  0, false);
2341  }
2342  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2343  not_merged_edges.merge(edges_to_remove);
2344 
2345  // remove edges adjacent corner_nodes execpt those on fixed edges
2346  Range fixed_edges_nodes;
2347  CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes, true);
2348  {
2349  Range ents_nodes_and_edges;
2350  ents_nodes_and_edges.merge(fixed_edges_nodes);
2351  ents_nodes_and_edges.merge(ends_nodes);
2352  ents_nodes_and_edges.merge(corner_nodes);
2353  ents_nodes_and_edges.merge(fixed_edges);
2354  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2355  0, false);
2356  }
2357  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2358  not_merged_edges.merge(edges_to_remove);
2359 
2360  // remove edges self connected to surface
2361  CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove, 0, false);
2362  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2363  not_merged_edges.merge(edges_to_remove);
2364 
2365  // remove edges self contented on surface skin
2366  {
2367  Range ents_nodes_and_edges;
2368  ents_nodes_and_edges.merge(surface_skin);
2369  ents_nodes_and_edges.merge(fixed_edges_nodes);
2370  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2371  0, false);
2372  }
2373  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2374  not_merged_edges.merge(edges_to_remove);
2375 
2376  // remove crack front nodes connected to the surface, except those short
2377  {
2378  Range ents_nodes_and_edges;
2379  ents_nodes_and_edges.merge(surface_front_nodes);
2380  ents_nodes_and_edges.merge(surface_front);
2381  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2382  ents_nodes_and_edges.merge(tets_skin_edges);
2383  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2384  tOL, false);
2385  }
2386  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2387  not_merged_edges.merge(edges_to_remove);
2388 
2389  // remove edges connecting crack front and fixed edges, except those short
2390  {
2391  Range ents_nodes_and_edges;
2392  ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
2393  ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2394  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2395  tOL, false);
2396  }
2397  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2398  not_merged_edges.merge(edges_to_remove);
2399 
2401  }
2402 
2403  private:
2404  MoFEMErrorCode removeSelfConectingEdges(const Range &ents,
2405  Range &edges_to_remove,
2406  const bool length,
2407  bool debug) const {
2408  moab::Interface &moab(mField.get_moab());
2410  // get nodes
2411  Range ents_nodes = ents.subset_by_type(MBVERTEX);
2412  if (ents_nodes.empty()) {
2413  CHKERR moab.get_connectivity(ents, ents_nodes, true);
2414  }
2415  // edges adj. to nodes
2416  Range ents_nodes_edges;
2417  CHKERR moab.get_adjacencies(ents_nodes, 1, false, ents_nodes_edges,
2418  moab::Interface::UNION);
2419  // nodes of adj. edges
2420  Range ents_nodes_edges_nodes;
2421  CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2422  true);
2423  // hanging nodes
2424  ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2425  Range ents_nodes_edges_nodes_edges;
2426  CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1, false,
2427  ents_nodes_edges_nodes_edges,
2428  moab::Interface::UNION);
2429  // remove edges adj. to hanging edges
2430  ents_nodes_edges =
2431  subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2432  ents_nodes_edges =
2433  subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2434  if (length > 0) {
2435  Range::iterator eit = ents_nodes_edges.begin();
2436  for (; eit != ents_nodes_edges.end();) {
2437 
2438  int num_nodes;
2439  const EntityHandle *conn;
2440  rval = moab.get_connectivity(*eit, conn, num_nodes, true);
2441  double coords[6];
2442  if (tH)
2443  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
2444  else
2445  CHKERR moab.get_coords(conn, num_nodes, coords);
2446 
2447  auto get_edge_length = [coords]() {
2449  &coords[0], &coords[1], &coords[2]);
2452  t_delta(i) = t_coords(i);
2453  ++t_coords;
2454  t_delta(i) -= t_coords(i);
2455  return sqrt(t_delta(i) * t_delta(i));
2456  };
2457 
2458  if (get_edge_length() < tOL) {
2459  eit = ents_nodes_edges.erase(eit);
2460  } else {
2461  eit++;
2462  }
2463  }
2464  }
2465  edges_to_remove.swap(ents_nodes_edges);
2466  if (debug) {
2467  CHKERR SaveData(moab)("edges_to_remove.vtk", edges_to_remove);
2468  }
2470  }
2471  };
2472 
2473  Range not_merged_edges;
2474  const double tol = 0.1;
2475  CHKERR Toplogy(m_field, th, tol * aveLength)
2476  .edgesToMerge(surface, tets, edges_to_merge);
2477  CHKERR Toplogy(m_field, th, tol * aveLength)
2478  .removeBadEdges(surface, tets, fixed_edges, corner_nodes, edges_to_merge,
2479  not_merged_edges);
2480  Toplogy::SetsMap sets_map;
2481  CHKERR Toplogy(m_field, th, tol * aveLength)
2482  .classifyVerts(surface, tets, fixed_edges, corner_nodes, sets_map);
2483  Range proc_tets;
2484  CHKERR Toplogy(m_field, th, tol * aveLength)
2485  .getProcTets(tets, edges_to_merge, proc_tets);
2486  out_tets = subtract(tets, proc_tets);
2487 
2488  if (bit_ptr) {
2489  Range all_out_ents = out_tets;
2490  for (int dd = 2; dd >= 0; dd--) {
2491  CHKERR moab.get_adjacencies(out_tets, dd, false, all_out_ents,
2492  moab::Interface::UNION);
2493  }
2494  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(all_out_ents,
2495  *bit_ptr);
2496  }
2497 
2498  int nb_nodes_merged = 0;
2499  LengthMapData_multi_index length_map;
2500  new_surf = surface;
2501 
2502  double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2503  for (int pp = 0; pp != nbMaxMergingCycles; ++pp) {
2504 
2505  int nb_nodes_merged_p = nb_nodes_merged;
2506  length_map.clear();
2507  min_pp = min_p;
2508  min_p = min;
2509  CHKERR LengthMap(m_field, th, aveLength)(proc_tets, edges_to_merge,
2510  length_map, ave);
2511  min = length_map.get<2>().begin()->qUality;
2512  if (pp == 0) {
2513  ave0 = ave;
2514  }
2515 
2516  int nn = 0;
2517  Range collapsed_edges;
2518  for (auto mit = length_map.get<0>().begin();
2519  mit != length_map.get<0>().end(); mit++, nn++) {
2520 
2521  if (!mit->skip) {
2522 
2523  auto get_father_and_mother =
2524  [&](EntityHandle &father, EntityHandle &mother, int &line_search) {
2526  int num_nodes;
2527  const EntityHandle *conn;
2528  CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes, true);
2529  int conn_type[2] = {0, 0};
2530  for (int nn = 0; nn != 2; nn++) {
2531  conn_type[nn] = 0;
2532  for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2533  sit != sets_map.rend(); sit++) {
2534  if (sit->second.find(conn[nn]) != sit->second.end()) {
2535  conn_type[nn] |= sit->first;
2536  }
2537  }
2538  }
2539  int type_father, type_mother;
2540  if (conn_type[0] > conn_type[1]) {
2541  father = conn[0];
2542  mother = conn[1];
2543  type_father = conn_type[0];
2544  type_mother = conn_type[1];
2545  } else {
2546  father = conn[1];
2547  mother = conn[0];
2548  type_father = conn_type[1];
2549  type_mother = conn_type[0];
2550  }
2551  if (type_father == type_mother) {
2552  line_search = lineSearchSteps;
2553  }
2555  };
2556 
2557  int line_search = 0;
2558  EntityHandle father, mother;
2559  CHKERR get_father_and_mother(father, mother, line_search);
2560  CHKERR MergeNodes(m_field, true, line_search, th,
2561  update_meshsets)(father, mother, proc_tets, new_surf,
2562  edges_to_merge, not_merged_edges);
2563 
2564  if (m_field.getInterface<NodeMergerInterface>()->getSucessMerge()) {
2565  Range adj_mother_tets;
2566  CHKERR moab.get_adjacencies(&mother, 1, 3, false, adj_mother_tets);
2567  Range adj_mother_tets_nodes;
2568  CHKERR moab.get_connectivity(adj_mother_tets, adj_mother_tets_nodes,
2569  true);
2570  Range adj_edges;
2571  CHKERR moab.get_adjacencies(adj_mother_tets_nodes, 1, false,
2572  adj_edges, moab::Interface::UNION);
2573  CHKERR moab.get_adjacencies(&father, 1, 1, false, adj_edges,
2574  moab::Interface::UNION);
2575  Range proc_edges;
2576  CHKERR moab.get_adjacencies(proc_tets, 1, true, proc_edges);
2577  adj_edges = intersect(proc_edges, adj_edges);
2578  for (Range::iterator ait = adj_edges.begin(); ait != adj_edges.end();
2579  ++ait) {
2580  auto miit = length_map.get<1>().find(*ait);
2581  if (miit != length_map.get<1>().end())
2582  (const_cast<LengthMapData &>(*miit)).skip = true;
2583  }
2584  nb_nodes_merged++;
2585  collapsed_edges.insert(mit->eDge);
2586  }
2587 
2588  if (nn > static_cast<int>(length_map.size() / fraction_level))
2589  break;
2590  if (mit->qUality > ave)
2591  break;
2592  }
2593  }
2594 
2595  Range adj_faces, adj_edges;
2596  CHKERR moab.get_adjacencies(proc_tets, 2, false, adj_faces,
2597  moab::Interface::UNION);
2598  new_surf = intersect(new_surf, adj_faces);
2599 
2600  PetscPrintf(m_field.get_comm(),
2601  "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e\n", pp,
2602  nb_nodes_merged, ave, min);
2603 
2604  if (debug) {
2605  // current surface and merged edges in step
2606  std::string name;
2607  name = "node_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2608  CHKERR SaveData(moab)(name, unite(new_surf, collapsed_edges));
2609  name = "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) +
2610  ".vtk";
2611  CHKERR SaveData(moab)(name, unite(new_surf, edges_to_merge));
2612  }
2613 
2614  if (nb_nodes_merged == nb_nodes_merged_p)
2615  break;
2616  if (min > 1e-2 && min == min_pp)
2617  break;
2618  if (min > ave0)
2619  break;
2620 
2621  CHKERR moab.get_adjacencies(proc_tets, 1, false, adj_edges,
2622  moab::Interface::UNION);
2623  edges_to_merge = intersect(edges_to_merge, adj_edges);
2624  CHKERR Toplogy(m_field, th, tol * aveLength)
2625  .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
2626  edges_to_merge, not_merged_edges);
2627  }
2628 
2629  if (bit_ptr) {
2630  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(proc_tets,
2631  *bit_ptr);
2632  }
2633  out_tets.merge(proc_tets);
2634 
2636 }
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
double maxLength
Maximal edge length.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
double aveLength
Average edge length.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
double tol
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
static const bool debug
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:99
#define CHKERR
Inline error check.
Definition: definitions.h:594
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
multi_index_container< LengthMapData, indexed_by< ordered_non_unique< member< LengthMapData, double, &LengthMapData::lEngth > >, hashed_unique< member< LengthMapData, EntityHandle, &LengthMapData::eDge > >, ordered_non_unique< member< LengthMapData, double, &LengthMapData::qUality > > > > LengthMapData_multi_index
multi_index_container< ParentChild, indexed_by< hashed_unique< member< ParentChild, EntityHandle, &ParentChild::pArent > >, hashed_non_unique< member< ParentChild, EntityHandle, &ParentChild::cHild > > > > ParentChildMap
Definition: NodeMerger.hpp:128

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

2642  {
2643  CoreInterface &m_field = cOre;
2645  Range tets_level;
2646  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2647  trim_bit, BitRefLevel().set(), MBTET, tets_level);
2648 
2649  Range edges_to_merge;
2650  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByParentType(
2651  trim_bit, cut_bit | trim_bit, MBEDGE, edges_to_merge);
2652 
2653  // get all entities not in database
2654  Range all_ents_not_in_database_before;
2655  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2656  all_ents_not_in_database_before);
2657 
2658  Range out_new_tets, out_new_surf;
2659  CHKERR mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
2660  corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
2661  th, update_meshsets, &bit, debug);
2662 
2663  // get all entities not in database after merge
2664  Range all_ents_not_in_database_after;
2665  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2666  all_ents_not_in_database_after);
2667 
2668  // delete hanging entities
2669  all_ents_not_in_database_after =
2670  subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
2671  Range meshsets;
2672  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
2673  false);
2674  for (auto m : meshsets)
2675  CHKERR m_field.get_moab().remove_entities(m,
2676  all_ents_not_in_database_after);
2677 
2678  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
2679 
2680  mergedVolumes.swap(out_new_tets);
2681  mergedSurfaces.swap(out_new_surf);
2683 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
MoFEMErrorCode mergeBadEdges(const int fraction_level, const Range &tets, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Range &merged_nodes, Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets=false, const BitRefLevel *bit_ptr=NULL, const bool debug=false)
Merge edges.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:594
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ mergeSurface()

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

merge surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 142 of file CutMeshInterface.cpp.

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

◆ mergeVolumes()

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

merge volume entities

Parameters
volumeentities which going to be added
Returns
error code

Definition at line 148 of file CutMeshInterface.cpp.

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

◆ moveMidNodesOnCutEdges()

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

projecting of mid edge nodes on new mesh on surface

Returns
error code

Definition at line 1341 of file CutMeshInterface.cpp.

1341  {
1343 
1344  CoreInterface &m_field = cOre;
1345  moab::Interface &moab = m_field.get_moab();
1347 
1348  // Range out_side_vertices;
1349  for (auto m : verticesOnCutEdges) {
1350  double dist = m.second.dIst;
1351  VectorDouble3 new_coors = m.second.rayPoint + dist * m.second.unitRayDir;
1352  if (th)
1353  CHKERR moab.tag_set_data(th, &m.first, 1, &new_coors[0]);
1354  else
1355  CHKERR moab.set_coords(&m.first, 1, &new_coors[0]);
1356  }
1357 
1359 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ moveMidNodesOnTrimmedEdges()

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

move trimmed edges mid nodes

Returns
error code

Definition at line 1361 of file CutMeshInterface.cpp.

1361  {
1362  CoreInterface &m_field = cOre;
1363  moab::Interface &moab = m_field.get_moab();
1365  // Range out_side_vertices;
1366  for (auto &v : verticesOnTrimEdges) {
1367  double dist = v.second.dIst;
1368  VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1369  if (th) {
1370  CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1371  } else {
1372  CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1373  }
1374  }
1376 }
map< EntityHandle, TreeData > verticesOnTrimEdges
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ projectZeroDistanceEnts() [1/2]

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

Definition at line 928 of file CutMeshInterface.cpp.

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

◆ projectZeroDistanceEnts() [2/2]

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

Definition at line 329 of file CutMeshInterface.hpp.

329  {
330  return zeroDistanceEnts;
331  }

◆ query_interface()

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

Implements MoFEM::UnknownInterface.

Definition at line 40 of file CutMeshInterface.cpp.

41  {
43  *iface = NULL;
44  if (uuid == IDD_MOFEMCutMesh) {
45  *iface = const_cast<CutMeshInterface *>(this);
47  }
48  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
50 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506
CutMeshInterface(const MoFEM::Core &core)
static const MOFEMuuid IDD_MOFEMCutMesh

◆ rebuildMeshWithTetGen()

MoFEMErrorCode MoFEM::CutMeshInterface::rebuildMeshWithTetGen ( vector< string > &  switches,
const BitRefLevel mesh_bit,
const BitRefLevel bit,
const Range &  surface,
const Range &  fixed_edges,
const Range &  corner_nodes,
Tag  th = NULL,
const bool  debug = false 
)
Examples:
mesh_cut.cpp.

Definition at line 2687 of file CutMeshInterface.cpp.

2690  {
2691  CoreInterface &m_field = cOre;
2692  moab::Interface &moab = m_field.get_moab();
2693  TetGenInterface *tetgen_iface;
2695  CHKERR m_field.getInterface(tetgen_iface);
2696 
2697  tetGenData.clear();
2698  moabTetGenMap.clear();
2699  tetGenMoabMap.clear();
2700 
2701  if (tetGenData.size() < 1) {
2702  tetGenData.push_back(new tetgenio);
2703  }
2704  tetgenio &in = tetGenData.back();
2705 
2706  struct BitEnts {
2707 
2708  CoreInterface &mField;
2709  const BitRefLevel &bIt;
2710  BitEnts(CoreInterface &m_field, const BitRefLevel &bit)
2711  : mField(m_field), bIt(bit) {}
2712 
2713  Range mTets;
2714  Range mTris;
2715  Range mEdges;
2716  Range mNodes;
2717 
2718  MoFEMErrorCode getLevelEnts() {
2720  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2721  bIt, BitRefLevel().set(), MBTET, mTets);
2722  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2723  bIt, BitRefLevel().set(), MBTRI, mTris);
2724  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2725  bIt, BitRefLevel().set(), MBEDGE, mEdges);
2726  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2727  bIt, BitRefLevel().set(), MBVERTEX, mNodes);
2729  }
2730 
2731  Range mSkin;
2732  Range mSkinNodes;
2733  Range mSkinEdges;
2734 
2735  MoFEMErrorCode getSkin() {
2736  moab::Interface &moab = mField.get_moab();
2738  Skinner skin(&moab);
2739  CHKERR skin.find_skin(0, mTets, false, mSkin);
2740  CHKERR mField.get_moab().get_connectivity(mSkin, mSkinNodes, true);
2741  CHKERR mField.get_moab().get_adjacencies(mSkin, 1, false, mSkinEdges,
2742  moab::Interface::UNION);
2744  }
2745  };
2746 
2747  struct SurfaceEnts {
2748 
2749  CoreInterface &mField;
2750  SurfaceEnts(CoreInterface &m_field) : mField(m_field) {}
2751 
2752  Range sNodes;
2753  Range sEdges;
2754  Range sVols;
2755  Range vNodes;
2756 
2757  MoFEMErrorCode getVolume(const BitEnts &bit_ents, const Range &tris) {
2758  moab::Interface &moab = mField.get_moab();
2760  CHKERR moab.get_connectivity(tris, sNodes, true);
2761  CHKERR moab.get_adjacencies(tris, 1, false, sEdges,
2762  moab::Interface::UNION);
2763  CHKERR moab.get_adjacencies(sNodes, 3, false, sVols,
2764  moab::Interface::UNION);
2765  sVols = intersect(sVols, bit_ents.mTets);
2766  CHKERR moab.get_connectivity(sVols, vNodes, true);
2768  }
2769 
2770  Range sSkin;
2771  Range sSkinNodes;
2772  Range vSkin;
2773  Range vSkinNodes;
2774  Range vSkinWithoutBodySkin;
2775  Range vSkinNodesWithoutBodySkin;
2776  Range vSkinOnBodySkin;
2777  Range vSkinOnBodySkinNodes;
2778 
2779  MoFEMErrorCode getSkin(const BitEnts &bit_ents, const Range &tris,
2780  const int levels = 3) {
2781  moab::Interface &moab = mField.get_moab();
2783  Skinner skin(&moab);
2784  rval = skin.find_skin(0, sVols, false, vSkin);
2785  for (int ll = 0; ll != levels; ll++) {
2786  CHKERR moab.get_adjacencies(vSkin, 3, false, sVols,
2787  moab::Interface::UNION);
2788  sVols = intersect(sVols, bit_ents.mTets);
2789  vSkin.clear();
2790  CHKERR skin.find_skin(0, sVols, false, vSkin);
2791  }
2792  vSkinWithoutBodySkin = subtract(vSkin, bit_ents.mSkin);
2793  vSkinOnBodySkin = intersect(vSkin, bit_ents.mSkin);
2794  CHKERR moab.get_connectivity(vSkinOnBodySkin, vSkinOnBodySkinNodes, true);
2795  CHKERR moab.get_connectivity(sVols, vNodes, true);
2796  CHKERR moab.get_connectivity(vSkin, vSkinNodes, true);
2797  vSkinNodesWithoutBodySkin = subtract(vSkinNodes, bit_ents.mSkinNodes);
2798  CHKERR skin.find_skin(0, tris, false, sSkin);
2799  CHKERR moab.get_connectivity(sSkin, sSkinNodes, true);
2800  tVols = sVols;
2802  }
2803 
2804  Range tVols;
2805 
2806  MoFEMErrorCode getTetsForRemesh(const BitEnts &bit_ents, Tag th = NULL) {
2807  moab::Interface &moab = mField.get_moab();
2809 
2810  Range tets_with_four_nodes_on_skin;
2811  CHKERR moab.get_adjacencies(vSkinOnBodySkinNodes, 3, false,
2812  tets_with_four_nodes_on_skin,
2813  moab::Interface::UNION);
2814  Range tets_nodes;
2815  CHKERR moab.get_connectivity(tets_with_four_nodes_on_skin, tets_nodes,
2816  true);
2817  tets_nodes = subtract(tets_nodes, vSkinOnBodySkinNodes);
2818  Range other_tets;
2819  CHKERR moab.get_adjacencies(tets_nodes, 3, false, other_tets,
2820  moab::Interface::UNION);
2821  tets_with_four_nodes_on_skin =
2822  subtract(tets_with_four_nodes_on_skin, other_tets);
2823  Range to_remove;
2824  for (Range::iterator tit = tets_with_four_nodes_on_skin.begin();
2825  tit != tets_with_four_nodes_on_skin.end(); tit++) {
2826  int num_nodes;
2827  const EntityHandle *conn;
2828  CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
2829  double coords[12];
2830  if (th) {
2831  CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
2832  } else {
2833  CHKERR moab.get_coords(conn, num_nodes, coords);
2834  }
2835  double quality = Tools::volumeLengthQuality(coords);
2836  if (quality < 1e-2) {
2837  to_remove.insert(*tit);
2838  }
2839  }
2840 
2841  sVols = subtract(sVols, to_remove);
2842 
2843  Skinner skin(&moab);
2844  vSkin.clear();
2845  CHKERR skin.find_skin(0, sVols, false, vSkin);
2846  Range m_skin;
2847  CHKERR
2848  skin.find_skin(0, subtract(bit_ents.mSkin, to_remove), false, m_skin);
2849  vSkinWithoutBodySkin = subtract(vSkin, m_skin);
2850  vSkinOnBodySkin = intersect(vSkin, m_skin);
2851  vNodes.clear();
2852  vSkinNodes.clear();
2853  vSkinOnBodySkinNodes.clear();
2854  CHKERR moab.get_connectivity(sVols, vNodes, true);
2855  CHKERR moab.get_connectivity(vSkinOnBodySkin, vSkinOnBodySkinNodes, true);
2856  CHKERR moab.get_connectivity(vSkin, vSkinNodes, true);
2858  }
2859  };
2860 
2861  BitEnts bit_ents(m_field, mesh_bit);
2862  CHKERR bit_ents.getLevelEnts();
2863  CHKERR bit_ents.getSkin();
2864  SurfaceEnts surf_ents(m_field);
2865  CHKERR surf_ents.getVolume(bit_ents, surface);
2866  CHKERR surf_ents.getSkin(bit_ents, surface);
2867  CHKERR surf_ents.getTetsForRemesh(bit_ents);
2868 
2869  map<int, Range> types_ents;
2870  types_ents[TetGenInterface::RIDGEVERTEX].merge(
2871  surf_ents.vSkinNodesWithoutBodySkin);
2872  // FREESEGVERTEX
2873  types_ents[TetGenInterface::FREESEGVERTEX].merge(surf_ents.sSkinNodes);
2874  types_ents[TetGenInterface::FREESEGVERTEX] =
2875  subtract(types_ents[TetGenInterface::FREESEGVERTEX],
2876  types_ents[TetGenInterface::RIDGEVERTEX]);
2877  // FREEFACETVERTEX
2878  types_ents[TetGenInterface::FREEFACETVERTEX].merge(surf_ents.sNodes);
2879  types_ents[TetGenInterface::FREEFACETVERTEX] =
2880  subtract(types_ents[TetGenInterface::FREEFACETVERTEX],
2881  types_ents[TetGenInterface::RIDGEVERTEX]);
2882  types_ents[TetGenInterface::FREEFACETVERTEX] =
2883  subtract(types_ents[TetGenInterface::FREEFACETVERTEX],
2884  types_ents[TetGenInterface::FREESEGVERTEX]);
2885  // FREEVOLVERTEX
2886  types_ents[TetGenInterface::FREEVOLVERTEX].merge(surf_ents.vNodes);
2887  types_ents[TetGenInterface::FREEVOLVERTEX] =
2888  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2889  types_ents[TetGenInterface::RIDGEVERTEX]);
2890  types_ents[TetGenInterface::FREEVOLVERTEX] =
2891  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2892  types_ents[TetGenInterface::FREESEGVERTEX]);
2893  types_ents[TetGenInterface::FREEVOLVERTEX] =
2894  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2895  types_ents[TetGenInterface::FREEFACETVERTEX]);
2896 
2897  Tag th_marker;
2898  // Clean markers
2899  rval = m_field.get_moab().tag_get_handle("TETGEN_MARKER", th_marker);
2900  if (rval == MB_SUCCESS) {
2901  CHKERR m_field.get_moab().tag_delete(th_marker);
2902  rval = MB_SUCCESS;
2903  }
2904 
2905  int def_marker = 0;
2906  CHKERR m_field.get_moab().tag_get_handle(
2907  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
2908  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
2909 
2910  // Mark surface with id = 1
2911  vector<int> markers(surface.size(), 1);
2912  CHKERR m_field.get_moab().tag_set_data(th_marker, surface, &*markers.begin());
2913  // Mark all side sets
2914  int shift = 1;
2915  map<int, int> id_shift_map; // each meshset has set unique bit
2917  (*cOre.getInterface<MeshsetsManager>()), SIDESET, it)) {
2918  int ms_id = it->getMeshsetId();
2919  id_shift_map[ms_id] = 1 << shift; // shift bit
2920  ++shift;
2921  Range sideset_faces;
2922  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
2923  ms_id, SIDESET, 2, sideset_faces, true);
2924  markers.resize(sideset_faces.size());
2925  CHKERR m_field.get_moab().tag_get_data(th_marker, sideset_faces,
2926  &*markers.begin());
2927  for (unsigned int ii = 0; ii < markers.size(); ii++) {
2928  markers[ii] |= id_shift_map[ms_id]; // add bit to marker
2929  }
2930  CHKERR m_field.get_moab().tag_set_data(th_marker, sideset_faces,
2931  &*markers.begin());
2932  }
2933  Range nodes_to_remove; // none
2934  markers.resize(nodes_to_remove.size());
2935  fill(markers.begin(), markers.end(), -1);
2936  CHKERR m_field.get_moab().tag_set_data(th_marker, nodes_to_remove,
2937  &*markers.begin());
2938 
2939  // nodes
2940  if (tetGenData.size() == 1) {
2941 
2942  Range ents_to_tetgen = surf_ents.sVols;
2943  CHKERR m_field.get_moab().get_connectivity(surf_ents.sVols, ents_to_tetgen,
2944  true);
2945  for (int dd = 2; dd >= 1; dd--) {
2946  CHKERR m_field.get_moab().get_adjacencies(
2947  surf_ents.sVols, dd, false, ents_to_tetgen, moab::Interface::UNION);
2948  }
2949 
2950  // Load mesh to TetGen data structures
2951  CHKERR tetgen_iface->inData(ents_to_tetgen, in, moabTetGenMap,
2952  tetGenMoabMap, th);
2953  CHKERR tetgen_iface->setGeomData(in, moabTetGenMap, tetGenMoabMap,
2954  types_ents);
2955  std::vector<pair<Range, int>> markers;
2956  for (Range::iterator tit = surface.begin(); tit != surface.end(); tit++) {
2957  Range facet;
2958  facet.insert(*tit);
2959  markers.push_back(pair<Range, int>(facet, 2));
2960  }
2961  for (Range::iterator tit = surf_ents.vSkinWithoutBodySkin.begin();
2962  tit != surf_ents.vSkinWithoutBodySkin.end(); tit++) {
2963  Range facet;
2964  facet.insert(*tit);
2965  markers.push_back(pair<Range, int>(facet, 1));
2966  }
2967  Range other_facets;
2968  other_facets = subtract(surf_ents.vSkin, surf_ents.vSkinWithoutBodySkin);
2969  for (Range::iterator tit = other_facets.begin(); tit != other_facets.end();
2970  tit++) {
2971  Range facet;
2972  facet.insert(*tit);
2973  markers.push_back(pair<Range, int>(facet, 0));
2974  }
2975  CHKERR tetgen_iface->setFaceData(markers, in, moabTetGenMap, tetGenMoabMap);
2976  }
2977 
2978  if (debug) {
2979  if (m_field.get_comm_rank() == 0) {
2980  char tetgen_in_file_name[] = "in";
2981  in.save_nodes(tetgen_in_file_name);
2982  in.save_elements(tetgen_in_file_name);
2983  in.save_faces(tetgen_in_file_name);
2984  in.save_edges(tetgen_in_file_name);
2985  in.save_poly(tetgen_in_file_name);
2986  }
2987  }
2988 
2989  // generate new mesh
2990  {
2991  vector<string>::iterator sw = switches.begin();
2992  for (int ii = 0; sw != switches.end(); sw++, ii++) {
2993  tetgenio &_in_ = tetGenData.back();
2994  tetGenData.push_back(new tetgenio);
2995  tetgenio &_out_ = tetGenData.back();
2996  char *s = const_cast<char *>(sw->c_str());
2997  CHKERR tetgen_iface->tetRahedralize(s, _in_, _out_);
2998  }
2999  }
3000  tetgenio &out = tetGenData.back();
3001  // save elems
3002  if (debug) {
3003  char tetgen_out_file_name[] = "out";
3004  out.save_nodes(tetgen_out_file_name);
3005  out.save_elements(tetgen_out_file_name);
3006  out.save_faces(tetgen_out_file_name);
3007  out.save_edges(tetgen_out_file_name);
3008  out.save_poly(tetgen_out_file_name);
3009  }
3010 
3011  CHKERR tetgen_iface->outData(in, out, moabTetGenMap, tetGenMoabMap, bit,
3012  false, false, true, th);
3013 
3014  Range rest_of_ents = subtract(bit_ents.mTets, surf_ents.tVols);
3015  for (int dd = 2; dd >= 0; dd--) {
3016  CHKERR moab.get_adjacencies(rest_of_ents.subset_by_dimension(3), dd, false,
3017  rest_of_ents, moab::Interface::UNION);
3018  }
3019  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(rest_of_ents,
3020  bit);
3021 
3022  Range tetgen_faces;
3023  map<int, Range> face_markers_map;
3024  CHKERR tetgen_iface->getTriangleMarkers(tetGenMoabMap, out, &tetgen_faces,
3025  &face_markers_map);
3026 
3027  tetgenSurfaces = face_markers_map[1];
3028  for (map<int, Range>::iterator mit = face_markers_map.begin();
3029  mit != face_markers_map.end(); mit++) {
3031  (*cOre.getInterface<MeshsetsManager>()), SIDESET, it)) {
3032  int msId = it->getMeshsetId();
3033  if (id_shift_map[msId] & mit->first) {
3034  EntityHandle meshset = it->getMeshset();
3035  CHKERR m_field.get_moab().add_entities(
3036  meshset, mit->second.subset_by_type(MBTRI));
3037  }
3038  }
3039  }
3040 
3041  auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
3042  Range tets_level; // test at level
3043  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3044  bit, BitRefLevel().set(), MBTET, tets_level);
3045  double min_q = 1;
3046  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
3047  return min_q;
3048  };
3049 
3050  PetscPrintf(PETSC_COMM_WORLD, "Min quality tetgen %6.4g\n",
3051  get_min_quality(bit, th));
3052 
3054 }
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:32
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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
map< EntityHandle, unsigned long > moabTetGenMap
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
static const bool debug
const Range & getVolume() const
boost::ptr_vector< tetgenio > tetGenData
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field...
#define CHKERR
Inline error check.
Definition: definitions.h:594
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
map< unsigned long, EntityHandle > tetGenMoabMap

◆ refCutTrimAndMerge()

MoFEMErrorCode MoFEM::CutMeshInterface::refCutTrimAndMerge ( int first_bit,
const int  fraction_level,
const int  ref_before_cut_levels,
const int  befor_trim_levels,
Tag  th,
const double  tol_cut,
const double  tol_cut_close,
const double  tol_trim,
const double  tol_trim_close,
Range &  fixed_edges,
Range &  corner_nodes,
const bool  update_meshsets = false,
const bool  debug = false 
)
Examples:
mesh_cut.cpp.

Definition at line 459 of file CutMeshInterface.cpp.

464  {
465  CoreInterface &m_field = cOre;
467 
468  std::vector<BitRefLevel> bit_levels;
469 
470  auto get_back_bit_levels = [&]() {
471  bit_levels.push_back(BitRefLevel());
472  bit_levels.back().set(first_bit + bit_levels.size() - 1);
473  return bit_levels.back();
474  };
475 
476  Range starting_volume = vOlume;
477  Range starting_volume_nodes;
478  CHKERR m_field.get_moab().get_connectivity(starting_volume,
479  starting_volume_nodes, true);
480  std::vector<double> staring_volume_coords(3 * starting_volume_nodes.size());
481  CHKERR m_field.get_moab().get_coords(starting_volume_nodes,
482  &*staring_volume_coords.begin());
483 
484  if (th) {
485  std::vector<double> staring_volume_th_coords(3 *
486  starting_volume_nodes.size());
487  CHKERR m_field.get_moab().tag_get_data(th, starting_volume_nodes,
488  &*staring_volume_th_coords.begin());
489  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
490  &*staring_volume_th_coords.begin());
491  }
492 
493  for (int ll = 0; ll != ref_before_cut_levels; ++ll) {
494  CHKERR findEdgesToCut(&fixed_edges, &corner_nodes, tol_trim, QUIET);
495  BitRefLevel bit = get_back_bit_levels();
496  CHKERR refineBeforeCut(bit, &fixed_edges, true);
497  if (th)
498  CHKERR setTagData(th);
499  }
500 
501  if (debug && ref_before_cut_levels)
502  CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByType(
503  bit_levels.back(), BitRefLevel().set(), MBTET,
504  "out_after_refine_before_cut.vtk", "VTK", "");
505 
506  first_bit += bit_levels.size();
507  CHKERR cutTrimAndMerge(first_bit, fraction_level, befor_trim_levels, th,
508  tol_cut, tol_cut_close, tol_trim, tol_trim_close,
509  fixed_edges, corner_nodes, update_meshsets, debug);
510 
511  if (th)
512  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
513  &*staring_volume_coords.begin());
514 
516 }
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
MoFEMErrorCode refineBeforeCut(const BitRefLevel &bit, Range *fixed_edges, const bool update_meshsets, const bool debug=false)
MoFEMErrorCode findEdgesToCut(Range *fixed_edges, Range *corner_nodes, const double low_tol, int verb=QUIET, const bool debug=false)
find edges to cut
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:594
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
MoFEMErrorCode cutTrimAndMerge(int &first_bit, const int fraction_level, const int before_trim_levels, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)

◆ refineBeforeCut()

MoFEMErrorCode MoFEM::CutMeshInterface::refineBeforeCut ( const BitRefLevel bit,
Range *  fixed_edges,
const bool  update_meshsets,
const bool  debug = false 
)

Definition at line 806 of file CutMeshInterface.cpp.

809  {
810  CoreInterface &m_field = cOre;
811  moab::Interface &moab = m_field.get_moab();
812  MeshRefinement *refiner;
813  BitRefManager *bit_ref_manager;
815 
816  if (cutEdges.size() != edgesToCut.size()) {
817  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
818  }
819  CHKERR m_field.getInterface(refiner);
820  CHKERR m_field.getInterface(bit_ref_manager);
821 
822  Range ref_edges_tets;
823  CHKERR moab.get_adjacencies(cutEdges, 3, false, ref_edges_tets,
824  moab::Interface::UNION);
825  CHKERR moab.get_adjacencies(zeroDistanceVerts, 3, false, ref_edges_tets,
826  moab::Interface::UNION);
827  Range ref_edges;
828  CHKERR moab.get_adjacencies(ref_edges_tets, 1, true, ref_edges,
829  moab::Interface::UNION);
830  CHKERR refiner->add_verices_in_the_middel_of_edges(ref_edges, bit);
831  CHKERR refiner->refine_TET(vOlume, bit, false, QUIET,
832  debug ? &cutEdges : NULL);
833 
834  auto update_range = [bit_ref_manager](Range *r_ptr) {
836  if (r_ptr) {
837  Range childs;
838  CHKERR bit_ref_manager->updateRange(*r_ptr, childs);
839  r_ptr->merge(childs);
840  }
842  };
843  CHKERR update_range(fixed_edges);
844  CHKERR update_range(&vOlume);
845 
846  if (update_meshsets)
847  CHKERR UpdateMeshsets()(cOre, bit);
848 
849  Range bit_tets;
850  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
851  bit, BitRefLevel().set(), MBTET, bit_tets);
852  vOlume = intersect(vOlume, bit_tets);
853 
854  Range bit_edges;
855  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
856  bit, BitRefLevel().set(), MBEDGE, bit_edges);
857  if (fixed_edges)
858  *fixed_edges = intersect(*fixed_edges, bit_edges);
859 
861 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:594
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
map< EntityHandle, TreeData > edgesToCut

◆ refineBeforeTrim()

MoFEMErrorCode MoFEM::CutMeshInterface::refineBeforeTrim ( const BitRefLevel bit,
Range *  fixed_edges,
const bool  update_meshsets,
const bool  debug = false 
)

Definition at line 863 of file CutMeshInterface.cpp.

866  {
867  CoreInterface &m_field = cOre;
868  moab::Interface &moab = m_field.get_moab();
869  MeshRefinement *refiner;
870  BitRefManager *bit_ref_manager;
872 
873  // refine before trim
874  Range vol_trim_edges;
875  CHKERR moab.get_adjacencies(trimEdges, 3, true, vol_trim_edges,
876  moab::Interface::UNION);
877  CHKERR moab.get_adjacencies(trimNewVertices, 3, true, vol_trim_edges,
878  moab::Interface::UNION);
879  vol_trim_edges = intersect(vol_trim_edges, cutNewVolumes);
880 
881  CHKERR m_field.getInterface(refiner);
882  CHKERR m_field.getInterface(bit_ref_manager);
883 
884  Range ref_edges;
885  CHKERR moab.get_adjacencies(vol_trim_edges, 1, true, ref_edges,
886  moab::Interface::UNION);
887  CHKERR refiner->add_verices_in_the_middel_of_edges(ref_edges, bit);
888  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
889  debug ? &cutEdges : NULL);
890 
891  auto update_range = [bit_ref_manager](Range *r_ptr) {
893  if (r_ptr) {
894  Range childs;
895  CHKERR bit_ref_manager->updateRange(*r_ptr, childs);
896  r_ptr->merge(childs);
897  }
899  };
900  CHKERR update_range(fixed_edges);
901  CHKERR update_range(&cutNewVolumes);
902  CHKERR update_range(&cutNewSurfaces);
903 
904  if (update_meshsets)
905  CHKERR UpdateMeshsets()(cOre, bit);
906 
907  Range bit_tets;
908  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
909  bit, BitRefLevel().set(), MBTET, bit_tets);
910  cutNewVolumes = intersect(cutNewVolumes, bit_tets);
911 
912  Range bit_faces;
913  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
914  bit, BitRefLevel().set(), MBTRI, bit_faces);
915  cutNewSurfaces = intersect(cutNewSurfaces, bit_faces);
916  cutNewVertices.clear();
917  CHKERR moab.get_connectivity(cutNewSurfaces, cutNewVertices, true);
918 
919  Range bit_edges;
920  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
921  bit, BitRefLevel().set(), MBEDGE, bit_edges);
922  if (fixed_edges)
923  *fixed_edges = intersect(*fixed_edges, bit_edges);
924 
926 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:594
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

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

1909  {
1910  CoreInterface &m_field = cOre;
1911  moab::Interface &moab = m_field.get_moab();
1912  PrismInterface *interface;
1914  CHKERR m_field.getInterface(interface);
1915  // Remove tris on surface front
1916  {
1917  Range front_tris;
1918  EntityHandle meshset;
1919  CHKERR moab.create_meshset(MESHSET_SET, meshset);
1920  CHKERR moab.add_entities(meshset, ents);
1921  CHKERR interface->findIfTringleHasThreeNodesOnInternalSurfaceSkin(
1922  meshset, split_bit, true, front_tris);
1923  CHKERR moab.delete_entities(&meshset, 1);
1924  ents = subtract(ents, front_tris);
1925  }
1926  // Remove entities on skin
1927  Range tets;
1928  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1929  split_bit, BitRefLevel().set(), MBTET, tets);
1930  // Remove entities on skin
1931  Skinner skin(&moab);
1932  Range tets_skin;
1933  rval = skin.find_skin(0, tets, false, tets_skin);
1934  ents = subtract(ents, tets_skin);
1935 
1937 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
#define CHKERR
Inline error check.
Definition: definitions.h:594
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ saveCutEdges()

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

Definition at line 3093 of file CutMeshInterface.cpp.

3093  {
3094  CoreInterface &m_field = cOre;
3095  moab::Interface &moab = m_field.get_moab();
3097  CHKERR SaveData(moab)(prefix + "out_vol.vtk", vOlume);
3098  CHKERR SaveData(moab)(prefix + "out_surface.vtk", sUrface);
3099  CHKERR SaveData(moab)(prefix + "out_cut_edges.vtk", cutEdges);
3100  CHKERR SaveData(moab)(prefix + "out_cut_volumes.vtk", cutVolumes);
3101  CHKERR SaveData(moab)(prefix + "out_cut_new_volumes.vtk", cutNewVolumes);
3102  CHKERR SaveData(moab)(prefix + "out_cut_new_surfaces.vtk", cutNewSurfaces);
3103  CHKERR SaveData(moab)(prefix + "out_cut_zero_distance_ents.vtk",
3106 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ saveTrimEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveTrimEdges ( )

Definition at line 3108 of file CutMeshInterface.cpp.

3108  {
3109  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
3111  CHKERR SaveData(moab)("out_trim_surface.vtk", sUrface);
3112  CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
3113  CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
3114  CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
3116 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
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:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

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

3076  {
3077  CoreInterface &m_field = cOre;
3078  moab::Interface &moab = m_field.get_moab();
3080  Range nodes;
3081  if (bit.none()) {
3082  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3083  } else {
3084  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3085  bit, mask, MBVERTEX, nodes);
3086  }
3087  std::vector<double> coords(3 * nodes.size());
3088  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
3089  CHKERR moab.set_coords(nodes, &coords[0]);
3091 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ setFront()

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

Definition at line 65 of file CutMeshInterface.cpp.

65  {
67  fRont = front;
69 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506

◆ setSurface()

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

set surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 71 of file CutMeshInterface.cpp.

71  {
73  sUrface = surface;
75 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506

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

3058  {
3059  CoreInterface &m_field = cOre;
3060  moab::Interface &moab = m_field.get_moab();
3062  Range nodes;
3063  if (bit.none()) {
3064  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3065  } else {
3066  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3067  bit, BitRefLevel().set(), MBVERTEX, nodes);
3068  }
3069  std::vector<double> coords(3 * nodes.size());
3070  CHKERR moab.get_coords(nodes, &coords[0]);
3071  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
3073 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

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

136  {
138  vOlume = volume;
140 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506

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

156  {
157  CoreInterface &m_field = cOre;
158  moab::Interface &moab = m_field.get_moab();
160 
161  // Get cutting surface skin
162  Skinner skin(&moab);
163  Range surface_skin;
164  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
165 
166  CHKERR snapSurfaceToEdges(surface_skin, fixed_edges, rel_tol, abs_tol, th,
167  debug);
168 
170 }
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)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

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

176  {
177  CoreInterface &m_field = cOre;
178  moab::Interface &moab = m_field.get_moab();
181 
182  auto get_point = [i](auto &t_w, auto &t_delta, auto t) {
184  t = std::max(0., std::min(1., t));
185  t_p(i) = t_w(i) + t * t_delta(i);
186  return t_p;
187  };
188 
189  auto get_distance = [i](auto &t_p, auto &t_n) {
190  FTensor::Tensor1<double, 3> t_dist_vector;
191  t_dist_vector(i) = t_p(i) - t_n(i);
192  return sqrt(t_dist_vector(i) * t_dist_vector(i));
193  };
194 
195  map<EntityHandle, double> map_verts_length;
196 
197  for (auto f : surface_edges) {
198  int num_nodes;
199  const EntityHandle *conn_skin;
200  CHKERR moab.get_connectivity(f, conn_skin, num_nodes, true);
201  VectorDouble6 coords_skin(6);
202  if (th)
203  CHKERR moab.tag_get_data(th, conn_skin, num_nodes, &coords_skin[0]);
204  else
205  CHKERR moab.get_coords(conn_skin, num_nodes, &coords_skin[0]);
207  &coords_skin[0], &coords_skin[1], &coords_skin[2]);
209  &coords_skin[3], &coords_skin[4], &coords_skin[5]);
210  FTensor::Tensor1<double, 3> t_skin_delta;
211  t_skin_delta(i) = t_n1(i) - t_n0(i);
212  const double skin_edge_length = sqrt(t_skin_delta(i) * t_skin_delta(i));
213  for (int nn = 0; nn != 2; ++nn) {
214  auto it = map_verts_length.find(conn_skin[nn]);
215  if (it != map_verts_length.end())
216  it->second = std::min(it->second, skin_edge_length);
217  else
218  map_verts_length[conn_skin[nn]] = skin_edge_length;
219  }
220  }
221 
222  for (auto m : map_verts_length) {
223 
225  if (th)
226  CHKERR moab.tag_get_data(th, &m.first, 1, &t_n(0));
227  else
228  CHKERR moab.get_coords(&m.first, 1, &t_n(0));
229 
230  double min_dist = rel_tol * m.second;
231  FTensor::Tensor1<double, 3> t_min_coords;
232  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
233  &t_n(0), fixed_edges, &min_dist, &t_min_coords(0));
234 
235  if (min_dist < rel_tol * m.second || min_dist < abs_tol) {
236  if (th)
237  CHKERR moab.tag_set_data(th, &m.first, 1, &t_min_coords(0));
238  else
239  CHKERR moab.set_coords(&m.first, 1, &t_min_coords(0));
240  }
241  }
242 
244 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:86
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:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

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

1941  {
1942  CoreInterface &m_field = cOre;
1943  moab::Interface &moab = m_field.get_moab();
1944  PrismInterface *interface;
1946  CHKERR m_field.getInterface(interface);
1947  EntityHandle meshset_volume;
1948  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
1949  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1950  split_bit, BitRefLevel().set(), MBTET, meshset_volume);
1951  EntityHandle meshset_surface;
1952  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
1953  CHKERR moab.add_entities(meshset_surface, ents);
1954  CHKERR interface->getSides(meshset_surface, split_bit, true);
1955  CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
1956  true);
1957  CHKERR moab.delete_entities(&meshset_volume, 1);
1958  CHKERR moab.delete_entities(&meshset_surface, 1);
1959  if (th) {
1960  Range prisms;
1961  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1962  bit, BitRefLevel().set(), MBPRISM, prisms);
1963  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
1964  int num_nodes;
1965  const EntityHandle *conn;
1966  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
1967  MatrixDouble data(3, 3);
1968  CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
1969  CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
1970  }
1971  }
1973 }
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ trimEdgesInTheMiddle()

MoFEMErrorCode MoFEM::CutMeshInterface::trimEdgesInTheMiddle ( const BitRefLevel  bit,
Tag  th = NULL,
const double  tol = 1e-4,
const bool  debug = false 
)

trim edges

Parameters
bitbit level of the trimmed mesh
Returns
error code

Definition at line 1812 of file CutMeshInterface.cpp.

1814  {
1815  CoreInterface &m_field = cOre;
1816  moab::Interface &moab = m_field.get_moab();
1817  MeshRefinement *refiner;
1818  const RefEntity_multiIndex *ref_ents_ptr;
1820 
1821  CHKERR m_field.getInterface(refiner);
1822  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1823  CHKERR refiner->add_verices_in_the_middel_of_edges(trimEdges, bit);
1824  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
1825  debug ? &trimEdges : NULL);
1826 
1827  trimNewVolumes.clear();
1828  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1829  bit, BitRefLevel().set(), MBTET, trimNewVolumes);
1830 
1831  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
1832  mit != edgesToTrim.end(); mit++) {
1833  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1834  boost::make_tuple(MBVERTEX, mit->first));
1835  if (vit ==
1836  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1837  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1838  "No vertex on trim edges, that make no sense");
1839  }
1840  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1841  if ((ref_ent->getBitRefLevel() & bit).any()) {
1842  EntityHandle vert = ref_ent->getRefEnt();
1843  trimNewVertices.insert(vert);
1844  verticesOnTrimEdges[vert] = mit->second;
1845  }
1846  }
1847 
1848  // Get faces which are trimmed
1849  trimNewSurfaces.clear();
1850  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1851  bit, BitRefLevel().set(), MBTRI, trimNewSurfaces);
1852 
1853  Range trim_new_surfaces_nodes;
1854  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
1855  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
1856  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
1857  Range faces_not_on_surface;
1858  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
1859  faces_not_on_surface, moab::Interface::UNION);
1860  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
1861 
1862  // Get surfaces which are not trimmed and add them to surface
1863  Range all_surfaces_on_bit_level;
1864  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1865  bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
1866  all_surfaces_on_bit_level =
1867  intersect(all_surfaces_on_bit_level, cutNewSurfaces);
1868  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
1869 
1870  Range trim_surface_edges;
1871  CHKERR moab.get_adjacencies(trimNewSurfaces, 1, false, trim_surface_edges,
1872  moab::Interface::UNION);
1873 
1874  // check of nodes are outside surface and if it are remove adjacent faces to
1875  // those nodes.
1876  Range check_verts;
1877  CHKERR moab.get_connectivity(trimNewSurfaces, check_verts, true);
1878  check_verts = subtract(check_verts, trimNewVertices);
1879  for (auto v : check_verts) {
1880 
1881  VectorDouble3 s(3);
1882  if (th) {
1883  CHKERR moab.tag_get_data(th, &v, 1, &s[0]);
1884  } else {
1885  CHKERR moab.get_coords(&v, 1, &s[0]);
1886  }
1887 
1888  VectorDouble3 p(3);
1889  EntityHandle facets_out;
1890  CHKERR treeSurfPtr->closest_to_location(&s[0], rootSetSurf, &p[0],
1891  facets_out);
1892  VectorDouble3 n(3);
1893  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
1894  n /= norm_2(n);
1895  VectorDouble3 delta = s - p;
1896  VectorDouble3 normal = inner_prod(delta, n) * n;
1897  if (norm_2(delta - normal) > tol * aveLength) {
1898  Range adj;
1899  CHKERR moab.get_adjacencies(&v, 1, 2, false, adj);
1900  trimNewSurfaces = subtract(trimNewSurfaces, adj);
1901  }
1902  }
1903 
1905 }
map< EntityHandle, TreeData > edgesToTrim
map< EntityHandle, TreeData > verticesOnTrimEdges
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
double aveLength
Average edge length.
double tol
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:594
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

Member Data Documentation

◆ aveLength

double MoFEM::CutMeshInterface::aveLength
private

Average edge length.

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

◆ cutNewSurfaces

Range MoFEM::CutMeshInterface::cutNewSurfaces
private

Definition at line 364 of file CutMeshInterface.hpp.

◆ cutNewVertices

Range MoFEM::CutMeshInterface::cutNewVertices
private

Definition at line 367 of file CutMeshInterface.hpp.

◆ cutNewVolumes

Range MoFEM::CutMeshInterface::cutNewVolumes
private

Definition at line 363 of file CutMeshInterface.hpp.

◆ cutVolumes

Range MoFEM::CutMeshInterface::cutVolumes
private

Definition at line 362 of file CutMeshInterface.hpp.

◆ edgesToCut

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

Definition at line 386 of file CutMeshInterface.hpp.

◆ edgesToTrim

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

Definition at line 388 of file CutMeshInterface.hpp.

◆ fRont

Range MoFEM::CutMeshInterface::fRont
private

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

◆ mergedSurfaces

Range MoFEM::CutMeshInterface::mergedSurfaces
private

Definition at line 375 of file CutMeshInterface.hpp.

◆ mergedVolumes

Range MoFEM::CutMeshInterface::mergedVolumes
private

Definition at line 374 of file CutMeshInterface.hpp.

◆ moabTetGenMap

map<EntityHandle, unsigned long> MoFEM::CutMeshInterface::moabTetGenMap
private

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

◆ sUrface

Range MoFEM::CutMeshInterface::sUrface
private

Definition at line 355 of file CutMeshInterface.hpp.

◆ tetGenData

boost::ptr_vector<tetgenio> MoFEM::CutMeshInterface::tetGenData
private

Definition at line 395 of file CutMeshInterface.hpp.

◆ tetGenMoabMap

map<unsigned long, EntityHandle> MoFEM::CutMeshInterface::tetGenMoabMap
private

Definition at line 394 of file CutMeshInterface.hpp.

◆ tetgenSurfaces

Range MoFEM::CutMeshInterface::tetgenSurfaces
private

Definition at line 377 of file CutMeshInterface.hpp.

◆ treeSurfPtr

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

Definition at line 358 of file CutMeshInterface.hpp.

◆ trimEdges

Range MoFEM::CutMeshInterface::trimEdges
private

Definition at line 372 of file CutMeshInterface.hpp.

◆ trimNewSurfaces

Range MoFEM::CutMeshInterface::trimNewSurfaces
private

Definition at line 371 of file CutMeshInterface.hpp.

◆ trimNewVertices

Range MoFEM::CutMeshInterface::trimNewVertices
private

Definition at line 370 of file CutMeshInterface.hpp.

◆ trimNewVolumes

Range MoFEM::CutMeshInterface::trimNewVolumes
private

Definition at line 369 of file CutMeshInterface.hpp.

◆ verticesOnCutEdges

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

Definition at line 387 of file CutMeshInterface.hpp.

◆ verticesOnTrimEdges

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

Definition at line 389 of file CutMeshInterface.hpp.

◆ vOlume

Range MoFEM::CutMeshInterface::vOlume
private

Definition at line 356 of file CutMeshInterface.hpp.

◆ zeroDistanceEnts

Range MoFEM::CutMeshInterface::zeroDistanceEnts
private

Definition at line 365 of file CutMeshInterface.hpp.

◆ zeroDistanceVerts

Range MoFEM::CutMeshInterface::zeroDistanceVerts
private

Definition at line 366 of file CutMeshInterface.hpp.


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