Salome HOME
bos #20256: [CEA 18523] Porting SMESH to int 64 bits
[modules/smesh.git] / src / SMESHUtils / SMESH_Offset.cxx
index 1b1f5204e179a0fbe31b1ba89975622e51419e4f..6d30cf6e2b31267dfe45a21561ef3e17578c7550 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2019  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2021  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -48,7 +48,7 @@ namespace
   //--------------------------------------------------------------------------------
   /*!
    * \brief Intersected face side storing a node created at this intersection
-   *        and a intersected face
+   *        and an intersected face
    */
   struct CutLink
   {
@@ -105,7 +105,7 @@ namespace
     int                     myIndex; // positive -> side index, negative -> State
     const SMDS_MeshElement* myFace;
 
-    enum State { _INTERNAL = -1, _COPLANAR = -2 };
+    enum State { _INTERNAL = -1, _COPLANAR = -2, _PENDING = -3 };
 
     void Set( const SMDS_MeshNode*    Node1,
               const SMDS_MeshNode*    Node2,
@@ -277,7 +277,7 @@ namespace
 
     static Standard_Integer HashCode(const CutFace& f, const Standard_Integer upper)
     {
-      return ::HashCode( f.myInitFace->GetID(), upper );
+      return ::HashCode( FromSmIdType<int>(f.myInitFace->GetID()), upper );
     }
     static Standard_Boolean IsEqual(const CutFace& f1, const CutFace& f2 )
     {
@@ -552,7 +552,7 @@ namespace
 
   bool getTranslatedPosition( const SMDS_MeshNode*         theNewNode,
                               const double                 theOffset,
-                              const double                 theTol,
+                              const double                 /*theTol*/,
                               const double                 theSign,
                               const std::vector< gp_XYZ >& theFaceNormals,
                               SMDS_Mesh&                   theSrcMesh,
@@ -712,6 +712,53 @@ namespace
     return useOneNormal;
   }
 
+  //================================================================================
+  /*!
+   * \brief Remove small faces
+   */
+  //================================================================================
+
+  void removeSmallFaces( SMDS_Mesh*                        theMesh,
+                         SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces,
+                         const double                      theTol2 )
+  {
+    std::vector< SMESH_NodeXYZ > points(3);
+    std::vector< const SMDS_MeshNode* > nodes(3);
+    for ( SMDS_ElemIteratorPtr faceIt = theMesh->elementsIterator(); faceIt->more(); )
+    {
+      const SMDS_MeshElement* face = faceIt->next();
+      points.assign( face->begin_nodes(), face->end_nodes() );
+
+      SMESH_NodeXYZ* prevN = & points.back();
+      for ( size_t i = 0; i < points.size(); ++i )
+      {
+        double dist2 = ( *prevN - points[ i ]).SquareModulus();
+        if ( dist2 < theTol2 )
+        {
+          const SMDS_MeshNode* nToRemove =
+            (*prevN)->GetID() > points[ i ]->GetID() ? prevN->Node() : points[ i ].Node();
+          const SMDS_MeshNode* nToKeep =
+            nToRemove == points[ i ].Node() ? prevN->Node() : points[ i ].Node();
+          for ( SMDS_ElemIteratorPtr fIt = nToRemove->GetInverseElementIterator(); fIt->more(); )
+          {
+            const SMDS_MeshElement* f = fIt->next();
+            if ( f == face )
+              continue;
+            nodes.assign( f->begin_nodes(), f->end_nodes() );
+            nodes[ f->GetNodeIndex( nToRemove )] = nToKeep;
+            theMesh->ChangeElementNodes( f, &nodes[0], nodes.size() );
+          }
+          theNew2OldFaces[ face->GetID() ].first = 0;
+          theMesh->RemoveFreeElement( face );
+          break;
+        }
+        prevN = & points[ i ];
+      }
+      continue;
+    }
+    return;
+  }
+
 } // namespace
 
 namespace SMESH_MeshAlgos
@@ -1239,6 +1286,8 @@ namespace SMESH_MeshAlgos
         if ( !link2.IntNode() ) link2.myIntNode = link1.myIntNode;
 
         cf1.AddPoint( link1, link2, myTol );
+        if ( l1n1 ) link1.Set( l1n1, l1n2, face2 );
+        if ( l2n1 ) link2.Set( l2n1, l2n2, face2 );
         cf2.AddPoint( link1, link2, myTol );
       }
       else
@@ -1524,7 +1573,8 @@ namespace SMESH_MeshAlgos
     SMESH_MeshAlgos::GetBarycentricCoords( p2D( p ),
                                            p2D( nodes[0] ), p2D( nodes[1] ), p2D( nodes[2] ),
                                            bc1, bc2 );
-    return ( 0. < bc1 && 0. < bc2 && bc1 + bc2 < 1. );
+    //return ( 0. < bc1 && 0. < bc2 && bc1 + bc2 < 1. );
+    return ( myTol < bc1 && myTol < bc2 && bc1 + bc2 + myTol < 1. );
   }
 
   //================================================================================
