v0.8.13
Classes | Public Member Functions | Public Attributes | Private Member Functions | 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 getOptions ()
 Get options from command line. More...
 
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 buildTree ()
 build tree More...
 
MoFEMErrorCode cutAndTrim (const BitRefLevel &bit_level1, const BitRefLevel &bit_level2, 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=true)
 
MoFEMErrorCode cutTrimAndMerge (const int fraction_level, const BitRefLevel &bit_level1, const BitRefLevel &bit_level2, const BitRefLevel &bit_level3, 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=0, int verb=0, const bool debug=false)
 find edges to cut More...
 
MoFEMErrorCode projectZeroDistanceEnts (Range *fixed_edges, Range *corner_nodes, const double low_tol=0, const int verb=QUIET, const bool debug=false)
 
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, int verb=0)
 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 & 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 ()
 
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
 
int nbMaxTrimSearchIterations
 

Private Member Functions

MoFEMErrorCode getRayForEdge (const EntityHandle ent, VectorAdaptor &ray_point, VectorAdaptor &unit_ray_dir, double &ray_length) const
 

Private Attributes

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

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

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

135  {
136  CoreInterface &m_field = cOre;
137  moab::Interface &moab = m_field.get_moab();
139  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
140  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
143 }
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
#define CHKERR
Inline error check.
Definition: definitions.h:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

◆ clearMap()

MoFEMErrorCode MoFEM::CutMeshInterface::clearMap ( )

Definition at line 40 of file CutMeshInterface.cpp.

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

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

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

◆ cutAndTrim()

MoFEMErrorCode MoFEM::CutMeshInterface::cutAndTrim ( const BitRefLevel bit_level1,
const BitRefLevel bit_level2,
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 = true 
)

Definition at line 156 of file CutMeshInterface.cpp.

160  {
161  CoreInterface &m_field = cOre;
163 
164  // cut mesh
165  CHKERR findEdgesToCut(fixed_edges, corner_nodes, tol_cut, QUIET, debug);
166  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_cut_close,
167  QUIET, debug);
168  CHKERR cutEdgesInMiddle(bit_level1, debug);
169  if (fixed_edges) {
170  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
171  *fixed_edges);
172  }
173  if (corner_nodes) {
174  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
175  *corner_nodes);
176  }
177  if (update_meshsets) {
178  CHKERR UpdateMeshsets()(cOre, bit_level1);
179  }
181 
182  auto get_min_quality = [&m_field](const BitRefLevel bit, Tag th) {
183  Range tets_level; // test at level
184  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
185  bit, BitRefLevel().set(), MBTET, tets_level);
186  double min_q = 1;
187  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
188  return min_q;
189  };
190 
191  PetscPrintf(PETSC_COMM_WORLD, "Min quality cut %6.4g\n",
192  get_min_quality(bit_level1, th));
193 
194  if(debug) {
196  }
197 
198  // trim mesh
199  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim);
200  CHKERR trimEdgesInTheMiddle(bit_level2, th, tol_trim_close, debug);
201  if (fixed_edges) {
202  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
203  *fixed_edges);
204  }
205  if (corner_nodes) {
206  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
207  *corner_nodes);
208  }
209  if (update_meshsets) {
210  CHKERR UpdateMeshsets()(cOre, bit_level2);
211  }
213 
214  PetscPrintf(PETSC_COMM_WORLD, "Min quality trim %3.2g\n",
215  get_min_quality(bit_level2, th));
216 
217  if(debug) {
219  }
220 
222 }
MoFEMErrorCode findEdgesToCut(Range *fixed_edges, Range *corner_nodes, const double low_tol=0, int verb=0, const bool debug=false)
find edges to cut
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, int verb=0)
Find edges to trimEdges.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
const Range & projectZeroDistanceEnts() const
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
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:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439
MoFEMErrorCode saveTrimEdges()
MoFEMErrorCode trimEdgesInTheMiddle(const BitRefLevel bit, Tag th=NULL, const double tol=1e-4, const bool debug=false)
trim edges
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 944 of file CutMeshInterface.cpp.

945  {
946  CoreInterface &m_field = cOre;
947  moab::Interface &moab = m_field.get_moab();
948  MeshRefinement *refiner;
949  const RefEntity_multiIndex *ref_ents_ptr;
951  if(cutEdges.size() != edgesToCut.size()) {
952  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
953  }
954  CHKERR m_field.getInterface(refiner);
955  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
956  CHKERR refiner->add_verices_in_the_middel_of_edges(cutEdges, bit);
957  CHKERR refiner->refine_TET(vOlume, bit, false, QUIET,
958  debug ? &cutEdges : NULL);
959  cutNewVolumes.clear();
960  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
961  bit, BitRefLevel().set(), MBTET, cutNewVolumes);
962  cutNewSurfaces.clear();
963  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
964  bit, BitRefLevel().set(), MBTRI, cutNewSurfaces);
965  // Find new vertices on cut edges
966  cutNewVertices.clear();
967  CHKERR moab.get_connectivity(zeroDistanceEnts, cutNewVertices, true);
969  for (map<EntityHandle, TreeData>::iterator mit = edgesToCut.begin();
970  mit != edgesToCut.end(); ++mit) {
971  RefEntity_multiIndex::index<
972  Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vit =
973  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
974  boost::make_tuple(MBVERTEX, mit->first));
975  if (vit ==
976  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
977  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
978  "No vertex on cut edges, that make no sense");
979  }
980  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
981  if ((ref_ent->getBitRefLevel() & bit).any()) {
982  EntityHandle vert = ref_ent->getRefEnt();
983  cutNewVertices.insert(vert);
984  verticesOnCutEdges[vert] = mit->second;
985  } else {
986  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
987  "Vertex has wrong bit ref level");
988  }
989  }
990  // Add zero distance entities faces
991  Range tets_skin;
992  Skinner skin(&moab);
993  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
994  cutNewSurfaces.merge(zeroDistanceEnts.subset_by_type(MBTRI));
995  // At that point cutNewSurfaces has all newly created faces, now take all
996  // nodes on those faces and subtract nodes on catted edges. Faces adjacent to
997  // nodes which left are not part of surface.
998  Range diff_verts;
999  CHKERR moab.get_connectivity(unite(cutNewSurfaces, zeroDistanceEnts),
1000  diff_verts, true);
1001  diff_verts = subtract(diff_verts, cutNewVertices);
1002  Range subtract_faces;
1003  CHKERR moab.get_adjacencies(diff_verts, 2, false, subtract_faces,
1004  moab::Interface::UNION);
1005  cutNewSurfaces = subtract(cutNewSurfaces, unite(subtract_faces, tets_skin));
1006  cutNewVertices.clear();
1007  CHKERR moab.get_connectivity(cutNewSurfaces, cutNewVertices, true);
1008 
1010 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:614
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:439
map< EntityHandle, TreeData > edgesToCut

◆ cutTrimAndMerge()

MoFEMErrorCode MoFEM::CutMeshInterface::cutTrimAndMerge ( const int  fraction_level,
const BitRefLevel bit_level1,
const BitRefLevel bit_level2,
const BitRefLevel bit_level3,
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 224 of file CutMeshInterface.cpp.

229  {
230  CoreInterface &m_field = cOre;
232  if(debug) {
233  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
234  "ents_not_in_database.vtk", "VTK", "");
235  }
236  CHKERR cutAndTrim(bit_level1, bit_level2, th, tol_cut, tol_cut_close,
237  tol_trim, tol_trim_close, &fixed_edges, &corner_nodes,
238  update_meshsets, debug);
239  if(debug) {
240  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
241  "cut_trim_ents_not_in_database.vtk", "VTK", "");
242  }
243 
244  CHKERR mergeBadEdges(fraction_level, bit_level2, bit_level1, bit_level3,
245  getNewTrimSurfaces(), fixed_edges, corner_nodes, th,
246  update_meshsets, debug);
247  // CHKERR removePathologicalFrontTris(bit_level3,
248  // const_cast<Range &>(getMergedSurfaces()));
249 
250  if(debug) {
251  CHKERR cOre.getInterface<BitRefManager>()->writeBitLevelByType(
252  bit_level3, BitRefLevel().set(), MBTET, "out_tets_merged.vtk", "VTK",
253  "");
254  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
255  "cut_trim_merge_ents_not_in_database.vtk", "VTK", "");
256  }
257 
258  auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
259  Range tets_level; // test at level
260  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
261  bit, BitRefLevel().set(), MBTET, tets_level);
262  double min_q = 1;
263  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
264  if (min_q < 0 && debug) {
265  CHKERR m_field.getInterface<Tools>()->writeTetsWithQuality(
266  "negative_tets.vtk", "VTK", "", tets_level, th);
267  }
268  return min_q;
269  };
270 
271  PetscPrintf(PETSC_COMM_WORLD, "Min quality node merge %6.4g\n",
272  get_min_quality(bit_level1, th));
273 
274  CHKERR
275  cOre.getInterface<BitRefManager>()->updateRange(fixed_edges, fixed_edges);
276  CHKERR cOre.getInterface<BitRefManager>()->updateRange(corner_nodes,
277  corner_nodes);
278 
280 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
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 cutAndTrim(const BitRefLevel &bit_level1, const BitRefLevel &bit_level2, 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=true)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
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:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

◆ findEdgesToCut()

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

find edges to cut

Parameters
verbverbosity level
Returns
error code

Definition at line 313 of file CutMeshInterface.cpp.

