v0.8.19
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, const bool debug=false)
 Find edges to trimEdges. More...
 
MoFEMErrorCode trimEdgesInTheMiddle (const BitRefLevel bit, Tag th=NULL, const double tol=1e-4, const bool debug=false)
 trim edges More...
 
MoFEMErrorCode moveMidNodesOnTrimmedEdges (Tag th=NULL)
 move trimmed edges mid nodes More...
 
MoFEMErrorCode removePathologicalFrontTris (const BitRefLevel split_bit, Range &ents)
 Remove pathological elements on surface internal front. More...
 
MoFEMErrorCode splitSides (const BitRefLevel split_bit, const BitRefLevel bit, const Range &ents, Tag th=NULL)
 split sides More...
 
MoFEMErrorCode mergeBadEdges (const int fraction_level, const Range &tets, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Range &merged_nodes, Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets=false, const BitRefLevel *bit_ptr=NULL, const bool debug=false)
 Merge edges. More...
 
MoFEMErrorCode mergeBadEdges (const int fraction_level, const BitRefLevel cut_bit, const BitRefLevel trim_bit, const BitRefLevel bit, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Tag th=NULL, const bool update_meshsets=false, const bool debug=false)
 Merge edges. More...
 
MoFEMErrorCode rebuildMeshWithTetGen (vector< string > &switches, const BitRefLevel &mesh_bit, const BitRefLevel &bit, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Tag th=NULL, const bool debug=false)
 
MoFEMErrorCode setTagData (Tag th, const BitRefLevel bit=BitRefLevel())
 set coords to tag More...
 
MoFEMErrorCode setCoords (Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
 set coords from tag More...
 
const Range & getVolume () const
 
const Range & getSurface () const
 
const Range & 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 47 of file CutMeshInterface.cpp.

48  : cOre(const_cast<Core &>(core)) {
49  lineSearchSteps = 10;
50  nbMaxMergingCycles = 200;
52 }

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

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

◆ clearMap()

MoFEMErrorCode MoFEM::CutMeshInterface::clearMap ( )

Definition at line 54 of file CutMeshInterface.cpp.

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

◆ copySurface()

MoFEMErrorCode MoFEM::CutMeshInterface::copySurface ( const Range &  surface,
Tag  th = NULL,
double shift = NULL,
double origin = NULL,
double transform = NULL,
const std::string  save_mesh = "" 
)

copy surface entities

Parameters
surfaceentities which going to be added
Returns
error code
Examples:
mesh_cut.cpp.

Definition at line 66 of file CutMeshInterface.cpp.

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

◆ cutAndTrim()

MoFEMErrorCode MoFEM::CutMeshInterface::cutAndTrim ( 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 164 of file CutMeshInterface.cpp.

168  {
169  CoreInterface &m_field = cOre;
170  moab::Interface &moab = m_field.get_moab();
172 
173  // cut mesh
174  CHKERR findEdgesToCut(fixed_edges, corner_nodes, tol_cut, QUIET, debug);
175  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_cut_close,
176  QUIET, debug);
177  CHKERR cutEdgesInMiddle(bit_level1, debug);
178  if (fixed_edges) {
179  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
180  *fixed_edges);
181  }
182  if (corner_nodes) {
183  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
184  *corner_nodes);
185  }
186  if (update_meshsets) {
187  CHKERR UpdateMeshsets()(cOre, bit_level1);
188  }
190 
191  auto get_min_quality = [&m_field](const BitRefLevel bit, Tag th) {
192  Range tets_level; // test at level
193  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
194  bit, BitRefLevel().set(), MBTET, tets_level);
195  double min_q = 1;
196  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
197  return min_q;
198  };
199 
200  PetscPrintf(PETSC_COMM_WORLD, "Min quality cut %6.4g\n",
201  get_min_quality(bit_level1, th));
202 
203  if(debug) {
205  }
206 
207  // trim mesh
208  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim, debug);
209  CHKERR trimEdgesInTheMiddle(bit_level2, th, tol_trim_close, debug);
210  if (fixed_edges) {
211  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
212  *fixed_edges);
213  }
214  if (corner_nodes) {
215  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
216  *corner_nodes);
217  }
218  if (update_meshsets) {
219  CHKERR UpdateMeshsets()(cOre, bit_level2);
220  }
222 
223  PetscPrintf(PETSC_COMM_WORLD, "Min quality trim %3.2g\n",
224  get_min_quality(bit_level2, th));
225 
226  if(debug) {
228  Range bit2_edges;
229  CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
230  bit_level2, BitRefLevel().set(), MBEDGE, bit2_edges);
231  CHKERR SaveData(moab)("trim_fixed_edges.vtk",
232  intersect(*fixed_edges, bit2_edges));
233  }
234 
236 }
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:475
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:594
MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, const bool debug=false)
Find edges to trimEdges.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
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:475
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:594
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
map< EntityHandle, TreeData > edgesToCut