@@ -1676,7 +1726,11 @@ namespace SMESH_MeshAlgos
 
     size_t limit = cf.myLinks.size() * cf.myLinks.size() * 2;
 
-    for ( size_t i1 = 3; i1 < cf.myLinks.size(); ++i1 )
+    size_t i1 = 3;
+    while ( cf.myLinks[i1-1].IsInternal() && i1 > 0 )
+      --i1;
+
+    for ( ; i1 < cf.myLinks.size(); ++i1 )
     {
       if ( !cf.myLinks[i1].IsInternal() )
         continue;
@@ -1879,8 +1933,8 @@ namespace SMESH_MeshAlgos
       for ( ; cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
       {
         const CutFace& cf = *cutFacesIt;
-        int index = cf.myInitFace->GetID(); // index in theNew2OldFaces
-        if ((int) theNew2OldFaces.size() <= index )
+        smIdType index = cf.myInitFace->GetID(); // index in theNew2OldFaces
+        if ((smIdType) theNew2OldFaces.size() <= index )
           theNew2OldFaces.resize( index + 1 );
         theNew2OldFaces[ index ] = std::make_pair( cf.myInitFace, index );
       }
@@ -1950,7 +2004,7 @@ namespace SMESH_MeshAlgos
       // erase loops that are cut off by face intersections
       cf.CutOffLoops( loopSet, theSign, myNormals, cutOffLinks, cutOffCoplanarLinks );
 
-      int index = cf.myInitFace->GetID(); // index in theNew2OldFaces
+      smIdType index = cf.myInitFace->GetID(); // index in theNew2OldFaces
 
       const SMDS_MeshElement* tria;
       for ( size_t iL = 0; iL < loopSet.myNbLoops; ++iL )
@@ -1968,6 +2022,8 @@ namespace SMESH_MeshAlgos
 #ifdef _DEBUG_
             std::cerr << "BAD tria" << std::endl;
             cf.Dump();
+#else
+            if ( i < 0 ) cf.Dump(); // avoid "CutFace::Dump() unused in release mode"
 #endif
             continue;
           }
@@ -1998,6 +2054,28 @@ namespace SMESH_MeshAlgos
         myMesh->RemoveFreeElement( f );
     }
 
+    // remove faces that are merged off
+    for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
+    {
+      const CutFace& cf = *cutFacesIt;
+      if ( !cf.myLinks.empty() || cf.myInitFace->IsNull() )
+        continue;
+
+      nodes.assign( cf.myInitFace->begin_nodes(), cf.myInitFace->end_nodes() );
+      for ( size_t i = 0; i < nodes.size(); ++i )
+      {
+        const SMDS_MeshNode* n = nodes[ i ];
+        while ( myRemove2KeepNodes.IsBound( n ))
+          n = myRemove2KeepNodes( n );
+        if ( n != nodes[ i ] && cf.myInitFace->GetNodeIndex( n ) >= 0 )
+        {
+          theNew2OldFaces[ cf.myInitFace->GetID() ].first = 0;
+          myMesh->RemoveFreeElement( cf.myInitFace );
+          break;
+        }
+      }
+    }
+
     // remove faces connected to cut off parts of cf.myInitFace
 
     nodes.resize(2);