316  {
317  CoreInterface &m_field = cOre;
318  moab::Interface &moab = m_field.get_moab();
320 
321  edgesToCut.clear();
322  cutEdges.clear();
323 
324  zeroDistanceVerts.clear();
325  zeroDistanceEnts.clear();
326  verticesOnCutEdges.clear();
327 
328  double ray_length;
329  double ray_point[3], unit_ray_dir[3];
330  VectorAdaptor vec_unit_ray_dir(
331  3, ublas::shallow_array_adaptor<double>(3, unit_ray_dir));
332  VectorAdaptor vec_ray_point(
333  3, ublas::shallow_array_adaptor<double>(3, ray_point));
334 
335  Tag th_dist;
336  rval = moab.tag_get_handle("DIST",th_dist);
337  if(rval == MB_SUCCESS) {
338  CHKERR moab.tag_delete(th_dist);
339  } else {
340  rval = MB_SUCCESS;
341  }
342  Tag th_dist_normal;
343  rval = moab.tag_get_handle("DIST_NORMAL", th_dist_normal);
344  if(rval == MB_SUCCESS) {
345  CHKERR moab.tag_delete(th_dist_normal);
346  } else {
347  rval = MB_SUCCESS;
348  }
349 
350  double def_val[] = {0, 0, 0};
351  CHKERR moab.tag_get_handle("DIST", 1, MB_TYPE_DOUBLE, th_dist,
352  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
353  CHKERR moab.tag_get_handle("DIST_NORMAL", 3, MB_TYPE_DOUBLE, th_dist_normal,
354  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
355 
356  Range vol_vertices;
357  CHKERR moab.get_connectivity(vOlume, vol_vertices, true);
358  for (auto v : vol_vertices) {
359  VectorDouble3 point_in(3);
360  CHKERR moab.get_coords(&v, 1, &point_in[0]);
361  VectorDouble3 point_out(3);
362  EntityHandle facets_out;
363  CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
364  &point_out[0], facets_out);
365  VectorDouble3 n(3);
366  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
367  VectorDouble3 delta = point_out - point_in;
368  double dist = norm_2(delta);
369  VectorDouble3 dist_normal = inner_prod(delta, n) * n;
370  CHKERR moab.tag_set_data(th_dist, &v, 1, &dist);
371  CHKERR moab.tag_set_data(th_dist_normal, &v, 1, &dist_normal[0]);
372  }
373 
374  auto get_normal_dist = [](const double *normal) {
375  FTensor::Tensor1<double, 3> t_n(normal[0], normal[1], normal[2]);
377  return sqrt(t_n(i) * t_n(i));
378  };
379 
380  auto get_edge_crossed = [&moab, get_normal_dist,
381  th_dist_normal](EntityHandle v0, EntityHandle v1) {
382  VectorDouble3 dist_normal0(3);
383  CHKERR moab.tag_get_data(th_dist_normal, &v0, 1, &dist_normal0[0]);
384  VectorDouble3 dist_normal1(3);
385  CHKERR moab.tag_get_data(th_dist_normal, &v1, 1, &dist_normal1[0]);
386  return (inner_prod(dist_normal0, dist_normal1) < 0);
387  };
388 
389  auto get_normal_dist_from_conn = [&moab, get_normal_dist,
390  th_dist_normal](EntityHandle v) {
391  double dist_normal[3];
392  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, dist_normal);
393  return get_normal_dist(dist_normal);
394  };
395 
396  auto project_node = [this, &moab, th_dist_normal](const EntityHandle v) {
398  VectorDouble3 dist_normal(3);
399  rval = moab.tag_get_data(th_dist_normal, &v, 1, &dist_normal[0]);
400  VectorDouble3 s0(3);
401  CHKERR moab.get_coords(&v, 1, &s0[0]);
402  double dist = norm_2(dist_normal);
403  verticesOnCutEdges[v].dIst = dist;
404  verticesOnCutEdges[v].lEngth = dist;
405  verticesOnCutEdges[v].unitRayDir =
406  dist > 0 ? dist_normal / dist : dist_normal;
407  verticesOnCutEdges[v].rayPoint = s0;
409  };
410 
411  auto not_project_node = [this, &moab](const EntityHandle v) {
413  VectorDouble3 s0(3);
414  CHKERR moab.get_coords(&v, 1, &s0[0]);
415  verticesOnCutEdges[v].dIst = 0;
416  verticesOnCutEdges[v].lEngth = 0;
417  verticesOnCutEdges[v].unitRayDir = s0;
418  verticesOnCutEdges[v].rayPoint = s0;
420  };
421 
422  auto check_if_is_on_fixed_edge = [fixed_edges](const EntityHandle e) {
423  if(fixed_edges) {
424  if(fixed_edges->find(e)!=fixed_edges->end()) {
425  return true;
426  } else {
427  return false;
428  }
429  } else {
430  return false;
431  }
432  };
433 
434  auto check_if_is_on_cornet_node = [corner_nodes](const EntityHandle v) {
435  if (corner_nodes) {
436  if (corner_nodes->find(v) != corner_nodes->end()) {
437  return true;
438  } else {
439  return false;
440  }
441  } else {
442  return false;
443  }
444  };
445 
446  Range vol_edges;
447  CHKERR moab.get_adjacencies(vOlume, 1, true, vol_edges,
448  moab::Interface::UNION);
449 
450  aveLength = 0;
451  maxLength = 0;
452  int nb_ave_length = 0;
453  for (auto e : vol_edges) {
454  int num_nodes;
455  const EntityHandle *conn;
456  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
457  double dist[num_nodes];
458  CHKERR moab.tag_get_data(th_dist, conn, num_nodes, dist);
459  CHKERR getRayForEdge(e, vec_ray_point, vec_unit_ray_dir, ray_length);
460  const double tol = ray_length * low_tol;
461  if (get_edge_crossed(conn[0], conn[1])) {
462  std::vector<double> distances_out;
463  std::vector<EntityHandle> facets_out;
464  CHKERR treeSurfPtr->ray_intersect_triangles(distances_out, facets_out,
465  rootSetSurf, tol, ray_point,
466  unit_ray_dir, &ray_length);
467  if (!distances_out.empty()) {
468  const auto dist_ptr =
469  std::min_element(distances_out.begin(), distances_out.end());
470  const double dist = *dist_ptr;
471  if (dist <= ray_length) {
472  aveLength += ray_length;
473  maxLength = fmax(maxLength, ray_length);
474  nb_ave_length++;
475  edgesToCut[e].dIst = dist;
476  edgesToCut[e].lEngth = ray_length;
477  edgesToCut[e].unitRayDir = vec_unit_ray_dir;
478  edgesToCut[e].rayPoint = vec_ray_point;
479  cutEdges.insert(e);
480  }
481  }
482  }
483 
484  if (fabs(dist[0]) < tol && fabs(dist[1]) < tol) {
485  aveLength += ray_length;
486  maxLength = fmax(maxLength, ray_length);
487  if (check_if_is_on_fixed_edge(e)) {
488  CHKERR not_project_node(conn[0]);
489  CHKERR not_project_node(conn[1]);
490  } else {
491  CHKERR project_node(conn[0]);
492  CHKERR project_node(conn[1]);
493  }
494  zeroDistanceEnts.insert(e);
495  }
496  }
497  aveLength /= nb_ave_length;
498 
499  Range cut_edges_verts;
500  CHKERR moab.get_connectivity(unite(cutEdges, zeroDistanceEnts),
501  cut_edges_verts, true);
502  vol_vertices = subtract(vol_vertices, cut_edges_verts);
503 
504  for (auto v : vol_vertices) {
505  double dist;
506  CHKERR moab.tag_get_data(th_dist, &v, 1, &dist);
507  const double tol = get_ave_edge_length(v, vol_edges, moab) * low_tol;
508  if(fabs(dist) < tol) {
509 
510  if (check_if_is_on_cornet_node(v)) {
511  CHKERR not_project_node(v);
512  } else {
513  CHKERR project_node(v);
514  }
515 
516  zeroDistanceVerts.insert(v);
517  }
518  }
519 
520  cutVolumes.clear();
521  // take all volumes adjacent to cut edges
522  CHKERR moab.get_adjacencies(cutEdges, 3, false, cutVolumes,
523  moab::Interface::UNION);
524  CHKERR moab.get_adjacencies(zeroDistanceVerts, 3, false, cutVolumes,
525  moab::Interface::UNION);
526  {
527  Range verts;
528  CHKERR moab.get_connectivity(unite(cutEdges, zeroDistanceEnts), verts,
529  true);
530  CHKERR moab.get_adjacencies(verts, 3, false, cutVolumes,
531  moab::Interface::UNION);
532  }
533  cutVolumes = intersect(cutVolumes,vOlume);
534 
535  // get edges on the cut volumes
536  Range edges;
537  CHKERR moab.get_adjacencies(cutVolumes, 1, false, edges,
538  moab::Interface::UNION);
539  edges = subtract(edges, cutEdges);
540 
541  // add to cut set edges which are cutted by extension of cutting surface
542  for (auto e : edges) {
543  int num_nodes;
544  const EntityHandle *conn;
545  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
546  const double tol = get_ave_edge_length(e, vol_edges, moab) * low_tol;
547  double dist_normal[2];
548  dist_normal[0] = get_normal_dist_from_conn(conn[0]);
549  dist_normal[1] = get_normal_dist_from_conn(conn[1]);
550  if (get_edge_crossed(conn[0], conn[1])) {
551  CHKERR getRayForEdge(e, vec_ray_point, vec_unit_ray_dir, ray_length);
552  double s =
553  fabs(dist_normal[0]) / (fabs(dist_normal[0]) + fabs(dist_normal[1]));
554  edgesToCut[e].dIst = s * ray_length;
555  edgesToCut[e].lEngth = ray_length;
556  edgesToCut[e].unitRayDir = vec_unit_ray_dir;
557  edgesToCut[e].rayPoint = vec_ray_point;
558  cutEdges.insert(e);
559  } else if (fabs(dist_normal[0]) < tol && fabs(dist_normal[1]) < tol) {
560  if (check_if_is_on_fixed_edge(e)) {
561  CHKERR not_project_node(conn[0]);
562  CHKERR not_project_node(conn[1]);
563  } else {
564  CHKERR project_node(conn[0]);
565  CHKERR project_node(conn[1]);
566  }
567  zeroDistanceEnts.insert(e);
568  }
569  }
570 
571  CHKERR moab.get_adjacencies(cutVolumes, 1, false, edges,
572  moab::Interface::UNION);
573  Range add_verts;
574  CHKERR moab.get_connectivity(edges, add_verts, true);
575  add_verts = subtract(add_verts, zeroDistanceVerts);
576  CHKERR moab.get_connectivity(unite(cutEdges, zeroDistanceEnts),
577  cut_edges_verts, true);
578  add_verts = subtract(add_verts, cut_edges_verts);
579 
580  for (auto v : add_verts) {
581  double dist_normal = get_normal_dist_from_conn(v);
582  const double tol = get_ave_edge_length(v, vol_edges, moab) * low_tol;
583  if (fabs(dist_normal) < tol) {
584 
585  if (check_if_is_on_cornet_node(v)) {
586  CHKERR not_project_node(v);
587  } else {
588  CHKERR project_node(v);
589  }
590 
591  zeroDistanceVerts.insert(v);
592  }
593  }
594 
595  for (auto f : zeroDistanceEnts) {
596  int num_nodes;
597  const EntityHandle *conn;
598  CHKERR moab.get_connectivity(f, conn, num_nodes, true);
599  Range adj_edges;
600  CHKERR moab.get_adjacencies(conn, num_nodes, 1, false, adj_edges,
601  moab::Interface::UNION);
602  for (auto e : adj_edges) {
603  cutEdges.erase(e);
604  edgesToCut.erase(e);
605  }
606  }
607 
608  for (auto v : zeroDistanceVerts) {
609  Range adj_edges;
610  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges,
611  moab::Interface::UNION);
612  for (auto e : adj_edges) {
613  cutEdges.erase(e);
614  edgesToCut.erase(e);
615  }
616  }
617 
618  if(debug) {
619  EntityHandle out_meshset;
620  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, out_meshset);
621  CHKERR m_field.get_moab().add_entities(out_meshset, cutEdges);
622  CHKERR m_field.get_moab().write_file("cut_edges.vtk", "VTK", "",
623  &out_meshset, 1);
624  CHKERR m_field.get_moab().delete_entities(&out_meshset, 1);
625  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, out_meshset);
626  CHKERR m_field.get_moab().add_entities(out_meshset, zeroDistanceEnts);
627  CHKERR m_field.get_moab().add_entities(out_meshset, zeroDistanceVerts);
628  CHKERR m_field.get_moab().write_file("zero_dist_ents.vtk", "VTK", "",
629  &out_meshset, 1);
630  CHKERR m_field.get_moab().delete_entities(&out_meshset, 1);
631  }
632 
634 }
double maxLength
Maximal edge length.
MoFEMErrorCode getRayForEdge(const EntityHandle ent, VectorAdaptor &ray_point, VectorAdaptor &unit_ray_dir, double &ray_length) const
static double get_ave_edge_length(const EntityHandle ent, const Range &vol_edges, moab::Interface &moab)
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
double aveLength
Average edge length.
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Common.hpp:234
double tol
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:614
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:439
map< EntityHandle, TreeData > edgesToCut
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Common.hpp:220

◆ findEdgesToTrim()

MoFEMErrorCode MoFEM::CutMeshInterface::findEdgesToTrim ( Range *  fixed_edges,
Range *  corner_nodes,
Tag  th = NULL,
const double  tol = 1e-4,
int  verb = 0 
)

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