◆ cutTrimAndMerge()

MoFEMErrorCode MoFEM::CutMeshInterface::cutTrimAndMerge ( 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 238 of file CutMeshInterface.cpp.

243  {
244  CoreInterface &m_field = cOre;
246  if(debug) {
247  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
248  "ents_not_in_database.vtk", "VTK", "");
249  }
250  CHKERR cutAndTrim(bit_level1, bit_level2, th, tol_cut, tol_cut_close,
251  tol_trim, tol_trim_close, &fixed_edges, &corner_nodes,
252  update_meshsets, debug);
253  if(debug) {
254  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
255  "cut_trim_ents_not_in_database.vtk", "VTK", "");
256  }
257 
258  CHKERR mergeBadEdges(fraction_level, bit_level2, bit_level1, bit_level3,
259  getNewTrimSurfaces(), fixed_edges, corner_nodes, th,
260  update_meshsets, debug);
261  // CHKERR removePathologicalFrontTris(bit_level3,
262  // const_cast<Range &>(getMergedSurfaces()));
263 
264  if(debug) {
265  CHKERR cOre.getInterface<BitRefManager>()->writeBitLevelByType(
266  bit_level3, BitRefLevel().set(), MBTET, "out_tets_merged.vtk", "VTK",
267  "");
268  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
269  "cut_trim_merge_ents_not_in_database.vtk", "VTK", "");
270  }
271 
272  auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
273  Range tets_level; // test at level
274  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
275  bit, BitRefLevel().set(), MBTET, tets_level);
276  double min_q = 1;
277  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
278  if (min_q < 0 && debug) {
279  CHKERR m_field.getInterface<Tools>()->writeTetsWithQuality(
280  "negative_tets.vtk", "VTK", "", tets_level, th);
281  }
282  return min_q;
283  };
284 
285  PetscPrintf(PETSC_COMM_WORLD, "Min quality node merge %6.4g\n",
286  get_min_quality(bit_level1, th));
287 
288  CHKERR
289  cOre.getInterface<BitRefManager>()->updateRange(fixed_edges, fixed_edges);
290  CHKERR cOre.getInterface<BitRefManager>()->updateRange(corner_nodes,
291  corner_nodes);
292 
294 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
MoFEMErrorCode mergeBadEdges(const int fraction_level, const Range &tets, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Range &merged_nodes, Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets=false, const BitRefLevel *bit_ptr=NULL, const bool debug=false)
Merge edges.
MoFEMErrorCode 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:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ findEdgesToCut()

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

find edges to cut

Parameters
verbverbosity level
Returns
error code

Definition at line 327 of file CutMeshInterface.cpp.

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

Find edges to trimEdges.

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

Parameters
verbverbosity level
Returns
error code

Definition at line 1052 of file CutMeshInterface.cpp.

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

◆ getRayForEdge()

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

Definition at line 1457 of file CutMeshInterface.cpp.

