]> SALOME platform Git repositories - modules/smesh.git/blob - src/SMESHUtils/SMESH_FillHole.cxx
Salome HOME
Merge branch 'master' into gni/adaptation
[modules/smesh.git] / src / SMESHUtils / SMESH_FillHole.cxx
1 // Copyright (C) 2007-2020  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
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.
10 //
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.
15 //
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
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22 // File      : SMESH_FillHole.cxx
23 // Created   : Tue Sep 26 15:11:17 2017
24 // Author    : Edward AGAPOV (eap)
25 //
26
27 #include "SMESH_MeshAlgos.hxx"
28
29 #include "SMESH_Comment.hxx"
30
31 #include "ObjectPool.hxx"
32 #include "SMDS_Mesh.hxx"
33 #include "SMESH_TypeDefs.hxx"
34
35 #include <Utils_SALOME_Exception.hxx>
36
37 #include <boost/intrusive/circular_list_algorithms.hpp>
38 //#include <boost/container/flat_map.hpp>
39
40 #include <Bnd_B3d.hxx>
41
42 namespace
43 {
44   bool isSmallAngle( double cos2 )
45   {
46     // cosine of min angle at which adjacent faces are considered overlapping
47     const double theMinCos2 = 0.996 * 0.996; // ~5 degrees
48     return ( cos2 > theMinCos2 );
49   }
50
51   struct BEdge;
52   typedef std::multimap< double, BEdge* >          TAngleMap;
53   typedef std::map< const SMDS_MeshElement*, int > TFaceIndMap;
54
55   //--------------------------------------------------------------------------------
56   /*!
57    * \brief Edge of a free border
58    */
59   struct BEdge
60   {
61     const SMDS_MeshNode*    myNode1;
62     const SMDS_MeshNode*    myNode2;
63     const SMDS_MeshElement* myFace;   // face adjacent to the border
64
65     gp_XYZ                  myFaceNorm;
66     gp_XYZ                  myDir;     // myNode1 -> myNode2
67     double                  myDirCoef; // 1. or -1, to make myDir oriented as myNodes in myFace
68     double                  myLength;  // between nodes
69     double                  myAngleWithPrev; // between myDir and -myPrev->myDir
70     double                  myMinMaxRatio; // of a possible triangle sides
71     TAngleMap::iterator     myAngleMapPos;
72     double                  myOverlapAngle;  // angle delta due to overlapping
73     const SMDS_MeshNode*    myNode1Shift;    // nodes created to avoid overlapping of faces
74     const SMDS_MeshNode*    myNode2Shift;
75
76     BEdge*                  myPrev; // neighbors in the border
77     BEdge*                  myNext;
78
79     BEdge(): myNode1Shift(0), myNode2Shift(0) {}
80     void   Init( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2,
81                  const SMDS_MeshElement* f=0,
82                  const SMDS_MeshNode* nf1=0, const SMDS_MeshNode* nf2=0 );
83     void   ComputeAngle( bool reverseAngle = false );
84     void   ShiftOverlapped( const SMDS_MeshNode*                  oppNode,
85                             const TFaceIndMap&                    capFaceWithBordInd,
86                             SMDS_Mesh&                            mesh,
87                             std::vector<const SMDS_MeshElement*>& newFaces);
88     void   MakeShiftfFaces( SMDS_Mesh&                            mesh,
89                             std::vector<const SMDS_MeshElement*>& newFaces,
90                             const bool                            isReverse );
91     gp_XYZ GetInFaceDir() const { return myFaceNorm ^ myDir * myDirCoef; }
92     double ShapeFactor()  const { return 0.5 * ( 1. - myMinMaxRatio ); }
93     void   InsertSelf(TAngleMap& edgesByAngle, bool isReverseFaces, bool reBind, bool useOverlap )
94     {
95       if ( reBind ) edgesByAngle.erase( myAngleMapPos );
96       double key = (( isReverseFaces ? 2 * M_PI - myAngleWithPrev : myAngleWithPrev )
97                     + myOverlapAngle * useOverlap
98                     + ShapeFactor() );
99       myAngleMapPos = edgesByAngle.insert( std::make_pair( key, this ));
100     }
101
102     // traits used by boost::intrusive::circular_list_algorithms
103     typedef BEdge         node;
104     typedef BEdge *       node_ptr;
105     typedef const BEdge * const_node_ptr;
106     static node_ptr get_next(const_node_ptr n)             {  return n->myNext;  }
107     static void     set_next(node_ptr n, node_ptr next)    {  n->myNext = next;  }
108     static node_ptr get_previous(const_node_ptr n)         {  return n->myPrev;  }
109     static void     set_previous(node_ptr n, node_ptr prev){  n->myPrev = prev;  }
110   };
111
112   //================================================================================
113   /*!
114    * \brief Initialize a border edge data
115    */
116   //================================================================================
117
118   void BEdge::Init( const SMDS_MeshNode*    n1,
119                     const SMDS_MeshNode*    n2,
120                     const SMDS_MeshElement* newFace, // new cap face
121                     const SMDS_MeshNode*    nf1,
122                     const SMDS_MeshNode*    nf2 )
123   {
124     myNode1  = n1;
125     myNode2  = n2;
126     myDir    = SMESH_NodeXYZ( n2 ) - SMESH_NodeXYZ( n1 );
127     myLength = myDir.Modulus();
128     if ( myLength > std::numeric_limits<double>::min() )
129       myDir /= myLength;
130
131     myFace = newFace;
132     if ( !myFace )
133     {
134       TIDSortedElemSet elemSet, avoidSet;
135       int ind1, ind2;
136       myFace = SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet, &ind1, &ind2 );
137       if ( !myFace )
138         throw SALOME_Exception( SMESH_Comment("No face sharing nodes #")
139                                 << myNode1->GetID() << " and #" << myNode2->GetID());
140       avoidSet.insert( myFace );
141       if ( SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet ))
142         throw SALOME_Exception( SMESH_Comment("No free border between nodes #")
143                                 << myNode1->GetID() << " and #" << myNode2->GetID());
144
145       myDirCoef = SMESH_MeshAlgos::IsRightOrder( myFace, myNode1, myNode2 ) ? 1. : -1.;
146     }
147
148     if (! SMESH_MeshAlgos::FaceNormal( myFace, myFaceNorm, /*normalized=*/false ))
149     {
150       SMDS_ElemIteratorPtr fIt = myNode1->GetInverseElementIterator( SMDSAbs_Face );
151       while ( fIt->more() )
152         if ( SMESH_MeshAlgos::FaceNormal( fIt->next(), myFaceNorm, /*normalized=*/false ))
153           break;
154     }
155
156     if ( newFace )
157     {
158       myFace    = 0;
159       myDirCoef = SMESH_MeshAlgos::IsRightOrder( newFace, nf1, nf2 ) ? 1. : -1.;
160       if ( myPrev->myNode2 == n1 )
161         myNode1Shift = myPrev->myNode2Shift;
162       if ( myNext->myNode1 == n2 )
163         myNode2Shift = myNext->myNode1Shift;
164     }
165     else if ( myDirCoef * myPrev->myDirCoef < 0 ) // different orientation of faces
166     {
167       myFaceNorm *= -1;
168       myDirCoef  *= -1;
169     }
170   }
171
172   //================================================================================
173   /*!
174    * \brief Compute myAngleWithPrev
175    */
176   //================================================================================
177
178   void BEdge::ComputeAngle( bool theReverseAngle )
179   {
180     double dot = myDir.Dot( myPrev->myDir.Reversed() );
181     if      ( dot >=  1 ) myAngleWithPrev = 0;
182     else if ( dot <= -1 ) myAngleWithPrev = M_PI;
183     else                  myAngleWithPrev = acos( dot );
184
185     bool isObtuse;
186     gp_XYZ inFaceDirNew = myDir - myPrev->myDir;
187     gp_XYZ   inFaceDir1 = myPrev->GetInFaceDir();
188     gp_XYZ   inFaceDir2 = this->GetInFaceDir();
189     double         dot1 = inFaceDirNew * inFaceDir1;
190     double         dot2 = inFaceDirNew * inFaceDir2;
191     bool     isOverlap1 = ( dot1 > 0 );
192     bool     isOverlap2 = ( dot2 > 0 );
193     if ( !myPrev->myFace )
194       isObtuse = isOverlap1;
195     else if  ( !myFace )
196       isObtuse = isOverlap2;
197     else
198     {
199       double dt1 = myDir.Dot( myPrev->myFaceNorm );
200       double dt2 = myPrev->myDir.Dot( myFaceNorm );
201       isObtuse = ( dt1 > 0 || dt2 < 0 ); // suppose face normals point outside the border
202       if ( theReverseAngle )
203         isObtuse = !isObtuse;
204     }
205     if ( isObtuse )
206     {
207       myAngleWithPrev = 2 * M_PI - myAngleWithPrev;
208     }
209
210     // if ( ! isObtuse )
211     //   isObtuse =
212     //     isSmallAngle( 1 - myDir.CrossSquareMagnitude( myPrev->myDir )); // edges co-directed
213
214     myOverlapAngle = 0;
215     //if ( !isObtuse )
216     {
217       // check if myFace and a triangle built on this and prev edges overlap
218       if ( isOverlap1 )
219       {
220         double cos2 = dot1 * dot1 / inFaceDirNew.SquareModulus() / inFaceDir1.SquareModulus();
221         myOverlapAngle += 1. * M_PI * cos2;
222       }
223       if ( isOverlap2 )
224       {
225         double cos2 = dot2 * dot2 / inFaceDirNew.SquareModulus() / inFaceDir2.SquareModulus();
226         myOverlapAngle += 1. * M_PI * cos2;
227       }
228     }
229
230     {
231       double len3 = SMESH_NodeXYZ( myPrev->myNode1 ).Distance( myNode2 );
232       double minLen = Min( myLength, Min( myPrev->myLength, len3 ));
233       double maxLen = Max( myLength, Max( myPrev->myLength, len3 ));
234       myMinMaxRatio = minLen / maxLen;
235     }
236   }
237
238   //================================================================================
239   /*!
240    * \brief Check if myFace is overlapped by a triangle formed by myNode's and a
241    *        given node. If so, create shifted nodes to avoid overlapping
242    */
243   //================================================================================
244
245   void BEdge::ShiftOverlapped( const SMDS_MeshNode*                  theOppNode,
246                                const TFaceIndMap&                    theCapFaceWithBordInd,
247                                SMDS_Mesh&                            theMesh,
248                                std::vector<const SMDS_MeshElement*>& theNewFaces )
249   {
250     if ( myNode1Shift && myNode2Shift )
251       return;
252
253     gp_XYZ inNewFaceDir = SMESH_NodeXYZ( theOppNode ) - SMESH_NodeXYZ( myNode1 );
254     double          dot = inNewFaceDir.Dot( myFaceNorm );
255     double         cos2 = dot * dot / myFaceNorm.SquareModulus() / inNewFaceDir.SquareModulus();
256     bool      isOverlap = ( isSmallAngle( 1 - cos2 ) && GetInFaceDir() * inNewFaceDir > 0 );
257
258     if ( isOverlap )
259     {
260       gp_XYZ shift = myFaceNorm / myLength / 4;
261       if ( myFace )
262         shift.Reverse();
263       if ( !myNode1Shift )
264       {
265         gp_XYZ p = SMESH_NodeXYZ( myNode1 ) + shift;
266         myNode1Shift = theMesh.AddNode( p.X(), p.Y(), p.Z() );
267         myPrev->myNode2Shift = myNode1Shift;
268       }
269       if ( !myNode2Shift )
270       {
271         gp_XYZ p = SMESH_NodeXYZ( myNode2 ) + shift;
272         myNode2Shift = theMesh.AddNode( p.X(), p.Y(), p.Z() );
273         myNext->myNode1Shift = myNode2Shift;
274       }
275
276       // MakeShiftfFaces() for already created cap faces
277       for ( int is2nd = 0; is2nd < 2; ++is2nd )
278       {
279         const SMDS_MeshNode* ns = is2nd ? myNode2Shift : myNode1Shift;
280         const SMDS_MeshNode* n  = is2nd ? myNode2 : myNode1;
281         if ( !ns ) continue;
282
283         SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
284         while ( fIt->more() )
285         {
286           const SMDS_MeshElement* f = fIt->next();
287           if ( !f->isMarked() ) continue;
288
289           TFaceIndMap::const_iterator f2i = theCapFaceWithBordInd.find( f );
290           if ( f2i == theCapFaceWithBordInd.end() )
291             continue;
292           const SMDS_MeshNode* nf1 = f->GetNode(  f2i->second );
293           const SMDS_MeshNode* nf2 = f->GetNode(( f2i->second+1 ) % f->NbNodes() );
294           if ( nf1 == n || nf2 == n )
295           {
296             BEdge tmpE;
297             tmpE.myPrev = tmpE.myNext = this;
298             tmpE.Init( nf1, nf2, f, nf1, nf2 );
299             if ( !tmpE.myNode1Shift && !tmpE.myNode2Shift )
300               tmpE.Init( nf2, nf1, f, nf2, nf1 );
301             tmpE.myFace = f;
302             tmpE.MakeShiftfFaces( theMesh, theNewFaces, tmpE.myDirCoef < 0 );
303           }
304           std::vector< const SMDS_MeshNode* > nodes( f->begin_nodes(), f->end_nodes() );
305           nodes[ f->GetNodeIndex( n ) ] = ns;
306           theMesh.ChangeElementNodes( f, &nodes[0], nodes.size() );
307         }
308       }
309     }
310   }
311
312   //================================================================================
313   /*!
314    * \brief Create a triangle
315    */
316   //================================================================================
317
318   const SMDS_MeshElement* MakeTria( SMDS_Mesh&           mesh,
319                                     const SMDS_MeshNode* n1,
320                                     const SMDS_MeshNode* n2,
321                                     const SMDS_MeshNode* n3,
322                                     const bool           isReverse )
323   {
324     if ( isReverse )
325       return mesh.AddFace( n1, n3, n2 );
326     return mesh.AddFace( n1, n2, n3 );
327   }
328
329   //================================================================================
330   /*!
331    * \brief Create a quadrangle
332    */
333   //================================================================================
334
335   // const SMDS_MeshElement* MakeQuad( SMDS_Mesh&           mesh,
336   //                                   const SMDS_MeshNode* n1,
337   //                                   const SMDS_MeshNode* n2,
338   //                                   const SMDS_MeshNode* n3,
339   //                                   const SMDS_MeshNode* n4,
340   //                                   const bool           isReverse )
341   // {
342   //   if ( isReverse )
343   //     return mesh.AddFace( n4, n3, n2, n1 );
344   //   return mesh.AddFace( n1, n2, n3, n4 );
345   // }
346
347   //================================================================================
348   /*!
349    * \brief Create faces on myNode* and myNode*Shift
350    */
351   //================================================================================
352
353   void BEdge::MakeShiftfFaces(SMDS_Mesh&                            mesh,
354                               std::vector<const SMDS_MeshElement*>& newFaces,
355                               const bool                            isReverse )
356   {
357     if ( !myFace )
358       return;
359     if ( myNode1Shift && myNode2Shift )
360     {
361       newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode2Shift, isReverse ));
362       newFaces.push_back( MakeTria( mesh, myNode1, myNode2Shift, myNode1Shift, isReverse ));
363     }
364     else if ( myNode1Shift )
365     {
366       newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode1Shift, isReverse ));
367     }
368     else if ( myNode2Shift )
369     {
370       newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode2Shift, isReverse ));
371     }
372   }
373
374 } // namespace
375
376 //================================================================================
377 /*!
378  * \brief Fill with 2D elements a hole defined by a TFreeBorder
379  */
380 //================================================================================
381
382 void SMESH_MeshAlgos::FillHole(const SMESH_MeshAlgos::TFreeBorder &  theFreeBorder,
383                                SMDS_Mesh&                            theMesh,
384                                std::vector<const SMDS_MeshElement*>& theNewFaces)
385 {
386   if ( theFreeBorder.size() < 4 ||                // at least 3 nodes
387        theFreeBorder[0] != theFreeBorder.back() ) // the hole must be closed
388     return;
389
390   // prepare data of the border
391
392   ObjectPool< BEdge > edgeAllocator;
393   boost::intrusive::circular_list_algorithms< BEdge > circularList;
394   BEdge* edge;
395   BEdge* edge0 = edgeAllocator.getNew();
396   BEdge* edgePrev = edge0;
397   circularList.init_header( edge0 );
398   edge0->Init( theFreeBorder[0], theFreeBorder[1], 0 );
399   Bnd_B3d box;
400   box.Add( SMESH_NodeXYZ( edge0->myNode1 ));
401   for ( size_t i = 2; i < theFreeBorder.size(); ++i )
402   {
403     edge = edgeAllocator.getNew();
404     circularList.link_after( edgePrev, edge );
405     edge->Init( theFreeBorder[i-1], theFreeBorder[i] );
406     edge->ComputeAngle();
407     edgePrev = edge;
408     box.Add( SMESH_NodeXYZ( edge->myNode1 ));
409   }
410   edge0->ComputeAngle();
411
412   // check if face normals point outside the border
413
414   gp_XYZ hSize = 0.5 * ( box.CornerMax() - box.CornerMin() );
415   const double hDelta = 1e-6 * hSize.Modulus();
416   hSize -= gp_XYZ( hDelta, hDelta, hDelta );
417   if ( hSize.X() < 0 ) hSize.SetX(hDelta);
418   if ( hSize.Y() < 0 ) hSize.SetY(hDelta);
419   if ( hSize.Z() < 0 ) hSize.SetZ(hDelta);
420   box.SetHSize( hSize ); // decrease the box by hDelta
421
422   size_t nbEdges = theFreeBorder.size() - 1;
423   edge = edge0;
424   int nbRev = 0, nbFrw = 0;
425   double angTol = M_PI - ( nbEdges - 2 ) * M_PI / nbEdges, sumDirCoeff = 0;
426   for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
427   {
428     if ( box.IsOut( SMESH_NodeXYZ( edge->myNode1 )) &&
429          edge->myOverlapAngle < 0.1 * M_PI )
430     {
431       nbRev += edge->myAngleWithPrev > M_PI + angTol;
432       nbFrw += edge->myAngleWithPrev < M_PI - angTol;
433     }
434     sumDirCoeff += edge->myDirCoef;
435
436     // unmark all adjacent faces, new faces will be marked
437     SMDS_ElemIteratorPtr fIt = edge->myNode1->GetInverseElementIterator( SMDSAbs_Face );
438     while ( fIt->more() )
439       fIt->next()->setIsMarked( false );
440   }
441   bool isReverseAngle = ( nbRev > nbFrw ); // true == face normals point inside the border
442   //std::cout << "nbRev="<< nbRev << ", nbFrw="<< nbFrw<<std::endl;
443
444   // sort border edges by myAngleWithPrev
445
446   TAngleMap edgesByAngle;
447   bool useOverlap = true; // to add BEdge.myOverlapAngle when filling edgesByAngle
448   edge = edge0;
449   for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
450     edge->InsertSelf( edgesByAngle, isReverseAngle, /*reBind=*/false, useOverlap );
451
452   // create triangles to fill the hole
453
454   //compare order of nodes in the edges with their order in faces
455   bool isReverse = sumDirCoeff > 0.5 * nbEdges;
456
457   // faces filling the hole (cap faces) and indices of border edges in them
458   TFaceIndMap capFaceWithBordInd;
459
460   theNewFaces.reserve( nbEdges - 2 );
461   std::vector< const SMDS_MeshNode* > nodes(3);
462   while ( edgesByAngle.size() > 2 )
463   {
464     TAngleMap::iterator a2e = edgesByAngle.begin();
465     edge = a2e->second;
466     if ( useOverlap &&
467          a2e->first - edge->ShapeFactor() > M_PI - angTol ) // all new triangles need shift
468     {
469       // re-sort the edges w/o overlap consideration
470       useOverlap = false;
471       nbEdges = edgesByAngle.size();
472       edgesByAngle.clear();
473       for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
474         edge->InsertSelf( edgesByAngle, isReverseAngle, /*reBind=*/false, useOverlap );
475       a2e = edgesByAngle.begin();
476     }
477     edge     = a2e->second;
478     edgePrev = edge->myPrev;
479
480     // create shift nodes and faces
481     edgePrev->ShiftOverlapped( edge->myNode2, capFaceWithBordInd, theMesh, theNewFaces );
482     edge->ShiftOverlapped( edgePrev->myNode1, capFaceWithBordInd, theMesh, theNewFaces );
483     edge    ->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
484     edgePrev->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
485
486     // make a cap face
487     //nodes.resize( 3 );
488     nodes[0] = edgePrev->myNode1Shift ? edgePrev->myNode1Shift : edgePrev->myNode1;
489     nodes[1] = edgePrev->myNode2Shift ? edgePrev->myNode2Shift : edgePrev->myNode2;
490     nodes[2] = edge->myNode2Shift     ? edge->myNode2Shift     : edge->myNode2;
491     theNewFaces.push_back( MakeTria( theMesh, nodes[0], nodes[1], nodes[2], isReverse ));
492     // std::cout << nodes[1]->GetID() << " "  << nodes[0]->GetID() << " "  << nodes[2]->GetID()
493     //           << " " << edge->myAngleWithPrev << std::endl;
494
495     // remember a border edge within the new cap face
496     theNewFaces.back()->setIsMarked( true );
497     if ( edgePrev->myFace )
498       capFaceWithBordInd.insert( std::make_pair( theNewFaces.back(), isReverse ? 2 : 0 ));
499     if ( edge->myFace )
500       capFaceWithBordInd.insert( std::make_pair( theNewFaces.back(), 1 ));
501
502     // remove edgePrev from the list and update <edge>
503     edgesByAngle.erase( edgePrev->myAngleMapPos );
504     circularList.unlink( edgePrev ); // remove edgePrev from the border
505
506     edge->Init( edgePrev->myNode1, edge->myNode2, theNewFaces.back(), nodes[0], nodes[2] );
507     edge->ComputeAngle( isReverseAngle );
508     edge->InsertSelf( edgesByAngle, /*isReverse=*/false, /*reBind=*/true, useOverlap );
509     edge->myNext->ComputeAngle( isReverseAngle );
510     edge->myNext->InsertSelf( edgesByAngle, /*isReverse=*/false, /*reBind=*/true, useOverlap );
511     // std::cout << "A " << edge->myNode1->GetID() << " " << edge->myAngleWithPrev
512     //           << " " << edge->myNext->myNode1->GetID() << " " << edge->myNext->myAngleWithPrev
513     //           << std::endl;
514   }
515   edge = edgesByAngle.begin()->second;
516   edge->        MakeShiftfFaces( theMesh, theNewFaces, isReverse );
517   edge->myNext->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
518 }