1058  {
1059  CoreInterface &m_field = cOre;
1060  moab::Interface &moab = m_field.get_moab();
1062 
1063  // takes edges on body skin
1064  Skinner skin(&moab);
1065  Range tets_skin;
1066  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
1067  // vertives on the skin
1068  Range tets_skin_verts;
1069  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1070  // edges on the skin
1071  Range tets_skin_edges;
1072  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1073  moab::Interface::UNION);
1074  // get edges on new surface
1075  Range edges;
1076  CHKERR moab.get_adjacencies(cutNewSurfaces, 1, false, edges,
1077  moab::Interface::UNION);
1078 
1079  Range cut_surface_edges_on_fixed_edges;
1080  if (fixed_edges) {
1081  cut_surface_edges_on_fixed_edges = intersect(edges, *fixed_edges);
1082  }
1083  Range cut_surface_edges_on_fixed_edges_verts;
1084  CHKERR moab.get_connectivity(cut_surface_edges_on_fixed_edges,
1085  cut_surface_edges_on_fixed_edges_verts, true);
1086  Range fixed_edges_verts_on_corners;
1087  if (fixed_edges) {
1088  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts_on_corners,
1089  true);
1090  }
1091  fixed_edges_verts_on_corners = subtract(
1092  fixed_edges_verts_on_corners, cut_surface_edges_on_fixed_edges_verts);
1093  if (corner_nodes) {
1094  fixed_edges_verts_on_corners.merge(*corner_nodes);
1095  }
1096 
1097  // clear data ranges
1098  trimEdges.clear();
1099  edgesToTrim.clear();
1100  verticesOnTrimEdges.clear();
1101  trimNewVertices.clear();
1102 
1103  // iterate over entities on new cut surface
1104  std::multimap<double, std::pair<EntityHandle, EntityHandle>> verts_map;
1105  for (auto e : edges) {
1106  // Get edge connectivity and coords
1107  int num_nodes;
1108  const EntityHandle *conn;
1109  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1110  double coords[3 * num_nodes];
1111  if (th) {
1112  CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
1113  } else {
1114  CHKERR moab.get_coords(conn, num_nodes, coords);
1115  }
1116  // Put edges coords into boost vectors
1117  auto get_s_adaptor = [&coords](const int n) {
1118  return VectorAdaptor(3,
1119  ublas::shallow_array_adaptor<double>(3, &coords[n]));
1120  };
1121  VectorAdaptor s0 = get_s_adaptor(0);
1122  VectorAdaptor s1 = get_s_adaptor(3);
1123  // get edge length
1124  double length = norm_2(s1 - s0);
1125 
1126  // Find point on surface closet to surface
1127  auto get_closets_delta = [this, &moab](const VectorAdaptor &s) {
1128  VectorDouble3 p(3);
1129  EntityHandle facets_out;
1130  // find closet point on the surface from first node
1131  CHKERR treeSurfPtr->closest_to_location(&s[0], rootSetSurf, &p[0],
1132  facets_out);
1133  VectorDouble3 n(3);
1134  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
1135  VectorDouble3 w = p - s;
1136  VectorDouble3 normal = inner_prod(w, n) * n;
1137  w -= normal;
1138  return w;
1139  };
1140 
1141  // Calculate deltas, i.e. vectors from edges to closet point on surface
1142  VectorDouble3 delta0(3), delta1(3);
1143  noalias(delta0) = get_closets_delta(s0);
1144  noalias(delta1) = get_closets_delta(s1);
1145 
1146  // moab.tag_set_data(th,&conn[0],1,&delta0[0]);
1147  // moab.tag_set_data(th,&conn[1],1,&delta1[0]);
1148  // Calculate distances
1149  double dist0 = norm_2(delta0);
1150  double dist1 = norm_2(delta1);
1151  double min_dist = fmin(dist0, dist1);
1152  double max_dist = fmax(dist0, dist1);
1153 
1154  // add edge to trim
1155  double dist;
1156  VectorDouble3 ray;
1157  VectorDouble3 trimmed_end;
1158  VectorDouble3 itersection_point;
1159 
1160  if (min_dist < 1e-6 * aveLength && max_dist >= 1e-6 * aveLength) {
1161  if (max_dist == dist0) {
1162  // move mid node in reference to node 0
1163  trimmed_end = s0;
1164  ray = s1 - trimmed_end;
1165  } else {
1166  // move node in reference to node 1
1167  trimmed_end = s1;
1168  ray = s0 - trimmed_end;
1169  }
1170 
1171  // Solve nonlinera problem of finding point on surface front
1172  auto closest_point_projection =
1173  [this, &moab](VectorDouble3 ray, VectorDouble3 trimmed_end,
1174  const int max_it, const double tol) {
1175  VectorDouble3 n(3), w(3), normal(3);
1176  double length = norm_2(ray);
1177  ray /= length;
1178  for (int ii = 0; ii != max_it; ii++) {
1179  EntityHandle facets_out;
1180  VectorDouble3 point_out(3);
1181  treeSurfPtr->closest_to_location(&trimmed_end[0], rootSetSurf,
1182  &point_out[0], facets_out);
1183  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
1184  noalias(w) = point_out - trimmed_end;
1185  noalias(normal) = inner_prod(w, n) * n;
1186  double s = inner_prod(ray, w - normal);
1187  trimmed_end += s * ray;
1188  // cerr << "s " << ii << " " << s << " " << norm_2(w) << endl;
1189  if ((s / length) < tol)
1190  break;
1191  }
1192  return trimmed_end;
1193  };
1194 
1195  itersection_point = closest_point_projection(
1196  ray, trimmed_end, nbMaxTrimSearchIterations, 1e-12);
1197 
1198  ray = itersection_point - trimmed_end;
1199  dist = norm_2(ray);
1200 
1201  if ((1 - dist / length) > 0) {
1202 
1203  // check if edges should be trimmed, i.e. if edge is trimmed at very
1204  // end simply move closed node rather than trim
1205  edgesToTrim[e].dIst = dist;
1206  edgesToTrim[e].lEngth = dist;
1207  edgesToTrim[e].unitRayDir = ray / dist;
1208  edgesToTrim[e].rayPoint = trimmed_end;
1209  trimEdges.insert(e);
1210 
1211  auto add_vert = [&verts_map, e](EntityHandle v, double dist) {
1212  verts_map.insert(
1213  std::pair<double, std::pair<EntityHandle, EntityHandle>>(
1214  dist, std::pair<EntityHandle, EntityHandle>(v, e)));
1215  };
1216 
1217  double dist0_to_intersection =
1218  norm_2(itersection_point - s0) / aveLength;
1219  double dist1_to_intersection =
1220  norm_2(itersection_point - s1) / aveLength;
1221  if (dist0_to_intersection < dist1_to_intersection) {
1222  add_vert(conn[0], dist0_to_intersection);
1223  } else {
1224  add_vert(conn[1], dist1_to_intersection);
1225  }
1226  }
1227  }
1228  }
1229 
1230  // EntityHandle meshset;
1231  // CHKERR moab.create_meshset(MESHSET_SET, meshset);
1232  // Tag th_aaa;
1233  // double def_val[] = {0,0,0};
1234  // CHKERR moab.tag_get_handle("AAAA", 3, MB_TYPE_DOUBLE, th_aaa,
1235  // MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
1236  // Range all_nodes;
1237  // CHKERR moab.get_entities_by_type(0,MBVERTEX,all_nodes);
1238  // std::vector<double> aaa(all_nodes.size()*3);
1239  // CHKERR moab.tag_get_data(th,all_nodes,&aaa[0]);
1240  // CHKERR moab.tag_set_data(th_aaa,all_nodes,&aaa[0]);
1241 
1242  for (auto m : verts_map) {
1243 
1244  if (m.first < tol) {
1245 
1246  EntityHandle v = m.second.first;
1247  if (verticesOnTrimEdges.find(v) != verticesOnTrimEdges.end()) {
1248  continue;
1249  }
1250 
1251  VectorDouble3 ray_point(3);
1252  if (th) {
1253  CHKERR moab.tag_get_data(th, &v, 1, &ray_point[0]);
1254  } else {
1255  CHKERR moab.get_coords(&v, 1, &ray_point[0]);
1256  }
1257 
1258  Range adj_edges;
1259  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges);
1260  adj_edges = intersect(adj_edges, edges);
1261  Range w = intersect(adj_edges, tets_skin_edges);
1262  if (!w.empty()) {
1263  adj_edges.swap(w);
1264  }
1265  if (fixed_edges) {
1266  Range r = intersect(adj_edges, *fixed_edges);
1267  if (!r.empty()) {
1268  adj_edges.swap(r);
1269  }
1270  }
1271  if (adj_edges.empty()) {
1272  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Imposible case");
1273  }
1274 
1275  EntityHandle e = m.second.second;
1276  if (adj_edges.find(e) == adj_edges.end()) {
1277  continue;
1278  }
1279 
1280  bool corner_node = false;
1281  if (fixed_edges_verts_on_corners.find(v) !=
1282  fixed_edges_verts_on_corners.end()) {
1283  corner_node = true;
1284  }
1285 
1286  if (corner_node) {
1287 
1288  if (edgesToTrim.find(e) != edgesToTrim.end()) {
1289  auto &m = edgesToTrim.at(e);
1290  verticesOnTrimEdges[v] = m;
1291  verticesOnTrimEdges[v].dIst = 0;
1292  trimNewVertices.insert(v);
1293  }
1294 
1295  } else {
1296 
1297  VectorDouble3 new_pos(3);
1298  new_pos.clear();
1299  if (edgesToTrim.find(e) != edgesToTrim.end()) {
1300  auto &r = edgesToTrim.at(e);
1301  noalias(new_pos) = r.rayPoint + r.dIst * r.unitRayDir;
1302  VectorDouble3 unit_ray_dir = ray_point - new_pos;
1303  double dist = norm_2(unit_ray_dir);
1304  unit_ray_dir /= dist;
1305 
1306  auto get_quality_change = [this, &m_field, &moab, &new_pos, v,
1307  th](const Range &adj_tets) {
1308  double q0 = 1;
1309  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0);
1310  double q = 1;
1311  for (auto t : adj_tets) {
1312  int num_nodes;
1313  const EntityHandle *conn;
1314  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes,
1315  true);
1316  VectorDouble12 coords(12);
1317  if (th) {
1318  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1319  } else {
1320  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1321  }
1322  // cerr << coords << endl;
1323  for (int n = 0; n != 4; ++n) {
1324  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1325  if (conn[n] == v) {
1326  noalias(n_coords) = new_pos;
1327  } else {
1328  auto m = verticesOnTrimEdges.find(conn[n]);
1329  if (m != verticesOnTrimEdges.end()) {
1330  auto r = m->second;
1331  noalias(n_coords) = r.rayPoint + r.dIst * r.unitRayDir;
1332  }
1333  }
1334  }
1335  q = std::min(q, Tools::volumeLengthQuality(&coords[0]));
1336  }
1337  return q / q0;
1338  };
1339 
1340  Range adj_tets;
1341  CHKERR moab.get_adjacencies(&v, 1, 3, false, adj_tets);
1342  adj_tets = intersect(adj_tets, cutNewVolumes);
1343  double q = get_quality_change(adj_tets);
1344  if (q > 0.75) {
1345  VectorDouble3 unit_ray_dir = new_pos - ray_point;
1346  double dist = norm_2(unit_ray_dir);
1347  unit_ray_dir /= dist;
1348  verticesOnTrimEdges[v].dIst = dist;
1349  verticesOnTrimEdges[v].lEngth = dist;
1350  verticesOnTrimEdges[v].unitRayDir = unit_ray_dir;
1351  verticesOnTrimEdges[v].rayPoint = ray_point;
1352  trimNewVertices.insert(v);
1353  // CHKERR moab.add_entities(meshset, &v, 1);
1354  // CHKERR moab.add_entities(meshset, &e, 1);
1355  // CHKERR moab.tag_set_data(th_aaa, &v, 1, &new_pos[0]);
1356  }
1357  }
1358  }
1359  }
1360  }
1361 
1362  // CHKERR moab.write_file("aaaa.vtk", "VTK", "", &meshset, 1);
1363 
1364  for (auto m : verticesOnTrimEdges) {
1365  EntityHandle v = m.first;
1366  Range adj_edges;
1367  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges);
1368  adj_edges = intersect(adj_edges, edges);
1369  for (auto e : adj_edges) {
1370  edgesToTrim.erase(e);
1371  trimEdges.erase(e);
1372  }
1373  }
1374 
1376 }
map< EntityHandle, TreeData > edgesToTrim
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:32
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:495
double aveLength
Average edge length.
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Common.hpp:223
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Common.hpp:234
double tol
#define CHKERR
Inline error check.
Definition: definitions.h:614
auto getVectorAdaptor
Get Vector adaptor.
Definition: Common.hpp:256
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Common.hpp:220

◆ getCutEdges()

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

Definition at line 291 of file CutMeshInterface.hpp.

291 { return cutEdges; }

◆ getCutVolumes()

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

Definition at line 292 of file CutMeshInterface.hpp.

292 { return cutVolumes; }

◆ getMergedSurfaces()

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

Definition at line 304 of file CutMeshInterface.hpp.

◆ getMergedVolumes()

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

Definition at line 303 of file CutMeshInterface.hpp.

303 { return mergedVolumes; }

◆ getNewCutSurfaces()

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

Definition at line 294 of file CutMeshInterface.hpp.

◆ getNewCutVertices()

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

Definition at line 295 of file CutMeshInterface.hpp.

◆ getNewCutVolumes()

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

Definition at line 293 of file CutMeshInterface.hpp.

293 { return cutNewVolumes; }

◆ getNewTrimSurfaces()

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

Definition at line 300 of file CutMeshInterface.hpp.

◆ getNewTrimVertices()

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

Definition at line 301 of file CutMeshInterface.hpp.

◆ getNewTrimVolumes()

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

Definition at line 299 of file CutMeshInterface.hpp.

◆ getOptions()

MoFEMErrorCode MoFEM::CutMeshInterface::getOptions ( )

Get options from command line.

Returns
error code

Definition at line 49 of file CutMeshInterface.hpp.

