1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // File : SMESH_FillHole.cxx
23 // Created : Tue Sep 26 15:11:17 2017
24 // Author : Edward AGAPOV (eap)
27 #include "SMESH_MeshAlgos.hxx"
29 #include "SMESH_Comment.hxx"
31 #include "SMESH_TypeDefs.hxx"
32 #include "SMDS_Mesh.hxx"
34 #include <Utils_SALOME_Exception.hxx>
36 #include <boost/intrusive/circular_list_algorithms.hpp>
37 #include <boost/container/flat_map.hpp>
39 #include <Bnd_B3d.hxx>
43 bool isSmallAngle( double cos2 )
45 // cosine of min angle at which adjacent faces are considered overlapping
46 const double theMinCos2 = 0.996 * 0.996; // ~5 degrees
47 return ( cos2 > theMinCos2 );
51 typedef std::multimap< double, BEdge* > TAngleMap;
52 typedef std::map< const SMDS_MeshElement*, int > TFaceIndMap;
54 //--------------------------------------------------------------------------------
56 * \brief Edge of a free border
60 const SMDS_MeshNode* myNode1;
61 const SMDS_MeshNode* myNode2;
62 const SMDS_MeshElement* myFace; // face adjacent to the border
65 gp_XYZ myDir; // myNode1 -> myNode2
66 double myDirCoef; // 1. or -1, to make myDir oriented as myNodes in myFace
67 double myLength; // between nodes
68 double myAngleWithPrev; // between myDir and -myPrev->myDir
69 TAngleMap::iterator myAngleMapPos;
70 double myOverlapAngle; // angle delta due to overlapping
71 const SMDS_MeshNode* myNode1Shift; // nodes created to avoid overlapping of faces
72 const SMDS_MeshNode* myNode2Shift;
74 BEdge* myPrev; // neighbors in the border
77 BEdge(): myNode1Shift(0), myNode2Shift(0) {}
78 void Init( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2,
79 const SMDS_MeshElement* f=0,
80 const SMDS_MeshNode* nf1=0, const SMDS_MeshNode* nf2=0 );
81 void ComputeAngle( bool reverseAngle = false );
82 void ShiftOverlapped( const SMDS_MeshNode* oppNode,
83 const TFaceIndMap& capFaceWithBordInd,
85 std::vector<const SMDS_MeshElement*>& newFaces);
86 void MakeShiftfFaces( SMDS_Mesh& mesh,
87 std::vector<const SMDS_MeshElement*>& newFaces,
88 const bool isReverse );
89 gp_XYZ GetInFaceDir() const { return myFaceNorm ^ myDir * myDirCoef; }
90 void InsertSelf(TAngleMap& edgesByAngle, bool isReverseFaces, bool reBind, bool useOverlap )
92 if ( reBind ) edgesByAngle.erase( myAngleMapPos );
93 double key = (( isReverseFaces ? 2 * M_PI - myAngleWithPrev : myAngleWithPrev )
94 + myOverlapAngle * useOverlap );
95 myAngleMapPos = edgesByAngle.insert( std::make_pair( key, this ));
98 // traits used by boost::intrusive::circular_list_algorithms
100 typedef BEdge * node_ptr;
101 typedef const BEdge * const_node_ptr;
102 static node_ptr get_next(const_node_ptr n) { return n->myNext; }
103 static void set_next(node_ptr n, node_ptr next) { n->myNext = next; }
104 static node_ptr get_previous(const_node_ptr n) { return n->myPrev; }
105 static void set_previous(node_ptr n, node_ptr prev){ n->myPrev = prev; }
108 //================================================================================
110 * \brief Initialize a border edge data
112 //================================================================================
114 void BEdge::Init( const SMDS_MeshNode* n1,
115 const SMDS_MeshNode* n2,
116 const SMDS_MeshElement* newFace, // new cap face
117 const SMDS_MeshNode* nf1,
118 const SMDS_MeshNode* nf2 )
122 myDir = SMESH_NodeXYZ( n2 ) - SMESH_NodeXYZ( n1 );
123 myLength = myDir.Modulus();
124 if ( myLength > std::numeric_limits<double>::min() )
130 TIDSortedElemSet elemSet, avoidSet;
132 myFace = SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet, &ind1, &ind2 );
134 throw SALOME_Exception( SMESH_Comment("No face sharing nodes #")
135 << myNode1->GetID() << " and #" << myNode2->GetID());
136 avoidSet.insert( myFace );
137 if ( SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet ))
138 throw SALOME_Exception( SMESH_Comment("No free border between nodes #")
139 << myNode1->GetID() << " and #" << myNode2->GetID());
141 myDirCoef = SMESH_MeshAlgos::IsRightOrder( myFace, myNode1, myNode2 ) ? 1. : -1.;
144 if (! SMESH_MeshAlgos::FaceNormal( myFace, myFaceNorm, /*normalized=*/false ))
146 SMDS_ElemIteratorPtr fIt = myNode1->GetInverseElementIterator( SMDSAbs_Face );
147 while ( fIt->more() )
148 if ( SMESH_MeshAlgos::FaceNormal( fIt->next(), myFaceNorm, /*normalized=*/false ))
155 myDirCoef = SMESH_MeshAlgos::IsRightOrder( newFace, nf1, nf2 ) ? 1. : -1.;
156 if ( myPrev->myNode2 == n1 )
157 myNode1Shift = myPrev->myNode2Shift;
158 if ( myNext->myNode1 == n2 )
159 myNode2Shift = myNext->myNode1Shift;
161 else if ( myDirCoef * myPrev->myDirCoef < 0 ) // different orientation of faces
169 //================================================================================
171 * \brief Compute myAngleWithPrev
173 //================================================================================
175 void BEdge::ComputeAngle( bool theReverseAngle )
177 myAngleWithPrev = ACos( myDir.Dot( myPrev->myDir.Reversed() ));
180 gp_XYZ inNewFaceDir = myDir - myPrev->myDir;
181 double dot1 = myDir.Dot( myPrev->myFaceNorm );
182 double dot2 = myPrev->myDir.Dot( myFaceNorm );
183 bool isOverlap1 = ( inNewFaceDir * myPrev->GetInFaceDir() > 0 );
184 bool isOverlap2 = ( inNewFaceDir * GetInFaceDir() > 0 );
185 if ( !myPrev->myFace )
186 isObtuse = isOverlap1;
188 isObtuse = isOverlap2;
191 isObtuse = ( dot1 > 0 || dot2 < 0 ); // suppose face normals point outside the border
192 if ( theReverseAngle )
193 isObtuse = !isObtuse;
197 myAngleWithPrev = 2 * M_PI - myAngleWithPrev;
202 // isSmallAngle( 1 - myDir.CrossSquareMagnitude( myPrev->myDir )); // edges co-directed
207 // check if myFace and a triangle built on this and prev edges overlap
210 double cos2 = dot1 * dot1 / myPrev->myFaceNorm.SquareModulus();
211 myOverlapAngle += 0.5 * M_PI * ( 1 - cos2 );
215 double cos2 = dot2 * dot2 / myFaceNorm.SquareModulus();
216 myOverlapAngle += 0.5 * M_PI * ( 1 - cos2 );
221 //================================================================================
223 * \brief Check if myFace is overlapped by a triangle formed by myNode's and a
224 * given node. If so, create shifted nodes to avoid overlapping
226 //================================================================================
228 void BEdge::ShiftOverlapped( const SMDS_MeshNode* theOppNode,
229 const TFaceIndMap& theCapFaceWithBordInd,
231 std::vector<const SMDS_MeshElement*>& theNewFaces )
233 if ( myNode1Shift && myNode2Shift )
236 gp_XYZ inNewFaceDir = SMESH_NodeXYZ( theOppNode ) - SMESH_NodeXYZ( myNode1 );
237 double dot = inNewFaceDir.Dot( myFaceNorm );
238 double cos2 = dot * dot / myFaceNorm.SquareModulus() / inNewFaceDir.SquareModulus();
239 bool isOverlap = ( isSmallAngle( 1 - cos2 ) && GetInFaceDir() * inNewFaceDir > 0 );
243 gp_XYZ shift = myFaceNorm / myLength / 4;
248 gp_XYZ p = SMESH_NodeXYZ( myNode1 ) + shift;
249 myNode1Shift = theMesh.AddNode( p.X(), p.Y(), p.Z() );
250 myPrev->myNode2Shift = myNode1Shift;
254 gp_XYZ p = SMESH_NodeXYZ( myNode2 ) + shift;
255 myNode2Shift = theMesh.AddNode( p.X(), p.Y(), p.Z() );
256 myNext->myNode1Shift = myNode2Shift;
259 // MakeShiftfFaces() for already created cap faces
260 for ( int is2nd = 0; is2nd < 2; ++is2nd )
262 const SMDS_MeshNode* ns = is2nd ? myNode2Shift : myNode1Shift;
263 const SMDS_MeshNode* n = is2nd ? myNode2 : myNode1;
266 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
267 while ( fIt->more() )
269 const SMDS_MeshElement* f = fIt->next();
270 if ( !f->isMarked() ) continue;
272 TFaceIndMap::const_iterator f2i = theCapFaceWithBordInd.find( f );
273 if ( f2i == theCapFaceWithBordInd.end() )
275 const SMDS_MeshNode* nf1 = f->GetNode( f2i->second );
276 const SMDS_MeshNode* nf2 = f->GetNode(( f2i->second+1 ) % f->NbNodes() );
277 if ( nf1 == n || nf2 == n )
280 tmpE.myPrev = tmpE.myNext = this;
281 tmpE.Init( nf1, nf2, f, nf1, nf2 );
282 if ( !tmpE.myNode1Shift && !tmpE.myNode2Shift )
283 tmpE.Init( nf2, nf1, f, nf2, nf1 );
285 tmpE.MakeShiftfFaces( theMesh, theNewFaces, tmpE.myDirCoef < 0 );
287 std::vector< const SMDS_MeshNode* > nodes( f->begin_nodes(), f->end_nodes() );
288 nodes[ f->GetNodeIndex( n ) ] = ns;
289 theMesh.ChangeElementNodes( f, &nodes[0], nodes.size() );
295 //================================================================================
297 * \brief Create a triangle
299 //================================================================================
301 const SMDS_MeshElement* MakeTria( SMDS_Mesh& mesh,
302 const SMDS_MeshNode* n1,
303 const SMDS_MeshNode* n2,
304 const SMDS_MeshNode* n3,
305 const bool isReverse )
308 return mesh.AddFace( n1, n3, n2 );
309 return mesh.AddFace( n1, n2, n3 );
312 //================================================================================
314 * \brief Create a quadrangle
316 //================================================================================
318 // const SMDS_MeshElement* MakeQuad( SMDS_Mesh& mesh,
319 // const SMDS_MeshNode* n1,
320 // const SMDS_MeshNode* n2,
321 // const SMDS_MeshNode* n3,
322 // const SMDS_MeshNode* n4,
323 // const bool isReverse )
326 // return mesh.AddFace( n4, n3, n2, n1 );
327 // return mesh.AddFace( n1, n2, n3, n4 );
330 //================================================================================
332 * \brief Create faces on myNode* and myNode*Shift
334 //================================================================================
336 void BEdge::MakeShiftfFaces(SMDS_Mesh& mesh,
337 std::vector<const SMDS_MeshElement*>& newFaces,
338 const bool isReverse )
342 if ( myNode1Shift && myNode2Shift )
344 newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode2Shift, isReverse ));
345 newFaces.push_back( MakeTria( mesh, myNode1, myNode2Shift, myNode1Shift, isReverse ));
347 else if ( myNode1Shift )
349 newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode1Shift, isReverse ));
351 else if ( myNode2Shift )
353 newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode2Shift, isReverse ));
359 //================================================================================
361 * \brief Fill with 2D elements a hole defined by a TFreeBorder
363 //================================================================================
365 void SMESH_MeshAlgos::FillHole(const SMESH_MeshAlgos::TFreeBorder & theFreeBorder,
367 std::vector<const SMDS_MeshElement*>& theNewFaces)
369 if ( theFreeBorder.size() < 4 || // at least 3 nodes
370 theFreeBorder[0] != theFreeBorder.back() ) // the hole must be closed
373 // prepare data of the border
375 ObjectPool< BEdge > edgeAllocator;
376 boost::intrusive::circular_list_algorithms< BEdge > circularList;
378 BEdge* edge0 = edgeAllocator.getNew();
379 BEdge* edgePrev = edge0;
380 circularList.init_header( edge0 );
381 edge0->Init( theFreeBorder[0], theFreeBorder[1], 0 );
383 box.Add( SMESH_NodeXYZ( edge0->myNode1 ));
384 for ( size_t i = 2; i < theFreeBorder.size(); ++i )
386 edge = edgeAllocator.getNew();
387 circularList.link_after( edgePrev, edge );
388 edge->Init( theFreeBorder[i-1], theFreeBorder[i] );
389 edge->ComputeAngle();
391 box.Add( SMESH_NodeXYZ( edge->myNode1 ));
393 edge0->ComputeAngle();
395 // check if face normals point outside the border
397 gp_XYZ hSize = 0.5 * ( box.CornerMax() - box.CornerMin() );
398 const double hDelta = 1e-6 * hSize.Modulus();
399 hSize -= gp_XYZ( hDelta, hDelta, hDelta );
400 if ( hSize.X() < 0 ) hSize.SetX(hDelta);
401 if ( hSize.Y() < 0 ) hSize.SetY(hDelta);
402 if ( hSize.Z() < 0 ) hSize.SetZ(hDelta);
403 box.SetHSize( hSize ); // decrease the box by hDelta
405 size_t nbEdges = theFreeBorder.size() - 1;
407 int nbRev = 0, nbFrw = 0;
408 double angTol = M_PI - ( nbEdges - 2 ) * M_PI / nbEdges, sumDirCoeff = 0;
409 for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
411 if ( box.IsOut( SMESH_NodeXYZ( edge->myNode1 )) &&
412 edge->myOverlapAngle < 0.1 * M_PI )
414 nbRev += edge->myAngleWithPrev > M_PI + angTol;
415 nbFrw += edge->myAngleWithPrev < M_PI - angTol;
417 sumDirCoeff += edge->myDirCoef;
419 // unmark all adjacent faces, new faces will be marked
420 SMDS_ElemIteratorPtr fIt = edge->myNode1->GetInverseElementIterator( SMDSAbs_Face );
421 while ( fIt->more() )
422 fIt->next()->setIsMarked( false );
424 bool isReverseAngle = ( nbRev > nbFrw ); // true == face normals point inside the border
425 //std::cout << "nbRev="<< nbRev << ", nbFrw="<< nbFrw<<std::endl;
427 // sort border edges by myAngleWithPrev
429 TAngleMap edgesByAngle;
430 bool useOverlap = true; // to add BEdge.myOverlapAngle when filling edgesByAngle
432 for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
433 edge->InsertSelf( edgesByAngle, isReverseAngle, /*reBind=*/false, useOverlap );
435 // create triangles to fill the hole
437 //compare order of nodes in the edges with their order in faces
438 bool isReverse = sumDirCoeff > 0.5 * nbEdges;
440 // faces filling the hole (cap faces) and indices of border edges in them
441 TFaceIndMap capFaceWithBordInd;
443 theNewFaces.reserve( nbEdges - 2 );
444 std::vector< const SMDS_MeshNode* > nodes(3);
445 while ( edgesByAngle.size() > 2 )
447 TAngleMap::iterator a2e = edgesByAngle.begin();
448 if ( useOverlap && a2e->first > M_PI - angTol ) // all new triangles need shift
453 nbEdges = edgesByAngle.size();
454 edgesByAngle.clear();
455 for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
456 edge->InsertSelf( edgesByAngle, isReverseAngle, /*reBind=*/false, useOverlap );
457 a2e = edgesByAngle.begin();
460 edgePrev = edge->myPrev;
462 // create shift nodes and faces
463 edgePrev->ShiftOverlapped( edge->myNode2, capFaceWithBordInd, theMesh, theNewFaces );
464 edge->ShiftOverlapped( edgePrev->myNode1, capFaceWithBordInd, theMesh, theNewFaces );
465 edge ->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
466 edgePrev->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
470 nodes[0] = edgePrev->myNode1Shift ? edgePrev->myNode1Shift : edgePrev->myNode1;
471 nodes[1] = edgePrev->myNode2Shift ? edgePrev->myNode2Shift : edgePrev->myNode2;
472 nodes[2] = edge->myNode2Shift ? edge->myNode2Shift : edge->myNode2;
473 theNewFaces.push_back( MakeTria( theMesh, nodes[0], nodes[1], nodes[2], isReverse ));
474 // std::cout << nodes[1]->GetID() << " " << nodes[0]->GetID() << " " << nodes[2]->GetID()
475 // << " " << edge->myAngleWithPrev << std::endl;
477 // remember a border edge within the new cap face
478 theNewFaces.back()->setIsMarked( true );
479 if ( edgePrev->myFace )
480 capFaceWithBordInd.insert( std::make_pair( theNewFaces.back(), isReverse ? 2 : 0 ));
482 capFaceWithBordInd.insert( std::make_pair( theNewFaces.back(), 1 ));
484 // remove edgePrev from the list and update <edge>
485 edgesByAngle.erase( edgePrev->myAngleMapPos );
486 circularList.unlink( edgePrev ); // remove edgePrev from the border
488 edge->Init( edgePrev->myNode1, edge->myNode2, theNewFaces.back(), nodes[0], nodes[2] );
489 edge->ComputeAngle( isReverseAngle );
490 edge->InsertSelf( edgesByAngle, /*isReverse=*/false, /*reBind=*/true, useOverlap );
491 edge->myNext->ComputeAngle( isReverseAngle );
492 edge->myNext->InsertSelf( edgesByAngle, /*isReverse=*/false, /*reBind=*/true, useOverlap );
493 // std::cout << "A " << edge->myNode1->GetID() << " " << edge->myAngleWithPrev
494 // << " " << edge->myNext->myNode1->GetID() << " " << edge->myNext->myAngleWithPrev
497 edge = edgesByAngle.begin()->second;
498 edge-> MakeShiftfFaces( theMesh, theNewFaces, isReverse );
499 edge->myNext->MakeShiftfFaces( theMesh, theNewFaces, isReverse );