@@ -2010,13 +2088,13 @@ namespace SMESH_MeshAlgos
       if ( nodes[0] != nodes[1] &&
            myMesh->GetElementsByNodes( nodes, faces ))
       {
-        if ( cutOffLinks[i].myFace &&
-             cutOffLinks[i].myIndex != EdgePart::_COPLANAR &&
-             faces.size() == 2 )
+        if ( // cutOffLinks[i].myFace &&
+            cutOffLinks[i].myIndex != EdgePart::_COPLANAR &&
+            faces.size() != 1 )
           continue;
         for ( size_t iF = 0; iF < faces.size(); ++iF )
         {
-          int index = faces[iF]->GetID();
+          smIdType index = faces[iF]->GetID();
           // if ( //faces[iF]->isMarked()         ||  // kept part of cutFace
           //      !theNew2OldFaces[ index ].first ) // already removed
           //   continue;
@@ -2045,13 +2123,22 @@ namespace SMESH_MeshAlgos
     for ( size_t i = 0; i < touchedFaces.size(); ++i )
     {
       const CutFace& cf = *touchedFaces[i];
+      if ( cf.myInitFace->IsNull() )
+        continue;
 
-      int index = cf.myInitFace->GetID(); // index in theNew2OldFaces
+      smIdType index = cf.myInitFace->GetID(); // index in theNew2OldFaces
       if ( !theNew2OldFaces[ index ].first )
         continue; // already cut off
 
+      cf.InitLinks();
       if ( !cf.ReplaceNodes( myRemove2KeepNodes ))
-        continue; // just keep as is
+      {
+        if ( cf.myLinks.size() == 3 &&
+             cf.myInitFace->GetNodeIndex( cf.myLinks[0].myNode1 ) >= 0 &&
+             cf.myInitFace->GetNodeIndex( cf.myLinks[1].myNode1 ) >= 0 &&
+             cf.myInitFace->GetNodeIndex( cf.myLinks[2].myNode1 ) >= 0 )
+          continue; // just keep as is
+      }
 
       if ( cf.myLinks.size() == 3 )
       {
@@ -2391,8 +2478,9 @@ namespace
           // if ( !myLinks[i].IsInternal() )
           //   myLinks[ i ].myFace = cutterFace;
           // else
-          myLinks[ i   ].ReplaceCoplanar( newEdge );
-          myLinks[ i+1 ].ReplaceCoplanar( newEdge );
+          myLinks[ i ].ReplaceCoplanar( newEdge );
+          if ( myLinks[i].IsInternal() && i+1 < myLinks.size() )
+            myLinks[ i+1 ].ReplaceCoplanar( newEdge );
           return;
         }
         i += myLinks[i].IsInternal();
@@ -2411,7 +2499,7 @@ namespace
    */
   //================================================================================
 
-  void CutFace::AddPoint( const CutLink& p1, const CutLink& p2, double tol ) const
+  void CutFace::AddPoint( const CutLink& p1, const CutLink& p2, double /*tol*/ ) const
   {
     if ( myInitFace->GetNodeIndex( p1.IntNode() ) >= 0 ||
          myInitFace->GetNodeIndex( p2.IntNode() ) >= 0 )
@@ -2670,7 +2758,7 @@ namespace
       // add links connecting internal loops with the boundary ones
 
       // find a pair of closest nodes
-      const SMDS_MeshNode *closestNode1, *closestNode2;
+      const SMDS_MeshNode *closestNode1 = 0, *closestNode2 = 0;
       double minDist = 1e100;
       for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE )
       {
@@ -2736,8 +2824,16 @@ namespace
     {
       theLoops.AddNewLoop();
       theLoops.AddEdge( myLinks[0] );
-      theLoops.AddEdge( myLinks[1] );
-      theLoops.AddEdge( myLinks[2] );
+      if ( myLinks[0].myNode2 == myLinks[1].myNode1 )
+      {
+        theLoops.AddEdge( myLinks[1] );
+        theLoops.AddEdge( myLinks[2] );
+      }
+      else
+      {
+        theLoops.AddEdge( myLinks[2] );
+        theLoops.AddEdge( myLinks[1] );
+      }
       return;
     }
 
@@ -2824,9 +2920,11 @@ namespace
                              const double                 theSign,
                              const std::vector< gp_XYZ >& theNormals,
                              std::vector< EdgePart >&     theCutOffLinks,
-                             TLinkMap&                    theCutOffCoplanarLinks) const
+                             TLinkMap&                    /*theCutOffCoplanarLinks*/) const
   {
     EdgePart sideEdge;
+    boost::container::flat_set< const SMDS_MeshElement* > checkedCoplanar;
+
     for ( size_t i = 0; i < myLinks.size(); ++i )
     {
       if ( !myLinks[i].myFace )
@@ -2874,6 +2972,21 @@ namespace
             loop = theLoops.GetLoopOf( twin );
           toErase = ( loop && !loop->myLinks.empty() );
         }
+
+        if ( toErase ) // do not erase if cutFace is connected to a co-planar cutFace
+        {
+          checkedCoplanar.clear();
+          for ( size_t iE = 0; iE < myLinks.size() &&  toErase; ++iE )
+          {
+            if ( !myLinks[iE].myFace || myLinks[iE].myIndex != EdgePart::_COPLANAR )
+              continue;
+            bool isAdded = checkedCoplanar.insert( myLinks[iE].myFace ).second;
+            if ( !isAdded )
+              continue;
+            toErase = ( SMESH_MeshAlgos::NbCommonNodes( myLinks[i ].myFace,
+                                                        myLinks[iE].myFace ) < 1 );
+          }
+        }
       }
 
       if ( toErase )
@@ -2885,8 +2998,8 @@ namespace
           {
             if ( !loop->myLinks[ iE ]->myFace              &&
                  !loop->myLinks[ iE ]->IsInternal()     )//   &&
-                 // !loop->myLinks[ iE ]->myNode1->isMarked() && // cut nodes are marked
-                 // !loop->myLinks[ iE ]->myNode2->isMarked() )
+              // !loop->myLinks[ iE ]->myNode1->isMarked() && // cut nodes are marked
+              // !loop->myLinks[ iE ]->myNode2->isMarked() )
             {
               int i = loop->myLinks[ iE ]->myIndex;
               sideEdge.Set( myInitFace->GetNode    ( i   ),
@@ -2962,7 +3075,9 @@ namespace
     if ( myIndex + e.myIndex == _COPLANAR + _INTERNAL )
     {
       //check if the faces are connected
-      int nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e.myFace, myFace ).size();
+      int nbCommonNodes = 0;
+      if ( e.myFace && myFace )
+        nbCommonNodes = SMESH_MeshAlgos::NbCommonNodes( e.myFace, myFace );
       bool toReplace = (( myIndex == _INTERNAL && nbCommonNodes > 1 ) ||
                         ( myIndex == _COPLANAR && nbCommonNodes < 2 ));
       if ( toReplace )
@@ -3069,7 +3184,7 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
     {
       p.Set( nodes[i] );
       double dist = ( pPrev - p ).SquareModulus();
-      if ( dist > std::numeric_limits<double>::min() )
+      if ( dist < minNodeDist && dist > std::numeric_limits<double>::min() )
         minNodeDist = dist;
       pPrev = p;
     }
@@ -3111,7 +3226,7 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
     for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
     {
       const SMDS_MeshElement* newFace = fIt->next();
-      const int             faceIndex = newFace->GetID();
+      const smIdType        faceIndex = newFace->GetID();
       const gp_XYZ&           oldNorm = normals[ faceIndex ];
       const gp_XYZ             newXYZ = oldXYZ + oldNorm * theOffset;
       if ( multiPos.empty() )
@@ -3160,7 +3275,7 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
       for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
       {
         const SMDS_MeshElement* newFace = fIt->next();
-        const int             faceIndex = newFace->GetID();
+        const smIdType        faceIndex = newFace->GetID();
         const gp_XYZ&           oldNorm = normals[ faceIndex ];
         if ( !SMESH_MeshAlgos::FaceNormal( newFace, faceNorm, /*normalize=*/false ) ||
              //faceNorm * moveVec < 0 )
@@ -3187,6 +3302,8 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
     // mark all new nodes located closer than theOffset from theSrcMesh
   }
 
+  removeSmallFaces( newMesh, theNew2OldFaces, tol*tol );
+
   // ==================================================
   // find self-intersections of new faces and fix them
   // ==================================================
@@ -3197,8 +3314,9 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
   Intersector intersector( newMesh, tol, normals );
 
   std::vector< const SMDS_MeshElement* > closeFaces;
-  std::vector< const SMDS_MeshNode* >    faceNodes;
+  std::vector< SMESH_NodeXYZ >           faceNodes;
   Bnd_B3d faceBox;
+
   for ( size_t iF = 1; iF < theNew2OldFaces.size(); ++iF )
   {
     const SMDS_MeshElement* newFace = theNew2OldFaces[iF].first;
@@ -3213,7 +3331,7 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
     closeFaces.clear();
     faceBox.Clear();
     for ( size_t i = 0; i < faceNodes.size(); ++i )
-      faceBox.Add( SMESH_NodeXYZ( faceNodes[i] ));
+      faceBox.Add( faceNodes[i] );
     faceBox.Enlarge( tol );
 
     fSearcher->GetElementsInBox( faceBox, SMDSAbs_Face, closeFaces );
@@ -3229,7 +3347,7 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
       // do not intersect connected faces if they have no concave nodes
       int nbCommonNodes = 0;
       for ( size_t iN = 0; iN < faceNodes.size(); ++iN )
-        nbCommonNodes += ( closeFace->GetNodeIndex( faceNodes[iN] ) >= 0 );
+        nbCommonNodes += ( closeFace->GetNodeIndex( faceNodes[iN].Node() ) >= 0 );
 
       if ( !isConcaveNode1 )
       {
@@ -3239,13 +3357,16 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
             break;
 
         if ( !isConcaveNode2 && nbCommonNodes > 0 )
-          continue;
+        {
+          if ( normals[ newFace->GetID() ] * normals[ closeFace->GetID() ] < 1.0 )
+            continue; // not co-planar
+        }
       }
 
       intersector.Cut( newFace, closeFace, nbCommonNodes );
     }
   }
-  intersector.MakeNewFaces( theNew2OldFaces, theNew2OldNodes, sign );
+  intersector.MakeNewFaces( theNew2OldFaces, theNew2OldNodes, sign, /*optimize=*/true );
 
   return newMesh;
 }