49  {
51  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "MOFEM Cut mesh options",
52  "none");
53 
54  CHKERR PetscOptionsInt("-cut_linesearch_steps",
55  "number of bisection steps which line search do to "
56  "find optimal merged nodes position",
57  "", lineSearchSteps, &lineSearchSteps, PETSC_NULL);
58 
59  CHKERR PetscOptionsInt("-cut_max_merging_cycles",
60  "number of maximal merging cycles", "",
62 
63  CHKERR PetscOptionsInt(
64  "-cut_max_trim_iterations", "number of maximal merging cycles", "",
66 
67  ierr = PetscOptionsEnd();
68  CHKERRG(ierr);
70  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:562
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
#define CHKERR
Inline error check.
Definition: definitions.h:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

◆ getRayForEdge()

MoFEMErrorCode MoFEM::CutMeshInterface::getRayForEdge ( const EntityHandle  ent,
VectorAdaptor ray_point,
VectorAdaptor unit_ray_dir,
double ray_length 
) const
private

Definition at line 1472 of file CutMeshInterface.cpp.

1475  {
1476  const CoreInterface &m_field = cOre;
1477  const moab::Interface &moab = m_field.get_moab();
1479  int num_nodes;
1480  const EntityHandle *conn;
1481  CHKERR moab.get_connectivity(ent, conn, num_nodes, true);
1482  double coords[6];
1483  CHKERR moab.get_coords(conn, num_nodes, coords);
1484  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
1485  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
1486  noalias(ray_point) = s0;
1487  noalias(unit_ray_dir) = s1 - s0;
1488  ray_length = norm_2(unit_ray_dir);
1489  unit_ray_dir /= ray_length;
1491 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Common.hpp:234
#define CHKERR
Inline error check.
Definition: definitions.h:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

◆ getSurface()

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

Definition at line 289 of file CutMeshInterface.hpp.

289 { return sUrface; }

◆ getTetgenSurfaces()

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

Definition at line 306 of file CutMeshInterface.hpp.

◆ getTreeSurfPtr()

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

Definition at line 312 of file CutMeshInterface.hpp.

312  {
313  return treeSurfPtr;
314  }
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr

◆ getTrimEdges()

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

Definition at line 298 of file CutMeshInterface.hpp.

298 { return trimEdges; }

◆ getVolume()

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

Definition at line 288 of file CutMeshInterface.hpp.

288 { 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 1612 of file CutMeshInterface.cpp.

1616  {
1617  CoreInterface &m_field = cOre;
1618  moab::Interface &moab = m_field.get_moab();
1620 
1621  /**
1622  * \brief Merge nodes
1623  */
1624  struct MergeNodes {
1625  CoreInterface &mField;
1626  const bool onlyIfImproveQuality;
1627  const int lineSearch;
1628  Tag tH;
1629  bool updateMehsets;
1630 
1631  MergeNodes(CoreInterface &m_field,
1632  const bool only_if_improve_quality, const int line_search,
1633  Tag th, bool update_mehsets)
1634  : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
1635  lineSearch(line_search), tH(th), updateMehsets(update_mehsets) {
1636  mField.getInterface(nodeMergerPtr);
1637  }
1638  NodeMergerInterface *nodeMergerPtr;
1639  MoFEMErrorCode operator()(EntityHandle father, EntityHandle mother,
1640  Range &proc_tets, Range &new_surf,
1641  Range &edges_to_merge, Range &not_merged_edges,
1642  bool add_child = true) const {
1643  moab::Interface &moab = mField.get_moab();
1645  const EntityHandle conn[] = {father, mother};
1646  Range vert_tets;
1647  CHKERR moab.get_adjacencies(conn, 2, 3, false, vert_tets,
1648  moab::Interface::UNION);
1649  vert_tets = intersect(vert_tets, proc_tets);
1650  Range out_tets;
1651  CHKERR nodeMergerPtr->mergeNodes(father, mother, out_tets, &vert_tets,
1652  onlyIfImproveQuality, 0, lineSearch, tH);
1653  out_tets.merge(subtract(proc_tets, vert_tets));
1654  proc_tets.swap(out_tets);
1655 
1656  if (add_child && nodeMergerPtr->getSucessMerge()) {
1657 
1658  NodeMergerInterface::ParentChildMap &parent_child_map =
1659  nodeMergerPtr->getParentChildMap();
1660 
1661  Range child_ents;
1662  NodeMergerInterface::ParentChildMap::iterator it;
1663  for (it = parent_child_map.begin(); it != parent_child_map.end();
1664  it++) {
1665  child_ents.insert(it->pArent);
1666  }
1667 
1668  Range new_surf_child_ents = intersect(new_surf, child_ents);
1669  new_surf = subtract(new_surf, new_surf_child_ents);
1670  Range child_surf_ents;
1671  CHKERR updateRangeByChilds(parent_child_map, new_surf_child_ents,
1672  child_surf_ents);
1673  new_surf.merge(child_surf_ents);
1674 
1675  Range edges_to_merge_child_ents = intersect(edges_to_merge, child_ents);
1676  edges_to_merge = subtract(edges_to_merge, edges_to_merge_child_ents);
1677  Range merged_child_edge_ents;
1678  CHKERR updateRangeByChilds(parent_child_map, edges_to_merge_child_ents,
1679  merged_child_edge_ents);
1680 
1681  Range not_merged_edges_child_ents =
1682  intersect(not_merged_edges, child_ents);
1683  not_merged_edges =
1684  subtract(not_merged_edges, not_merged_edges_child_ents);
1685  Range not_merged_child_edge_ents;
1686  CHKERR updateRangeByChilds(parent_child_map, not_merged_edges_child_ents,
1687  not_merged_child_edge_ents);
1688 
1689  merged_child_edge_ents =
1690  subtract(merged_child_edge_ents, not_merged_child_edge_ents);
1691  edges_to_merge.merge(merged_child_edge_ents);
1692  not_merged_edges.merge(not_merged_child_edge_ents);
1693 
1694  if (updateMehsets) {
1696  (*mField.getInterface<MeshsetsManager>()), cubit_it)) {
1697  EntityHandle cubit_meshset = cubit_it->meshset;
1698  Range parent_ents;
1699  CHKERR moab.get_entities_by_handle(cubit_meshset, parent_ents, true);
1700  Range child_ents;
1701  CHKERR updateRangeByChilds(parent_child_map, parent_ents, child_ents);
1702  CHKERR moab.add_entities(cubit_meshset, child_ents);
1703  }
1704  }
1705  }
1707  }
1708 
1709  private:
1710  MoFEMErrorCode updateRangeByChilds(
1711  const NodeMergerInterface::ParentChildMap &parent_child_map,
1712  const Range &parents, Range &childs) const {
1714  NodeMergerInterface::ParentChildMap::nth_index<0>::type::iterator it;
1715  for (Range::const_iterator eit = parents.begin(); eit != parents.end();
1716  eit++) {
1717  it = parent_child_map.get<0>().find(*eit);
1718  if (it == parent_child_map.get<0>().end())
1719  continue;
1720  childs.insert(it->cHild);
1721  }
1723  }
1724  };
1725 
1726  /**
1727  * \brief Calculate edge length
1728  */
1729  struct LengthMap {
1730  Tag tH;
1731  CoreInterface &mField;
1732  moab::Interface &moab;
1733  const double maxLength;
1734  LengthMap(CoreInterface &m_field, Tag th, double max_length)
1735  : tH(th), mField(m_field), moab(m_field.get_moab()),
1736  maxLength(max_length) {
1737  gammaL = 1.;
1738  gammaQ = 3.;
1739  acceptedThrasholdMergeQuality = 0.5;
1740  }
1741 
1742  double
1743  gammaL; ///< Controls importance of length when ranking edges for merge
1744  double
1745  gammaQ; ///< Controls importance of quality when ranking edges for merge
1746  double acceptedThrasholdMergeQuality; ///< Do not merge quality if quality
1747  ///< above accepted thrashold
1748 
1749  MoFEMErrorCode operator()(const Range &tets, const Range &edges,
1750  LengthMapData_multi_index &length_map,
1751  double &ave) const {
1752  int num_nodes;
1753  const EntityHandle *conn;
1754  double coords[6];
1756  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
1757  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
1758  VectorDouble3 delta(3);
1759  for (auto edge : edges) {
1760  CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
1761  Range adj_tet;
1762  CHKERR moab.get_adjacencies(conn, num_nodes, 3, false, adj_tet);
1763  adj_tet = intersect(adj_tet, tets);
1764  if (tH) {
1765  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
1766  } else {
1767  CHKERR moab.get_coords(conn, num_nodes, coords);
1768  }
1769  double q = 1;
1770  auto abs_min = [](double a, double b) {
1771  return std::min(fabs(a), fabs(b));
1772  };
1773  CHKERR mField.getInterface<Tools>()->minTetsQuality(adj_tet, q, tH,
1774  abs_min);
1775  if (q != q)
1776  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1777  "Quality not a number");
1778  if (fabs(q) > acceptedThrasholdMergeQuality)
1779  continue;
1780  noalias(delta) = (s0 - s1) / maxLength;
1781  double dot = inner_prod(delta, delta);
1782  double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
1783  length_map.insert(LengthMapData(val, q, edge));
1784  }
1785  ave = 0;
1786  for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
1787  length_map.get<0>().begin();
1788  mit != length_map.get<0>().end(); mit++) {
1789  ave += mit->qUality;
1790  }
1791  ave /= length_map.size();
1793  }
1794  };
1795 
1796  /**
1797  * \brief Topological relations
1798  */
1799  struct Toplogy {
1800 
1801  CoreInterface &mField;
1802  Tag tH;
1803  const double tOL;
1804  Toplogy(CoreInterface &m_field, Tag th, const double tol)
1805  : mField(m_field), tH(th), tOL(tol) {}
1806 
1807  enum TYPE {
1808  FREE = 0,
1809  SKIN = 1 << 0,
1810  SURFACE = 1 << 1,
1811  SURFACE_SKIN = 1 << 2,
1812  FRONT_ENDS = 1 << 3,
1813  FIX_EDGES = 1 << 4,
1814  FIX_CORNERS = 1 << 5
1815  };
1816 
1817  typedef map<int, Range> SetsMap;
1818 
1819  MoFEMErrorCode classifyVerts(const Range &surface, const Range &tets,
1820  const Range &fixed_edges,
1821  const Range &corner_nodes,
1822  SetsMap &sets_map) const {
1823  moab::Interface &moab(mField.get_moab());
1824  Skinner skin(&moab);
1826 
1827  sets_map[FIX_CORNERS].merge(corner_nodes);
1828  Range fixed_verts;
1829  CHKERR moab.get_connectivity(fixed_edges, fixed_verts, true);
1830  sets_map[FIX_EDGES].swap(fixed_verts);
1831 
1832  Range tets_skin;
1833  CHKERR skin.find_skin(0, tets, false, tets_skin);
1834  Range tets_skin_edges;
1835  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1836  moab::Interface::UNION);
1837 
1838  // surface skin
1839  Range surface_skin;
1840  CHKERR skin.find_skin(0, surface, false, surface_skin);
1841  Range front_in_the_body;
1842  front_in_the_body = subtract(surface_skin, tets_skin_edges);
1843  Range front_ends;
1844  CHKERR skin.find_skin(0, front_in_the_body, false, front_ends);
1845  sets_map[FRONT_ENDS].swap(front_ends);
1846 
1847  Range surface_skin_verts;
1848  CHKERR moab.get_connectivity(surface_skin, surface_skin_verts, true);
1849  sets_map[SURFACE_SKIN].swap(surface_skin_verts);
1850 
1851  // surface
1852  Range surface_verts;
1853  CHKERR moab.get_connectivity(surface, surface_verts, true);
1854  sets_map[SURFACE].swap(surface_verts);
1855 
1856  // skin
1857  Range tets_skin_verts;
1858  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1859  sets_map[SKIN].swap(tets_skin_verts);
1860 
1861  Range tets_verts;
1862  CHKERR moab.get_connectivity(tets, tets_verts, true);
1863  sets_map[FREE].swap(tets_verts);
1864 
1866  }
1867 
1868  MoFEMErrorCode getProcTets(const Range &tets, const Range &edges_to_merge,
1869  Range &proc_tets) const {
1870  moab::Interface &moab(mField.get_moab());
1872  Range edges_to_merge_verts;
1873  CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts, true);
1874  Range edges_to_merge_verts_tets;
1875  CHKERR moab.get_adjacencies(edges_to_merge_verts, 3, false,
1876  edges_to_merge_verts_tets,
1877  moab::Interface::UNION);
1878  edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
1879  proc_tets.swap(edges_to_merge_verts_tets);
1881  }
1882 
1883  MoFEMErrorCode edgesToMerge(const Range &surface, const Range &tets,
1884  Range &edges_to_merge) const {
1885  moab::Interface &moab(mField.get_moab());
1887 
1888  Range surface_verts;
1889  CHKERR moab.get_connectivity(surface, surface_verts, true);
1890  Range surface_verts_edges;
1891  CHKERR moab.get_adjacencies(surface_verts, 1, false, surface_verts_edges,
1892  moab::Interface::UNION);
1893  edges_to_merge.merge(surface_verts_edges);
1894  Range tets_edges;
1895  CHKERR moab.get_adjacencies(tets, 1, false, tets_edges,
1896  moab::Interface::UNION);
1897  edges_to_merge = intersect(edges_to_merge, tets_edges);
1899  }
1900 
1901  MoFEMErrorCode removeBadEdges(const Range &surface, const Range &tets,
1902  const Range &fixed_edges,
1903  const Range &corner_nodes,
1904  Range &edges_to_merge,
1905  Range &not_merged_edges) {
1906  moab::Interface &moab(mField.get_moab());
1908 
1909  // find skin
1910  Skinner skin(&moab);
1911  Range tets_skin;
1912  CHKERR skin.find_skin(0, tets, false, tets_skin);
1913  Range surface_skin;
1914  CHKERR skin.find_skin(0, surface, false, surface_skin);
1915 
1916  // end nodes
1917  Range tets_skin_edges;
1918  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1919  moab::Interface::UNION);
1920  Range surface_front;
1921  surface_front = subtract(surface_skin, tets_skin_edges);
1922  Range surface_front_nodes;
1923  CHKERR moab.get_connectivity(surface_front, surface_front_nodes, true);
1924  Range ends_nodes;
1925  CHKERR skin.find_skin(0, surface_front, false, ends_nodes);
1926 
1927  // remove bad merges
1928 
1929  // get surface and body skin verts
1930  Range surface_edges;
1931  CHKERR moab.get_adjacencies(surface, 1, false, surface_edges,
1932  moab::Interface::UNION);
1933  // get nodes on the surface
1934  Range surface_edges_verts;
1935  CHKERR moab.get_connectivity(surface_edges, surface_edges_verts, true);
1936  // get vertices on the body skin
1937  Range tets_skin_edges_verts;
1938  CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts, true);
1939 
1940  Range edges_to_remove;
1941 
1942  // remove edges self connected to body skin
1943  {
1944  Range ents_nodes_and_edges;
1945  ents_nodes_and_edges.merge(tets_skin_edges_verts);
1946  ents_nodes_and_edges.merge(tets_skin_edges);
1947  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
1948  0, false);
1949  }
1950  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1951  not_merged_edges.merge(edges_to_remove);
1952 
1953  // remove edges self connected to surface
1954  {
1955  Range ents_nodes_and_edges;
1956  ents_nodes_and_edges.merge(surface_edges_verts);
1957  ents_nodes_and_edges.merge(surface_edges);
1958  ents_nodes_and_edges.merge(tets_skin_edges_verts);
1959  ents_nodes_and_edges.merge(tets_skin_edges);
1960  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
1961  0, false);
1962  }
1963  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1964  not_merged_edges.merge(edges_to_remove);
1965 
1966  // remove edges adjacent corner_nodes execpt those on fixed edges
1967  Range fixed_edges_nodes;
1968  CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes, true);
1969  {
1970  Range ents_nodes_and_edges;
1971  ents_nodes_and_edges.merge(fixed_edges_nodes);
1972  ents_nodes_and_edges.merge(ends_nodes);
1973  ents_nodes_and_edges.merge(corner_nodes);
1974  ents_nodes_and_edges.merge(fixed_edges);
1975  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
1976  0, false);
1977  }
1978  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1979  not_merged_edges.merge(edges_to_remove);
1980 
1981  // remove edges self connected to surface
1982  CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove, 0, false);
1983  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1984  not_merged_edges.merge(edges_to_remove);
1985 
1986  // remove edges self contented on surface skin
1987  {
1988  Range ents_nodes_and_edges;
1989  ents_nodes_and_edges.merge(surface_skin);
1990  ents_nodes_and_edges.merge(fixed_edges_nodes);
1991  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
1992  0, false);
1993  }
1994  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1995  not_merged_edges.merge(edges_to_remove);
1996 
1997  // remove crack front nodes connected to the surface
1998  {
1999  Range ents_nodes_and_edges;
2000  ents_nodes_and_edges.merge(surface_front_nodes);
2001  ents_nodes_and_edges.merge(surface_front);
2002  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2003  ents_nodes_and_edges.merge(tets_skin_edges);
2004  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2005  0, false);
2006  }
2007  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2008  not_merged_edges.merge(edges_to_remove);
2009 
2010  // remove edges connecting crack front and fixed edges, except those
2011  {
2012  Range ents_nodes_and_edges;
2013  ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
2014  ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2015  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2016  tOL, false);
2017  }
2018  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2019  not_merged_edges.merge(edges_to_remove);
2020 
2021 
2023  }
2024 
2025  private:
2026  MoFEMErrorCode removeSelfConectingEdges(const Range &ents,
2027  Range &edges_to_remove,
2028  const bool length,
2029  bool debug) const {
2030  moab::Interface &moab(mField.get_moab());
2032  // get nodes
2033  Range ents_nodes = ents.subset_by_type(MBVERTEX);
2034  if (ents_nodes.empty()) {
2035  CHKERR moab.get_connectivity(ents, ents_nodes, true);
2036  }
2037  // edges adj. to nodes
2038  Range ents_nodes_edges;
2039  CHKERR moab.get_adjacencies(ents_nodes, 1, false, ents_nodes_edges,
2040  moab::Interface::UNION);
2041  // nodes of adj. edges
2042  Range ents_nodes_edges_nodes;
2043  CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2044  true);
2045  // hanging nodes
2046  ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2047  Range ents_nodes_edges_nodes_edges;
2048  CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1, false,
2049  ents_nodes_edges_nodes_edges,
2050  moab::Interface::UNION);
2051  // remove edges adj. to hanging edges
2052  ents_nodes_edges =
2053  subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2054  ents_nodes_edges =
2055  subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2056  if(length>0) {
2057  Range::iterator eit = ents_nodes_edges.begin();
2058  for (; eit != ents_nodes_edges.end();) {
2059 
2060  int num_nodes;
2061  const EntityHandle *conn;
2062  rval = moab.get_connectivity(*eit, conn, num_nodes, true);
2063  double coords[6];
2064  if(tH) {
2065  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
2066  } else {
2067  CHKERR moab.get_coords(conn, num_nodes, coords);
2068  }
2069 
2070  auto get_edge_length = [coords]() {
2072  &coords[0], &coords[1], &coords[2]);
2075  t_delta(i) = t_coords(i);
2076  ++t_coords;
2077  t_delta(i) -= t_coords(i);
2078  return sqrt(t_delta(i) * t_delta(i));
2079  };
2080 
2081  if (get_edge_length() < tOL) {
2082  eit = ents_nodes_edges.erase(eit);
2083  } else {
2084  eit++;
2085  }
2086  }
2087  }
2088  edges_to_remove.swap(ents_nodes_edges);
2089  if (debug) {
2090  EntityHandle meshset;
2091  CHKERR moab.create_meshset(MESHSET_SET, meshset);
2092  CHKERR moab.add_entities(meshset, edges_to_remove);
2093  CHKERR moab.write_file("edges_to_remove.vtk", "VTK", "", &meshset, 1);
2094  CHKERR moab.delete_entities(&meshset, 1);
2095  }
2097  }
2098  };
2099 
2100  Range not_merged_edges;
2101  const double tol = 0.05;
2102  CHKERR Toplogy(m_field, th, tol * aveLength)
2103  .edgesToMerge(surface, tets, edges_to_merge);
2104  CHKERR Toplogy(m_field, th, tol * aveLength)
2105  .removeBadEdges(surface, tets, fixed_edges, corner_nodes, edges_to_merge,
2106  not_merged_edges);
2107  Toplogy::SetsMap sets_map;
2108  CHKERR Toplogy(m_field, th, tol * aveLength)
2109  .classifyVerts(surface, tets, fixed_edges, corner_nodes, sets_map);
2110  Range proc_tets;
2111  CHKERR Toplogy(m_field, th, tol * aveLength)
2112  .getProcTets(tets, edges_to_merge, proc_tets);
2113  out_tets = subtract(tets, proc_tets);
2114  if (bit_ptr) {
2115  for (int dd = 2; dd >= 0; dd--) {
2116  CHKERR moab.get_adjacencies(out_tets.subset_by_dimension(3), dd, false,
2117  out_tets, moab::Interface::UNION);
2118  }
2119  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(out_tets,
2120  *bit_ptr);
2121  }
2122 
2123  int nb_nodes_merged = 0;
2124  LengthMapData_multi_index length_map;
2125  new_surf = surface;
2126 
2127  double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2128  for (int pp = 0; pp != nbMaxMergingCycles; ++pp) {
2129 
2130  int nb_nodes_merged_p = nb_nodes_merged;
2131  length_map.clear();
2132  min_pp = min_p;
2133  min_p = min;
2134  CHKERR LengthMap(m_field, th, aveLength)(proc_tets, edges_to_merge,
2135  length_map, ave);
2136  min = length_map.get<0>().begin()->qUality;
2137  if(pp == 0) {
2138  ave0 = ave;
2139  }
2140 
2141  int nn = 0;
2142  Range collapsed_edges;
2143  for (LengthMapData_multi_index::nth_index<0>::type::iterator
2144  mit = length_map.get<0>().begin();
2145  mit != length_map.get<0>().end(); mit++, nn++) {
2146  // cerr << mit->lEngth << endl; //" " << mit->eDge << endl;
2147  if(mit->skip) continue;
2148  int num_nodes;
2149  const EntityHandle *conn;
2150  CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes, true);
2151  int conn_type[2] = {0, 0};
2152  for (int nn = 0; nn != 2; nn++) {
2153  conn_type[nn] = 0;
2154  for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2155  sit != sets_map.rend(); sit++) {
2156  if (sit->second.find(conn[nn]) != sit->second.end()) {
2157  conn_type[nn] |= sit->first;
2158  }
2159  }
2160  }
2161  int type_father, type_mother;
2162  EntityHandle father, mother;
2163  if (conn_type[0] > conn_type[1]) {
2164  father = conn[0];
2165  mother = conn[1];
2166  type_father = conn_type[0];
2167  type_mother = conn_type[1];
2168  } else {
2169  father = conn[1];
2170  mother = conn[0];
2171  type_father = conn_type[1];
2172  type_mother = conn_type[0];
2173  }
2174  int line_search = 0;
2175  if (type_father == type_mother) {
2176  line_search = lineSearchSteps;
2177  }
2178 
2179  CHKERR MergeNodes(m_field, true, line_search, th,
2180  update_meshsets)(father, mother, proc_tets, new_surf,
2181  edges_to_merge, not_merged_edges);
2182 
2183  if (m_field.getInterface<NodeMergerInterface>()->getSucessMerge()) {
2184  Range adj_mother_tets;
2185  CHKERR moab.get_adjacencies(&mother, 1, 3, false, adj_mother_tets);
2186  Range adj_mother_tets_nodes;
2187  CHKERR moab.get_connectivity(adj_mother_tets, adj_mother_tets_nodes,
2188  true);
2189  Range adj_edges;
2190  CHKERR moab.get_adjacencies(adj_mother_tets_nodes, 1, false, adj_edges,
2191  moab::Interface::UNION);
2192  CHKERR moab.get_adjacencies(&father, 1, 1, false, adj_edges,
2193  moab::Interface::UNION);
2194  for (Range::iterator ait = adj_edges.begin(); ait != adj_edges.end();
2195  ait++) {
2196  LengthMapData_multi_index::nth_index<1>::type::iterator miit =
2197  length_map.get<1>().find(*ait);
2198  if (miit != length_map.get<1>().end()) {
2199  (const_cast<LengthMapData &>(*miit)).skip = true;
2200  }
2201  }
2202  nb_nodes_merged++;
2203  collapsed_edges.insert(mit->eDge);
2204  }
2205 
2206  if (nn > static_cast<int>(length_map.size() / fraction_level))
2207  break;
2208  if (mit->qUality > ave)
2209  break;
2210  }
2211 
2212  Range adj_faces, adj_edges;
2213  CHKERR moab.get_adjacencies(proc_tets, 2, false, adj_faces,
2214  moab::Interface::UNION);
2215  new_surf = intersect(new_surf, adj_faces);
2216 
2217  PetscPrintf(m_field.get_comm(),
2218  "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e\n", pp,
2219  nb_nodes_merged, ave, min);
2220 
2221  if (debug) {
2222  // current surface and merged edges in step
2223  EntityHandle meshset;
2224  std::string name;
2225  CHKERR moab.create_meshset(MESHSET_SET, meshset);
2226  CHKERR moab.add_entities(meshset, new_surf);
2227  CHKERR moab.add_entities(meshset, collapsed_edges);
2228  name = "node_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2229  CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1);
2230  CHKERR moab.delete_entities(&meshset, 1);
2231  // edges to merge
2232  CHKERR moab.create_meshset(MESHSET_SET, meshset);
2233  CHKERR moab.add_entities(meshset, edges_to_merge);
2234  name = "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) +
2235  ".vtk";
2236  CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1);
2237  CHKERR moab.delete_entities(&meshset, 1);
2238  }
2239 
2240  if (nb_nodes_merged == nb_nodes_merged_p)
2241  break;
2242  if (min > 1e-2 && min == min_pp)
2243  break;
2244  if (min > ave0)
2245  break;
2246 
2247  CHKERR moab.get_adjacencies(proc_tets, 1, false, adj_edges,
2248  moab::Interface::UNION);
2249  edges_to_merge = intersect(edges_to_merge, adj_edges);
2250  CHKERR Toplogy(m_field, th, tol * aveLength)
2251  .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
2252  edges_to_merge, not_merged_edges);
2253  }
2254 
2255  if (bit_ptr) {
2256  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(proc_tets,
2257  *bit_ptr);
2258  }
2259  out_tets.merge(proc_tets);
2260 
2262 }
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Common.hpp:60
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:519
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
double aveLength
Average edge length.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:526
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Common.hpp:234
double tol
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
static const bool debug
multi_index_container< LengthMapData, indexed_by< ordered_non_unique< member< LengthMapData, double, &LengthMapData::lEngth > >, hashed_unique< member< LengthMapData, EntityHandle, &LengthMapData::eDge > > >> LengthMapData_multi_index
#define CHKERR
Inline error check.
Definition: definitions.h:614
#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:439
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:150
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Common.hpp:220

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