1460  {
1461  const CoreInterface &m_field = cOre;
1462  const moab::Interface &moab = m_field.get_moab();
1464  int num_nodes;
1465  const EntityHandle *conn;
1466  CHKERR moab.get_connectivity(ent, conn, num_nodes, true);
1467  double coords[6];
1468  CHKERR moab.get_coords(conn, num_nodes, coords);
1469  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
1470  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
1471  noalias(ray_point) = s0;
1472  noalias(unit_ray_dir) = s1 - s0;
1473  ray_length = norm_2(unit_ray_dir);
1474  unit_ray_dir /= ray_length;
1476 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Common.hpp:234
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

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

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

2239  {
2240  CoreInterface &m_field = cOre;
2242  Range tets_level;
2243  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2244  trim_bit, BitRefLevel().set(), MBTET, tets_level);
2245  CHKERRG(ierr);
2246 
2247  Range edges_to_merge;
2248  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByParentType(
2249  trim_bit, trim_bit | cut_bit, MBEDGE, edges_to_merge);
2250  CHKERRG(ierr);
2251  edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
2252 
2253  // get all entities not in database
2254  Range all_ents_not_in_database_before;
2255  ierr = cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2256  all_ents_not_in_database_before);
2257  CHKERRG(ierr);
2258 
2259  Range out_new_tets, out_new_surf;
2260  ierr = mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
2261  corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
2262  th, update_meshsets, &bit, debug);
2263  CHKERRG(ierr);
2264 
2265  // get all entities not in database after merge
2266  Range all_ents_not_in_database_after;
2267  ierr = cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2268  all_ents_not_in_database_after);
2269  CHKERRG(ierr);
2270  // delete hanging entities
2271  all_ents_not_in_database_after =
2272  subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
2273  Range meshsets;
2274  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
2275  false);
2276  for (Range::iterator mit = meshsets.begin(); mit != meshsets.end(); mit++) {
2277  CHKERR m_field.get_moab().remove_entities(*mit,
2278  all_ents_not_in_database_after);
2279  }
2280  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
2281 
2282  mergedVolumes.swap(out_new_tets);
2283  mergedSurfaces.swap(out_new_surf);
2285 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
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:542
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506
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:594

◆ mergeSurface()

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

merge surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 131 of file CutMeshInterface.cpp.

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

◆ mergeVolumes()

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

merge volume entities

Parameters
volumeentities which going to be added
Returns
error code

Definition at line 137 of file CutMeshInterface.cpp.

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

◆ moveMidNodesOnCutEdges()

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

projecting of mid edge nodes on new mesh on surface

Returns
error code

Definition at line 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:499
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
VectorBoundedArray< double, 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 (auto &v : verticesOnTrimEdges) {
1041  double dist = v.second.dIst;
1042  VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1043  if (th) {
1044  CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1045  } else {
1046  CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1047  }
1048  }
1050 }
map< EntityHandle, TreeData > verticesOnTrimEdges
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
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 641 of file CutMeshInterface.cpp.

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

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

◆ rebuildMeshWithTetGen()

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

Definition at line 2289 of file CutMeshInterface.cpp.

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

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

1506  {
1507  CoreInterface &m_field = cOre;
1508  moab::Interface &moab = m_field.get_moab();
1509  PrismInterface *interface;
1511  CHKERR m_field.getInterface(interface);
1512  // Remove tris on surface front
1513  {
1514  Range front_tris;
1515  EntityHandle meshset;
1516  CHKERR moab.create_meshset(MESHSET_SET, meshset);
1517  CHKERR moab.add_entities(meshset, ents);
1518  CHKERR interface->findIfTringleHasThreeNodesOnInternalSurfaceSkin(
1519  meshset, split_bit, true, front_tris);
1520  CHKERR moab.delete_entities(&meshset, 1);
1521  ents = subtract(ents, front_tris);
1522  }
1523  // Remove entities on skin
1524  Range tets;
1525  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1526  split_bit,BitRefLevel().set(),MBTET,tets
1527  );
1528  // Remove entities on skin
1529  Skinner skin(&moab);
1530  Range tets_skin;
1531  rval = skin.find_skin(0, tets, false, tets_skin);
1532  ents = subtract(ents, tets_skin);
1533 
1535 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
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:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ saveCutEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveCutEdges ( )