2268  {
2269  CoreInterface &m_field = cOre;
2271  Range tets_level;
2272  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2273  trim_bit, BitRefLevel().set(), MBTET, tets_level);
2274  CHKERRG(ierr);
2275 
2276  Range edges_to_merge;
2277  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByParentType(
2278  trim_bit, trim_bit | cut_bit, MBEDGE, edges_to_merge);
2279  CHKERRG(ierr);
2280  edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
2281 
2282  // get all entities not in database
2283  Range all_ents_not_in_database_before;
2284  ierr = cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2285  all_ents_not_in_database_before);
2286  CHKERRG(ierr);
2287 
2288  Range out_new_tets, out_new_surf;
2289  ierr = mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
2290  corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
2291  th, update_meshsets, &bit, debug);
2292  CHKERRG(ierr);
2293 
2294  // get all entities not in database after merge
2295  Range all_ents_not_in_database_after;
2296  ierr = cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2297  all_ents_not_in_database_after);
2298  CHKERRG(ierr);
2299  // delete hanging entities
2300  all_ents_not_in_database_after =
2301  subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
2302  Range meshsets;
2303  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
2304  false);
2305  for (Range::iterator mit = meshsets.begin(); mit != meshsets.end(); mit++) {
2306  CHKERR m_field.get_moab().remove_entities(*mit,
2307  all_ents_not_in_database_after);
2308  }
2309  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
2310 
2311  mergedVolumes.swap(out_new_tets);
2312  mergedSurfaces.swap(out_new_surf);
2314 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:519
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.
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:562
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:526
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
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:614

◆ mergeSurface()

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

merge surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 123 of file CutMeshInterface.cpp.

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

◆ mergeVolumes()

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

merge volume entities

Parameters
volumeentities which going to be added
Returns
error code

Definition at line 129 of file CutMeshInterface.cpp.

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

◆ moveMidNodesOnCutEdges()

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

projecting of mid edge nodes on new mesh on surface

Returns
error code

Definition at line 1012 of file CutMeshInterface.cpp.

1012  {
1014 
1015  CoreInterface &m_field = cOre;
1016  moab::Interface &moab = m_field.get_moab();
1018 
1019  // Range out_side_vertices;
1020  for (map<EntityHandle, TreeData>::iterator mit = verticesOnCutEdges.begin();
1021  mit != verticesOnCutEdges.end(); mit++) {
1022  double dist = mit->second.dIst;
1023  VectorDouble3 new_coors =
1024  mit->second.rayPoint + dist * mit->second.unitRayDir;
1025  if (th) {
1026  CHKERR moab.tag_set_data(th, &mit->first, 1, &new_coors[0]);
1027  } else {
1028  CHKERR moab.set_coords(&mit->first, 1, &new_coors[0]);
1029  }
1030  }
1031 
1033 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:519
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
#define CHKERR
Inline error check.
Definition: definitions.h:614
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:439
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Common.hpp:220

◆ moveMidNodesOnTrimmedEdges()

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

move trimmed edges mid nodes

Returns
error code

Definition at line 1035 of file CutMeshInterface.cpp.

1035  {
1036  CoreInterface &m_field = cOre;
1037  moab::Interface &moab = m_field.get_moab();
1039  // Range out_side_vertices;
1040  for (map<EntityHandle, TreeData>::iterator mit = verticesOnTrimEdges.begin();
1041  mit != verticesOnTrimEdges.end(); mit++) {
1042  double dist = mit->second.dIst;
1043  // cout << s << " " << mit->second.dIst << " " << mit->second.lEngth <<
1044  // endl;
1045  VectorDouble3 new_coors =
1046  mit->second.rayPoint + dist * mit->second.unitRayDir;
1047  if (th) {
1048  CHKERR moab.tag_set_data(th, &mit->first, 1, &new_coors[0]);
1049  } else {
1050  CHKERR moab.set_coords(&mit->first, 1, &new_coors[0]);
1051  }
1052  }
1054 }
map< EntityHandle, TreeData > verticesOnTrimEdges
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
#define CHKERR
Inline error check.
Definition: definitions.h:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Common.hpp:220

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

640  {
641  CoreInterface &m_field = cOre;
642  moab::Interface &moab = m_field.get_moab();
644 
645  // Get entities on body skin
646  Skinner skin(&moab);
647  Range tets_skin;
648  rval = skin.find_skin(0, vOlume, false, tets_skin);
649  Range tets_skin_edges;
650  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
651  moab::Interface::UNION);
652  Range tets_skin_verts;
653  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
654 
655  // Get entities in volume
656  Range vol_faces, vol_edges, vol_nodes;
657  CHKERR moab.get_adjacencies(vOlume, 2, false, vol_faces,
658  moab::Interface::UNION);
659  CHKERR moab.get_adjacencies(vOlume, 1, false, vol_edges,
660  moab::Interface::UNION);
661  CHKERR moab.get_connectivity(vOlume, vol_nodes, true);
662 
663  // Get nodes on cut edges
664  Range cut_edge_verts;
665  CHKERR moab.get_connectivity(cutEdges, cut_edge_verts, true);
666 
667  Range fixed_edges_nodes;
668  if(fixed_edges) {
669  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_nodes, true);
670  }
671 
672  // Get faces and edges
673  Range cut_edges_faces;
674  CHKERR moab.get_adjacencies(cut_edge_verts, 2, true, cut_edges_faces,
675  moab::Interface::UNION);
676  cut_edges_faces = intersect(cut_edges_faces, vol_faces);
677  Range cut_edges_faces_verts;
678  CHKERR moab.get_connectivity(cut_edges_faces, cut_edges_faces_verts, true);
679  cut_edges_faces_verts = subtract(cut_edges_faces_verts, cut_edge_verts);
680  Range to_remove_cut_edges_faces;
681  CHKERR moab.get_adjacencies(cut_edges_faces_verts, 2, true,
682  to_remove_cut_edges_faces,
683  moab::Interface::UNION);
684  cut_edges_faces = subtract(cut_edges_faces, to_remove_cut_edges_faces);
685  cut_edges_faces.merge(cutEdges);
686 
687  Tag th_dist_normal;
688  CHKERR moab.tag_get_handle("DIST_NORMAL", th_dist_normal);
689 
690  auto get_quality_change =
691  [this, &m_field,
692  &moab](const Range &adj_tets,
693  map<EntityHandle, TreeData> vertices_on_cut_edges) {
694  vertices_on_cut_edges.insert(verticesOnCutEdges.begin(),
695  verticesOnCutEdges.end());
696  double q0 = 1;
697  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0);
698  double q = 1;
699  for (auto t : adj_tets) {
700  int num_nodes;
701  const EntityHandle *conn;
702  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
703  VectorDouble12 coords(12);
704  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
705  // cerr << coords << endl;
706  for (int n = 0; n != 4; ++n) {
707  bool ray_found = false;
708  auto mit = vertices_on_cut_edges.find(conn[n]);
709  if (mit != vertices_on_cut_edges.end()) {
710  ray_found = true;
711  }
712  if (ray_found) {
713  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
714  double dist = mit->second.dIst;
715  noalias(n_coords) =
716  mit->second.rayPoint + dist * mit->second.unitRayDir;
717  }
718  }
719  q = std::min(q, Tools::volumeLengthQuality(&coords[0]));
720  }
721  if (std::isnormal(q))
722  return q / q0;
723  else
724  return -1.;
725  };
726 
727  auto get_conn = [&moab](const EntityHandle &e,
728  const EntityHandle *&conn, int &num_nodes) {
730  EntityType type = moab.type_from_handle(e);
731  if (type == MBVERTEX) {
732  conn = &e;
733  num_nodes = 1;
734  } else {
735  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
736  }
738  };
739 
740  auto get_normal_dist = [](const double *normal) {
741  FTensor::Tensor1<double, 3> t_n(normal[0], normal[1], normal[2]);
743  return sqrt(t_n(i) * t_n(i));
744  };
745 
746  auto get_normal_dist_from_conn = [&moab, get_normal_dist,
747  th_dist_normal](EntityHandle v) {
748  double dist_normal[3];
749  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, dist_normal);
750  return get_normal_dist(dist_normal);
751  };
752 
753  auto project_node = [&moab, th_dist_normal](
754  const EntityHandle v,
755  map<EntityHandle, TreeData> &vertices_on_cut_edges) {
757  VectorDouble3 dist_normal(3);
758  rval = moab.tag_get_data(th_dist_normal, &v, 1, &dist_normal[0]);
759  VectorDouble3 s0(3);
760  CHKERR moab.get_coords(&v, 1, &s0[0]);
761  double dist = norm_2(dist_normal);
762  vertices_on_cut_edges[v].dIst = dist;
763  vertices_on_cut_edges[v].lEngth = dist;
764  vertices_on_cut_edges[v].unitRayDir =
765  dist > 0 ? dist_normal / dist : dist_normal;
766  vertices_on_cut_edges[v].rayPoint = s0;
768  };
769 
770  for (int d = 2; d >= 0; --d) {
771 
772  Range ents;
773  if (d > 0)
774  ents = cut_edges_faces.subset_by_dimension(d);
775  else
776  ents = cut_edge_verts;
777 
778  // make list of entities
779  multimap<double, EntityHandle> ents_to_check;
780  for (auto f : ents) {
781  int num_nodes;
782  const EntityHandle *conn;
783  CHKERR get_conn(f, conn, num_nodes);
784  VectorDouble3 dist(3);
785  for (int n = 0; n != num_nodes; ++n) {
786  dist[n] = get_normal_dist_from_conn(conn[n]);
787  }
788  double max_dist = 0;
789  for (int n = 0; n != num_nodes; ++n) {
790  max_dist = std::max(max_dist, fabs(dist[n]));
791  }
792  if (max_dist < low_tol * get_ave_edge_length(f, vol_edges, moab)) {
793  ents_to_check.insert(std::pair<double, EntityHandle>(max_dist, f));
794  }
795  }
796 
797  double ray_point[3], unit_ray_dir[3];
798  VectorAdaptor vec_unit_ray_dir(
799  3, ublas::shallow_array_adaptor<double>(3, unit_ray_dir));
800  VectorAdaptor vec_ray_point(
801  3, ublas::shallow_array_adaptor<double>(3, ray_point));
802 
803  for (auto m : ents_to_check) {
804 
805  EntityHandle f = m.second;
806 
807  int num_nodes;
808  const EntityHandle *conn;
809  CHKERR get_conn(f, conn, num_nodes);
810  VectorDouble9 coords(9);
811  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
812 
813  Range adj_tets;
814  CHKERR moab.get_adjacencies(conn, num_nodes, 3, false, adj_tets,
815  moab::Interface::UNION);
816  adj_tets = intersect(adj_tets, vOlume);
817 
818  map<EntityHandle, TreeData> vertices_on_cut_edges;
819  for (int n = 0; n != num_nodes; ++n) {
820  const EntityHandle node = conn[n];
821  CHKERR project_node(node, vertices_on_cut_edges);
822  }
823  if (static_cast<int>(vertices_on_cut_edges.size()) != num_nodes) {
824  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
825  "Data inconsistency");
826  }
827 
828  double q = get_quality_change(adj_tets, vertices_on_cut_edges);
829  if (q > 0.75) {
830  bool check_q_again = false;
831  for (auto &m : vertices_on_cut_edges) {
832  EntityHandle node = m.first;
833  if (tets_skin_verts.find(node) != tets_skin_verts.end()) {
834 
835  check_q_again = true;
836 
837  // check if node is at the corner
838  bool zero_disp_node = false;
839  if (corner_nodes) {
840  if (corner_nodes->find(node) != corner_nodes->end()) {
841  zero_disp_node = true;
842  }
843  }
844 
845  // check node is on the fixed edge
846  Range adj_edges;
847  CHKERR moab.get_adjacencies(&node, 1, 1, false, adj_edges);
848  adj_edges = intersect(adj_edges, tets_skin_edges);
849  if (fixed_edges) {
850  Range e;
851  // check if node is on fixed edge
852  e = intersect(adj_edges, *fixed_edges);
853  if (!e.empty()) {
854  adj_edges.swap(e);
855  }
856  // check if split edge is fixed edge
857  e = intersect(adj_edges, cutEdges);
858  if (!e.empty()) {
859  adj_edges.swap(e);
860  } else {
861  zero_disp_node = true;
862  }
863  }
864 
865  VectorDouble3 s0(3);
866  CHKERR moab.get_coords(&node, 1, &s0[0]);
867 
868  if (zero_disp_node) {
869  VectorDouble3 z(3);
870  z.clear();
871  m.second.dIst = 0;
872  m.second.lEngth = 0;
873  m.second.unitRayDir = z;
874  m.second.rayPoint = s0;
875  } else {
876  if (adj_edges.empty()) {
877  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
878  "Data inconsistency");
879  }
880  for (auto e : adj_edges) {
881  if (edgesToCut.find(e) != edgesToCut.end()) {
882  auto d = edgesToCut.at(e);
883  VectorDouble3 new_pos = d.rayPoint + d.dIst * d.unitRayDir;
884  VectorDouble3 ray = new_pos - s0;
885  double dist0 = norm_2(ray);
886  m.second.dIst = dist0;
887  m.second.lEngth = dist0;
888  m.second.unitRayDir = dist0 > 0 ? ray / dist0 : ray;
889  m.second.rayPoint = s0;
890  break;
891  } else {
892  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
893  "Data inconsistency");
894  }
895  }
896  }
897  }
898  }
899 
900  if (check_q_again) {
901  q = get_quality_change(adj_tets, vertices_on_cut_edges);
902  }
903  if (q > 0.75) {
904  verticesOnCutEdges.insert(vertices_on_cut_edges.begin(),
905  vertices_on_cut_edges.end());
906  EntityHandle type = moab.type_from_handle(f);
907  if (type == MBVERTEX) {
908  zeroDistanceVerts.insert(f);
909  } else {
910  zeroDistanceEnts.insert(f);
911  }
912  }
913  }
914  }
915 
916  }
917 
918  for (auto f : unite(zeroDistanceEnts, zeroDistanceVerts)) {
919  int num_nodes;
920  const EntityHandle *conn;
921  CHKERR get_conn(f, conn, num_nodes);
922  Range adj_edges;
923  CHKERR moab.get_adjacencies(conn, num_nodes, 1, false, adj_edges,
924  moab::Interface::UNION);
925  for (auto e : adj_edges) {
926  cutEdges.erase(e);
927  edgesToCut.erase(e);
928  }
929  }
930 
931  if (debug) {
932  EntityHandle out_meshset;
933  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, out_meshset);
934  CHKERR m_field.get_moab().add_entities(out_meshset, zeroDistanceEnts);
935  CHKERR m_field.get_moab().add_entities(out_meshset, zeroDistanceVerts);
936  CHKERR m_field.get_moab().write_file("projected_zero_distance_ents.vtk",
937  "VTK", "", &out_meshset, 1);
938  CHKERR m_field.get_moab().delete_entities(&out_meshset, 1);
939  }
940 
942 }
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:32
VectorBoundedArray< double, 9 > VectorDouble9
Definition: Common.hpp:222
static double get_ave_edge_length(const EntityHandle ent, const Range &vol_edges, moab::Interface &moab)
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:495
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Common.hpp:223
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Common.hpp:234
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:614
auto getVectorAdaptor
Get Vector adaptor.
Definition: Common.hpp:256
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:439
map< EntityHandle, TreeData > edgesToCut
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Common.hpp:220

◆ projectZeroDistanceEnts() [2/2]

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

Definition at line 296 of file CutMeshInterface.hpp.

◆ query_interface()

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

Implements MoFEM::UnknownInterface.

Definition at line 21 of file CutMeshInterface.cpp.

22  {
24  *iface = NULL;
25  if (uuid == IDD_MOFEMCutMesh) {
26  *iface = const_cast<CutMeshInterface *>(this);
28  }
29  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
31 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:519
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:526
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 2318 of file CutMeshInterface.cpp.

2321  {
2322  CoreInterface &m_field = cOre;
2323  moab::Interface &moab = m_field.get_moab();
2324  TetGenInterface *tetgen_iface;
2326  CHKERR m_field.getInterface(tetgen_iface);
2327 
2328  tetGenData.clear();
2329  moabTetGenMap.clear();
2330  tetGenMoabMap.clear();
2331 
2332  if (tetGenData.size() < 1) {
2333  tetGenData.push_back(new tetgenio);
2334  }
2335  tetgenio &in = tetGenData.back();
2336 
2337  struct BitEnts {
2338 
2339  CoreInterface &mField;
2340  const BitRefLevel &bIt;
2341  BitEnts(CoreInterface &m_field, const BitRefLevel &bit)
2342  : mField(m_field), bIt(bit) {}
2343 
2344  Range mTets;
2345  Range mTris;
2346  Range mEdges;
2347  Range mNodes;
2348 
2349  MoFEMErrorCode getLevelEnts() {
2351  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2352  bIt, BitRefLevel().set(), MBTET, mTets);
2353  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2354  bIt, BitRefLevel().set(), MBTRI, mTris);
2355  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2356  bIt, BitRefLevel().set(), MBEDGE, mEdges);
2357  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2358  bIt, BitRefLevel().set(), MBVERTEX, mNodes);
2360  }
2361 
2362  Range mSkin;
2363  Range mSkinNodes;
2364  Range mSkinEdges;
2365 
2366  MoFEMErrorCode getSkin() {
2367  moab::Interface &moab = mField.get_moab();
2369  Skinner skin(&moab);
2370  CHKERR skin.find_skin(0, mTets, false, mSkin);
2371  CHKERR mField.get_moab().get_connectivity(mSkin, mSkinNodes, true);
2372  CHKERR mField.get_moab().get_adjacencies(mSkin, 1, false, mSkinEdges,
2373  moab::Interface::UNION);
2375  }
2376  };
2377 
2378  struct SurfaceEnts {
2379 
2380  CoreInterface &mField;
2381  SurfaceEnts(CoreInterface &m_field) : mField(m_field) {}
2382 
2383  Range sNodes;
2384  Range sEdges;
2385  Range sVols;
2386  Range vNodes;
2387 
2388  MoFEMErrorCode getVolume(const BitEnts &bit_ents, const Range &tris) {
2389  moab::Interface &moab = mField.get_moab();
2391  CHKERR moab.get_connectivity(tris, sNodes, true);
2392  CHKERR moab.get_adjacencies(tris, 1, false, sEdges,
2393  moab::Interface::UNION);
2394  CHKERR moab.get_adjacencies(sNodes, 3, false, sVols,
2395  moab::Interface::UNION);
2396  sVols = intersect(sVols, bit_ents.mTets);
2397  CHKERR moab.get_connectivity(sVols, vNodes, true);
2399  }
2400 
2401  Range sSkin;
2402  Range sSkinNodes;
2403  Range vSkin;
2404  Range vSkinNodes;
2405  Range vSkinWithoutBodySkin;
2406  Range vSkinNodesWithoutBodySkin;
2407  Range vSkinOnBodySkin;
2408  Range vSkinOnBodySkinNodes;
2409 
2410  MoFEMErrorCode getSkin(const BitEnts &bit_ents, const Range &tris,
2411  const int levels = 3) {
2412  moab::Interface &moab = mField.get_moab();
2414  Skinner skin(&moab);
2415  rval = skin.find_skin(0, sVols, false, vSkin);
2416  for (int ll = 0; ll != levels; ll++) {
2417  CHKERR moab.get_adjacencies(vSkin, 3, false, sVols,
2418  moab::Interface::UNION);
2419  sVols = intersect(sVols, bit_ents.mTets);
2420  vSkin.clear();
2421  CHKERR skin.find_skin(0, sVols, false, vSkin);
2422  }
2423  vSkinWithoutBodySkin = subtract(vSkin, bit_ents.mSkin);
2424  vSkinOnBodySkin = intersect(vSkin, bit_ents.mSkin);
2425  CHKERR moab.get_connectivity(vSkinOnBodySkin, vSkinOnBodySkinNodes, true);
2426  CHKERR moab.get_connectivity(sVols, vNodes, true);
2427  CHKERR moab.get_connectivity(vSkin, vSkinNodes, true);
2428  vSkinNodesWithoutBodySkin = subtract(vSkinNodes, bit_ents.mSkinNodes);
2429  CHKERR skin.find_skin(0, tris, false, sSkin);
2430  CHKERR moab.get_connectivity(sSkin, sSkinNodes, true);
2431  tVols = sVols;
2433  }
2434 
2435  Range tVols;
2436 
2437  MoFEMErrorCode getTetsForRemesh(const BitEnts &bit_ents, Tag th = NULL) {
2438  moab::Interface &moab = mField.get_moab();
2440 
2441  Range tets_with_four_nodes_on_skin;
2442  rval = moab.get_adjacencies(vSkinOnBodySkinNodes, 3, false,
2443  tets_with_four_nodes_on_skin,
2444  moab::Interface::UNION);
2445  Range tets_nodes;
2446  CHKERR moab.get_connectivity(tets_with_four_nodes_on_skin, tets_nodes,
2447  true);
2448  tets_nodes = subtract(tets_nodes, vSkinOnBodySkinNodes);
2449  Range other_tets;
2450  CHKERR moab.get_adjacencies(tets_nodes, 3, false, other_tets,
2451  moab::Interface::UNION);
2452  tets_with_four_nodes_on_skin =
2453  subtract(tets_with_four_nodes_on_skin, other_tets);
2454  Range to_remove;
2455  for (Range::iterator tit = tets_with_four_nodes_on_skin.begin();
2456  tit != tets_with_four_nodes_on_skin.end(); tit++) {
2457  int num_nodes;
2458  const EntityHandle *conn;
2459  CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
2460  double coords[12];
2461  if (th) {
2462  CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
2463  } else {
2464  CHKERR moab.get_coords(conn, num_nodes, coords);
2465  }
2466  double quality = Tools::volumeLengthQuality(coords);
2467  if (quality < 1e-2) {
2468  to_remove.insert(*tit);
2469  }
2470  }
2471 
2472  sVols = subtract(sVols, to_remove);
2473 
2474  Skinner skin(&moab);
2475  vSkin.clear();
2476  CHKERR skin.find_skin(0, sVols, false, vSkin);
2477  Range m_skin;
2478  CHKERR
2479  skin.find_skin(0, subtract(bit_ents.mSkin, to_remove), false, m_skin);
2480  vSkinWithoutBodySkin = subtract(vSkin, m_skin);
2481  vSkinOnBodySkin = intersect(vSkin, m_skin);
2482  vNodes.clear();
2483  vSkinNodes.clear();
2484  vSkinOnBodySkinNodes.clear();
2485  CHKERR moab.get_connectivity(sVols, vNodes, true);
2486  CHKERR moab.get_connectivity(vSkinOnBodySkin, vSkinOnBodySkinNodes, true);
2487  CHKERR moab.get_connectivity(vSkin, vSkinNodes, true);
2489  }
2490  };
2491 
2492  BitEnts bit_ents(m_field, mesh_bit);
2493  CHKERR bit_ents.getLevelEnts();
2494  CHKERR bit_ents.getSkin();
2495  SurfaceEnts surf_ents(m_field);
2496  CHKERR surf_ents.getVolume(bit_ents, surface);
2497  CHKERR surf_ents.getSkin(bit_ents, surface);
2498  CHKERR surf_ents.getTetsForRemesh(bit_ents);
2499 
2500  map<int, Range> types_ents;
2501  types_ents[TetGenInterface::RIDGEVERTEX].merge(
2502  surf_ents.vSkinNodesWithoutBodySkin);
2503  // FREESEGVERTEX
2504  types_ents[TetGenInterface::FREESEGVERTEX].merge(surf_ents.sSkinNodes);
2505  types_ents[TetGenInterface::FREESEGVERTEX] =
2506  subtract(types_ents[TetGenInterface::FREESEGVERTEX],
2507  types_ents[TetGenInterface::RIDGEVERTEX]);
2508  // FREEFACETVERTEX
2509  types_ents[TetGenInterface::FREEFACETVERTEX].merge(surf_ents.sNodes);
2510  types_ents[TetGenInterface::FREEFACETVERTEX] =
2511  subtract(types_ents[TetGenInterface::FREEFACETVERTEX],
2512  types_ents[TetGenInterface::RIDGEVERTEX]);
2513  types_ents[TetGenInterface::FREEFACETVERTEX] =
2514  subtract(types_ents[TetGenInterface::FREEFACETVERTEX],
2515  types_ents[TetGenInterface::FREESEGVERTEX]);
2516  // FREEVOLVERTEX
2517  types_ents[TetGenInterface::FREEVOLVERTEX].merge(surf_ents.vNodes);
2518  types_ents[TetGenInterface::FREEVOLVERTEX] =
2519  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2520  types_ents[TetGenInterface::RIDGEVERTEX]);
2521  types_ents[TetGenInterface::FREEVOLVERTEX] =
2522  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2523  types_ents[TetGenInterface::FREESEGVERTEX]);
2524  types_ents[TetGenInterface::FREEVOLVERTEX] =
2525  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2526  types_ents[TetGenInterface::FREEFACETVERTEX]);
2527 
2528  Tag th_marker;
2529  // Clean markers
2530  rval = m_field.get_moab().tag_get_handle("TETGEN_MARKER", th_marker);
2531  if(rval == MB_SUCCESS) {
2532  CHKERR m_field.get_moab().tag_delete(th_marker);
2533  rval = MB_SUCCESS;
2534  }
2535 
2536  int def_marker = 0;
2537  CHKERR m_field.get_moab().tag_get_handle(
2538  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
2539  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
2540 
2541  // Mark surface with id = 1
2542  vector<int> markers(surface.size(), 1);
2543  CHKERR m_field.get_moab().tag_set_data(th_marker, surface, &*markers.begin());
2544  // Mark all side sets
2545  int shift = 1;
2546  map<int, int> id_shift_map; // each meshset has set unique bit
2548  (*cOre.getInterface<MeshsetsManager>()), SIDESET, it)) {
2549  int ms_id = it->getMeshsetId();
2550  id_shift_map[ms_id] = 1 << shift; // shift bit
2551  ++shift;
2552  Range sideset_faces;
2553  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
2554  ms_id, SIDESET, 2, sideset_faces, true);
2555  markers.resize(sideset_faces.size());
2556  CHKERR m_field.get_moab().tag_get_data(th_marker, sideset_faces,
2557  &*markers.begin());
2558  for (unsigned int ii = 0; ii < markers.size(); ii++) {
2559  markers[ii] |= id_shift_map[ms_id]; // add bit to marker
2560  }
2561  CHKERR m_field.get_moab().tag_set_data(th_marker, sideset_faces,
2562  &*markers.begin());
2563  }
2564  Range nodes_to_remove; // none
2565  markers.resize(nodes_to_remove.size());
2566  fill(markers.begin(), markers.end(), -1);
2567  CHKERR m_field.get_moab().tag_set_data(th_marker, nodes_to_remove,
2568  &*markers.begin());
2569 
2570  // nodes
2571  if (tetGenData.size() == 1) {
2572 
2573  Range ents_to_tetgen = surf_ents.sVols;
2574  CHKERR m_field.get_moab().get_connectivity(surf_ents.sVols, ents_to_tetgen,
2575  true);
2576  for (int dd = 2; dd >= 1; dd--) {
2577  CHKERR m_field.get_moab().get_adjacencies(
2578  surf_ents.sVols, dd, false, ents_to_tetgen, moab::Interface::UNION);
2579  }
2580 
2581  // Load mesh to TetGen data structures
2582  CHKERR tetgen_iface->inData(ents_to_tetgen, in, moabTetGenMap,
2583  tetGenMoabMap, th);
2584  CHKERR tetgen_iface->setGeomData(in, moabTetGenMap, tetGenMoabMap,
2585  types_ents);
2586  std::vector<pair<Range, int> > markers;
2587  for (Range::iterator tit = surface.begin(); tit != surface.end(); tit++) {
2588  Range facet;
2589  facet.insert(*tit);
2590  markers.push_back(pair<Range, int>(facet, 2));
2591  }
2592  for (Range::iterator tit = surf_ents.vSkinWithoutBodySkin.begin();
2593  tit != surf_ents.vSkinWithoutBodySkin.end(); tit++) {
2594  Range facet;
2595  facet.insert(*tit);
2596  markers.push_back(pair<Range, int>(facet, 1));
2597  }
2598  Range other_facets;
2599  other_facets = subtract(surf_ents.vSkin, surf_ents.vSkinWithoutBodySkin);
2600  for (Range::iterator tit = other_facets.begin(); tit != other_facets.end();
2601  tit++) {
2602  Range facet;
2603  facet.insert(*tit);
2604  markers.push_back(pair<Range, int>(facet, 0));
2605  }
2606  CHKERR tetgen_iface->setFaceData(markers, in, moabTetGenMap, tetGenMoabMap);
2607  }
2608 
2609  if (debug) {
2610  if (m_field.get_comm_rank() == 0) {
2611  char tetgen_in_file_name[] = "in";
2612  in.save_nodes(tetgen_in_file_name);
2613  in.save_elements(tetgen_in_file_name);
2614  in.save_faces(tetgen_in_file_name);
2615  in.save_edges(tetgen_in_file_name);
2616  in.save_poly(tetgen_in_file_name);
2617  }
2618  }
2619 
2620  // generate new mesh
2621  {
2622  vector<string>::iterator sw = switches.begin();
2623  for (int ii = 0; sw != switches.end(); sw++, ii++) {
2624  tetgenio &_in_ = tetGenData.back();
2625  tetGenData.push_back(new tetgenio);
2626  tetgenio &_out_ = tetGenData.back();
2627  char *s = const_cast<char *>(sw->c_str());
2628  CHKERR tetgen_iface->tetRahedralize(s, _in_, _out_);
2629  }
2630  }
2631  tetgenio &out = tetGenData.back();
2632  // save elems
2633  if (debug) {
2634  char tetgen_out_file_name[] = "out";
2635  out.save_nodes(tetgen_out_file_name);
2636  out.save_elements(tetgen_out_file_name);
2637  out.save_faces(tetgen_out_file_name);
2638  out.save_edges(tetgen_out_file_name);
2639  out.save_poly(tetgen_out_file_name);
2640  }
2641 
2642  CHKERR tetgen_iface->outData(in, out, moabTetGenMap, tetGenMoabMap, bit,
2643  false, false);
2644 
2645  Range rest_of_ents = subtract(bit_ents.mTets, surf_ents.tVols);
2646  for (int dd = 2; dd >= 0; dd--) {
2647  CHKERR moab.get_adjacencies(rest_of_ents.subset_by_dimension(3), dd, false,
2648  rest_of_ents, moab::Interface::UNION);
2649  }
2650  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(rest_of_ents,
2651  bit);
2652 
2653  Range tetgen_faces;
2654  map<int, Range> face_markers_map;
2655  CHKERR tetgen_iface->getTriangleMarkers(tetGenMoabMap, out, &tetgen_faces,
2656  &face_markers_map);
2657 
2658  tetgenSurfaces = face_markers_map[1];
2659  for (map<int, Range>::iterator mit = face_markers_map.begin();
2660  mit != face_markers_map.end(); mit++) {
2662  (*cOre.getInterface<MeshsetsManager>()), SIDESET, it)) {
2663  int msId = it->getMeshsetId();
2664  if (id_shift_map[msId] & mit->first) {
2665  EntityHandle meshset = it->getMeshset();
2666  CHKERR m_field.get_moab().add_entities(
2667  meshset, mit->second.subset_by_type(MBTRI));
2668  }
2669  }
2670  }
2671 
2673 }
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:32
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Common.hpp:60
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:519
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
map< EntityHandle, unsigned long > moabTetGenMap
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:526
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
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:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439
map< unsigned long, EntityHandle > tetGenMoabMap

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