Definition at line 2684 of file CutMeshInterface.cpp.

2684  {
2685  CoreInterface &m_field = cOre;
2686  moab::Interface &moab = m_field.get_moab();
2688  CHKERR SaveData(moab)("out_vol.vtk", vOlume);
2689  CHKERR SaveData(moab)("out_surface.vtk", sUrface);
2690  CHKERR SaveData(moab)("out_cut_edges.vtk", cutEdges);
2691  CHKERR SaveData(moab)("out_cut_volumes.vtk", cutVolumes);
2692  CHKERR SaveData(moab)("out_cut_new_volumes.vtk", cutNewVolumes);
2693  CHKERR SaveData(moab)("out_cut_new_surfaces.vtk", cutNewSurfaces);
2694  CHKERR SaveData(moab)("out_cut_zero_distance_ents.vtk", zeroDistanceEnts);
2696 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ saveTrimEdges()

MoFEMErrorCode MoFEM::CutMeshInterface::saveTrimEdges ( )

Definition at line 2698 of file CutMeshInterface.cpp.

2698  {
2699  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
2701  CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
2702  CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
2703  CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
2705 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ setCoords()

MoFEMErrorCode MoFEM::CutMeshInterface::setCoords ( Tag  th,
const BitRefLevel  bit = BitRefLevel(),
const BitRefLevel  mask = BitRefLevel().set() 
)

set coords from tag

Parameters
thtag handle
Returns
error code
Examples:
mesh_cut.cpp.

Definition at line 2666 of file CutMeshInterface.cpp.

2667  {
2668  CoreInterface &m_field = cOre;
2669  moab::Interface &moab = m_field.get_moab();
2671  Range nodes;
2672  if(bit.none()) {
2673  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
2674  } else {
2675  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2676  bit, mask, MBVERTEX, nodes);
2677  }
2678  std::vector<double> coords(3 * nodes.size());
2679  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
2680  CHKERR moab.set_coords(nodes, &coords[0]);
2682 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
#define CHKERR
Inline error check.
Definition: definitions.h:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ setSurface()

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

set surface entities

Parameters
surfaceentities which going to be added
Returns
error code

Definition at line 60 of file CutMeshInterface.cpp.

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

◆ setTagData()

MoFEMErrorCode MoFEM::CutMeshInterface::setTagData ( Tag  th,
const BitRefLevel  bit = BitRefLevel() 
)

set coords to tag

Parameters
thtag handle
Returns
error code
Examples:
mesh_cut.cpp.

Definition at line 2648 of file CutMeshInterface.cpp.

2648  {
2649  CoreInterface &m_field = cOre;
2650  moab::Interface &moab = m_field.get_moab();
2652  Range nodes;
2653  if(bit.none()) {
2654  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
2655  } else {
2656  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2657  bit,BitRefLevel().set(),MBVERTEX,nodes
2658  );
2659  }
2660  std::vector<double> coords(3 * nodes.size());
2661  CHKERR moab.get_coords(nodes, &coords[0]);
2662  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
2664 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
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:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ setVolume()

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

set volume entities

Parameters
volumeentities which going to be added
Returns
error code
Examples:
mesh_cut.cpp.

Definition at line 125 of file CutMeshInterface.cpp.

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

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

1539  {
1540  CoreInterface &m_field = cOre;
1541  moab::Interface &moab = m_field.get_moab();
1542  PrismInterface *interface;
1544  CHKERR m_field.getInterface(interface);
1545  EntityHandle meshset_volume;
1546  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
1547  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1548  split_bit, BitRefLevel().set(), MBTET, meshset_volume);
1549  EntityHandle meshset_surface;
1550  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
1551  CHKERR moab.add_entities(meshset_surface, ents);
1552  CHKERR interface->getSides(meshset_surface, split_bit, true);
1553  CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
1554  true);
1555  CHKERR moab.delete_entities(&meshset_volume, 1);
1556  CHKERR moab.delete_entities(&meshset_surface, 1);
1557  if (th) {
1558  Range prisms;
1559  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1560  bit, BitRefLevel().set(), MBPRISM, prisms);
1561  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
1562  int num_nodes;
1563  const EntityHandle *conn;
1564  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
1565  MatrixDouble data(3, 3);
1566  CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
1567  // cerr << data << endl;
1568  CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
1569  }
1570  }
1572 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
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:594
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405