1521  {
1522  CoreInterface &m_field = cOre;
1523  moab::Interface &moab = m_field.get_moab();
1524  PrismInterface *interface;
1526  CHKERR m_field.getInterface(interface);
1527  // Remove tris on surface front
1528  {
1529  Range front_tris;
1530  EntityHandle meshset;
1531  CHKERR moab.create_meshset(MESHSET_SET, meshset);
1532  CHKERR moab.add_entities(meshset, ents);
1533  CHKERR interface->findIfTringleHasThreeNodesOnInternalSurfaceSkin(
1534  meshset, split_bit, true, front_tris);
1535  CHKERR moab.delete_entities(&meshset, 1);
1536  ents = subtract(ents, front_tris);
1537  }
1538  // Remove entities on skin
1539  Range tets;
1540  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1541  split_bit,BitRefLevel().set(),MBTET,tets
1542  );
1543  // Remove entities on skin
1544  Skinner skin(&moab);
1545  Range tets_skin;
1546  rval = skin.find_skin(0, tets, false, tets_skin);
1547  ents = subtract(ents, tets_skin);
1548 
1550 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
#define CHKERR
Inline error check.
Definition: definitions.h:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

◆ saveCutEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveCutEdges ( )

Definition at line 2731 of file CutMeshInterface.cpp.

2731  {
2732  CoreInterface &m_field = cOre;
2733  moab::Interface &moab = m_field.get_moab();
2735  CHKERR SaveData(moab)("out_vol.vtk", vOlume);
2736  CHKERR SaveData(moab)("out_surface.vtk", sUrface);
2737  CHKERR SaveData(moab)("out_cut_edges.vtk", cutEdges);
2738  CHKERR SaveData(moab)("out_cut_volumes.vtk", cutVolumes);
2739  CHKERR SaveData(moab)("out_cut_new_volumes.vtk", cutNewVolumes);
2740  CHKERR SaveData(moab)("out_cut_new_surfaces.vtk", cutNewSurfaces);
2741  CHKERR SaveData(moab)("out_cut_zero_distance_ents.vtk", zeroDistanceEnts);
2743 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
#define CHKERR
Inline error check.
Definition: definitions.h:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

◆ saveTrimEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveTrimEdges ( )

Definition at line 2745 of file CutMeshInterface.cpp.

2745  {
2746  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
2748  CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
2749  CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
2750  CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
2752 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
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:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

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

2696  {
2697  CoreInterface &m_field = cOre;
2698  moab::Interface &moab = m_field.get_moab();
2700  Range nodes;
2701  if(bit.none()) {
2702  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
2703  } else {
2704  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2705  bit, mask, MBVERTEX, nodes);
2706  }
2707  std::vector<double> coords(3 * nodes.size());
2708  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
2709  CHKERR moab.set_coords(nodes, &coords[0]);
2711 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
#define CHKERR
Inline error check.
Definition: definitions.h:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

◆ setSurface()

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

set surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 46 of file CutMeshInterface.cpp.

46  {
48  sUrface = surface;
50 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:519
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:526

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

2677  {
2678  CoreInterface &m_field = cOre;
2679  moab::Interface &moab = m_field.get_moab();
2681  Range nodes;
2682  if(bit.none()) {
2683  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
2684  } else {
2685  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2686  bit,BitRefLevel().set(),MBVERTEX,nodes
2687  );
2688  }
2689  std::vector<double> coords(3 * nodes.size());
2690  CHKERR moab.get_coords(nodes, &coords[0]);
2691  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
2693 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
#define CHKERR
Inline error check.
Definition: definitions.h:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

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

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

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

1554  {
1555  CoreInterface &m_field = cOre;
1556  moab::Interface &moab = m_field.get_moab();
1557  PrismInterface *interface;
1559  CHKERR m_field.getInterface(interface);
1560  EntityHandle meshset_volume;
1561  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
1562  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1563  split_bit, BitRefLevel().set(), MBTET, meshset_volume);
1564  EntityHandle meshset_surface;
1565  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
1566  CHKERR moab.add_entities(meshset_surface, ents);
1567  CHKERR interface->getSides(meshset_surface, split_bit, true);
1568  CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
1569  true);
1570  CHKERR moab.delete_entities(&meshset_volume, 1);
1571  CHKERR moab.delete_entities(&meshset_surface, 1);
1572  if (th) {
1573  Range prisms;
1574  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1575  bit, BitRefLevel().set(), MBPRISM, prisms);
1576  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
1577  int num_nodes;
1578  const EntityHandle *conn;
1579  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
1580  MatrixDouble data(3, 3);
1581  CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
1582  // cerr << data << endl;
1583  CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
1584  }
1585  }
1587 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:495
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Common.hpp:212
#define CHKERR
Inline error check.
Definition: definitions.h:614
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:439

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

1380  {
1381  CoreInterface &m_field = cOre;
1382  moab::Interface &moab = m_field.get_moab();
1383  MeshRefinement *refiner;
1384  const RefEntity_multiIndex *ref_ents_ptr;
1386 
1387  CHKERR m_field.getInterface(refiner);
1388  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1389  CHKERR refiner->add_verices_in_the_middel_of_edges(trimEdges, bit);
1390  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
1391  debug ? &trimEdges : NULL);
1392 
1393  trimNewVolumes.clear();
1394  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1395  bit, BitRefLevel().set(), MBTET, trimNewVolumes);
1396 
1397  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
1398  mit != edgesToTrim.end(); mit++) {
1399  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1400  boost::make_tuple(MBVERTEX, mit->first));
1401  if (vit ==
1402  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1403  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1404  "No vertex on trim edges, that make no sense");
1405  }
1406  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1407  if ((ref_ent->getBitRefLevel() & bit).any()) {
1408  EntityHandle vert = ref_ent->getRefEnt();
1409  trimNewVertices.insert(vert);
1410  verticesOnTrimEdges[vert] = mit->second;
1411  }
1412  }
1413 
1414  // Get faces which are trimmed
1415  trimNewSurfaces.clear();
1416  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1417  bit, BitRefLevel().set(), MBTRI, trimNewSurfaces);
1418 
1419  Range trim_new_surfaces_nodes;
1420  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
1421  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
1422  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
1423  Range faces_not_on_surface;
1424  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
1425  faces_not_on_surface, moab::Interface::UNION);
1426  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
1427 
1428  // Get surfaces which are not trimmed and add them to surface
1429  Range all_surfaces_on_bit_level;
1430  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1431  bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
1432  all_surfaces_on_bit_level =
1433  intersect(all_surfaces_on_bit_level, cutNewSurfaces);
1434  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
1435 
1436  Range trim_surface_edges;
1437  CHKERR moab.get_adjacencies(trimNewSurfaces, 1, false, trim_surface_edges,
1438  moab::Interface::UNION);
1439 
1440  // check of nodes are outside surface and if it are remove adjacent faces to
1441  // those nodes.
1442  Range check_verts;
1443  CHKERR moab.get_connectivity(trimNewSurfaces, check_verts, true);
1444  check_verts = subtract(check_verts, trimNewVertices);
1445  for (auto v : check_verts) {
1446 
1447  VectorDouble3 s(3);
1448  if (th) {
1449  CHKERR moab.tag_get_data(th, &v, 1, &s[0]);
1450  } else {
1451  CHKERR moab.get_coords(&v, 1, &s[0]);
1452  }
1453 
1454  VectorDouble3 p(3);
1455  EntityHandle facets_out;
1456  CHKERR treeSurfPtr->closest_to_location(&s[0], rootSetSurf, &p[0],
1457  facets_out);
1458  VectorDouble3 n(3);
1459  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
1460  VectorDouble3 delta = s - p;
1461  VectorDouble3 normal = inner_prod(delta, n) * n;
1462  if (norm_2(delta - normal) > tol * aveLength) {
1463  Range adj;
1464  CHKERR moab.get_adjacencies(&v, 1, 2, false, adj);
1465  trimNewSurfaces = subtract(trimNewSurfaces, adj);
1466  }
1467  }
1468 
1470 }
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:495
double aveLength
Average edge length.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
double tol
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:614
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:439
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Common.hpp:220

Member Data Documentation

◆ aveLength

double MoFEM::CutMeshInterface::aveLength
private

Average edge length.

Find if segment in on the plain

Parameters
s0segment first point
s1segment second point
x0point on the plain
nnormal on the plain
sintersect point
Returns
1 - intersect, 2 - segment on the plain, 0 - no intersect

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

◆ cutNewSurfaces

Range MoFEM::CutMeshInterface::cutNewSurfaces
private

Definition at line 328 of file CutMeshInterface.hpp.

◆ cutNewVertices

Range MoFEM::CutMeshInterface::cutNewVertices
private

Definition at line 331 of file CutMeshInterface.hpp.

◆ cutNewVolumes

Range MoFEM::CutMeshInterface::cutNewVolumes
private

Definition at line 327 of file CutMeshInterface.hpp.

◆ cutVolumes

Range MoFEM::CutMeshInterface::cutVolumes
private

Definition at line 326 of file CutMeshInterface.hpp.

◆ edgesToCut

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

Definition at line 350 of file CutMeshInterface.hpp.

◆ edgesToTrim

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

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

◆ mergedSurfaces

Range MoFEM::CutMeshInterface::mergedSurfaces
private

Definition at line 339 of file CutMeshInterface.hpp.

◆ mergedVolumes

Range MoFEM::CutMeshInterface::mergedVolumes
private

Definition at line 338 of file CutMeshInterface.hpp.

◆ moabTetGenMap

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

Definition at line 357 of file CutMeshInterface.hpp.

◆ nbMaxMergingCycles

int MoFEM::CutMeshInterface::nbMaxMergingCycles

Definition at line 42 of file CutMeshInterface.hpp.

◆ nbMaxTrimSearchIterations

int MoFEM::CutMeshInterface::nbMaxTrimSearchIterations

Definition at line 43 of file CutMeshInterface.hpp.

◆ rootSetSurf

EntityHandle MoFEM::CutMeshInterface::rootSetSurf
private

Definition at line 323 of file CutMeshInterface.hpp.

◆ sUrface

Range MoFEM::CutMeshInterface::sUrface
private

Definition at line 319 of file CutMeshInterface.hpp.

◆ tetGenData

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

Definition at line 359 of file CutMeshInterface.hpp.

◆ tetGenMoabMap

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

Definition at line 358 of file CutMeshInterface.hpp.

◆ tetgenSurfaces

Range MoFEM::CutMeshInterface::tetgenSurfaces
private

Definition at line 341 of file CutMeshInterface.hpp.

◆ treeSurfPtr

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

Definition at line 322 of file CutMeshInterface.hpp.

◆ trimEdges

Range MoFEM::CutMeshInterface::trimEdges
private

Definition at line 336 of file CutMeshInterface.hpp.

◆ trimNewSurfaces

Range MoFEM::CutMeshInterface::trimNewSurfaces
private

Definition at line 335 of file CutMeshInterface.hpp.

◆ trimNewVertices

Range MoFEM::CutMeshInterface::trimNewVertices
private

Definition at line 334 of file CutMeshInterface.hpp.

◆ trimNewVolumes

Range MoFEM::CutMeshInterface::trimNewVolumes
private

Definition at line 333 of file CutMeshInterface.hpp.

◆ verticesOnCutEdges

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

Definition at line 351 of file CutMeshInterface.hpp.

◆ verticesOnTrimEdges

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

Definition at line 353 of file CutMeshInterface.hpp.

◆ vOlume

Range MoFEM::CutMeshInterface::vOlume
private

Definition at line 320 of file CutMeshInterface.hpp.

◆ zeroDistanceEnts

Range MoFEM::CutMeshInterface::zeroDistanceEnts
private

Definition at line 329 of file CutMeshInterface.hpp.

◆ zeroDistanceVerts

Range MoFEM::CutMeshInterface::zeroDistanceVerts
private

Definition at line 330 of file CutMeshInterface.hpp.


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