◆ trimEdgesInTheMiddle()

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

trim edges

Parameters
bitbit level of the trimmed mesh
Returns
error code

Definition at line 1363 of file CutMeshInterface.cpp.

1365  {
1366  CoreInterface &m_field = cOre;
1367  moab::Interface &moab = m_field.get_moab();
1368  MeshRefinement *refiner;
1369  const RefEntity_multiIndex *ref_ents_ptr;
1371 
1372  CHKERR m_field.getInterface(refiner);
1373  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1374  CHKERR refiner->add_verices_in_the_middel_of_edges(trimEdges, bit);
1375  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
1376  debug ? &trimEdges : NULL);
1377 
1378  trimNewVolumes.clear();
1379  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1380  bit, BitRefLevel().set(), MBTET, trimNewVolumes);
1381 
1382  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
1383  mit != edgesToTrim.end(); mit++) {
1384  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1385  boost::make_tuple(MBVERTEX, mit->first));
1386  if (vit ==
1387  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1388  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1389  "No vertex on trim edges, that make no sense");
1390  }
1391  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1392  if ((ref_ent->getBitRefLevel() & bit).any()) {
1393  EntityHandle vert = ref_ent->getRefEnt();
1394  trimNewVertices.insert(vert);
1395  verticesOnTrimEdges[vert] = mit->second;
1396  }
1397  }
1398 
1399  // Get faces which are trimmed
1400  trimNewSurfaces.clear();
1401  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1402  bit, BitRefLevel().set(), MBTRI, trimNewSurfaces);
1403 
1404  Range trim_new_surfaces_nodes;
1405  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
1406  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
1407  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
1408  Range faces_not_on_surface;
1409  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
1410  faces_not_on_surface, moab::Interface::UNION);
1411  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
1412 
1413  // Get surfaces which are not trimmed and add them to surface
1414  Range all_surfaces_on_bit_level;
1415  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1416  bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
1417  all_surfaces_on_bit_level =
1418  intersect(all_surfaces_on_bit_level, cutNewSurfaces);
1419  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
1420 
1421  Range trim_surface_edges;
1422  CHKERR moab.get_adjacencies(trimNewSurfaces, 1, false, trim_surface_edges,
1423  moab::Interface::UNION);
1424 
1425  // check of nodes are outside surface and if it are remove adjacent faces to
1426  // those nodes.
1427  Range check_verts;
1428  CHKERR moab.get_connectivity(trimNewSurfaces, check_verts, true);
1429  check_verts = subtract(check_verts, trimNewVertices);
1430  for (auto v : check_verts) {
1431 
1432  VectorDouble3 s(3);
1433  if (th) {
1434  CHKERR moab.tag_get_data(th, &v, 1, &s[0]);
1435  } else {
1436  CHKERR moab.get_coords(&v, 1, &s[0]);
1437  }
1438 
1439  VectorDouble3 p(3);
1440  EntityHandle facets_out;
1441  CHKERR treeSurfPtr->closest_to_location(&s[0], rootSetSurf, &p[0],
1442  facets_out);
1443  VectorDouble3 n(3);
1444  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
1445  VectorDouble3 delta = s - p;
1446  VectorDouble3 normal = inner_prod(delta, n) * n;
1447  if (norm_2(delta - normal) > tol * aveLength) {
1448  Range adj;
1449  CHKERR moab.get_adjacencies(&v, 1, 2, false, adj);
1450  trimNewSurfaces = subtract(trimNewSurfaces, adj);
1451  }
1452  }
1453 
1455 }
map< EntityHandle, TreeData > edgesToTrim
map< EntityHandle, TreeData > verticesOnTrimEdges
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
double aveLength
Average edge length.
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:594
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
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: