Salome HOME
bos #32738 [CEA] Scaled Jacobian quality mesh measure for volumetric elements.
[modules/smesh.git] / src / SMESHUtils / SMESH_Offset.cxx
index 71f53213ded3a5c564c73528fa19379518f835bc..ee14065d0bda9bc37b90298077a0b830f888decc 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2019  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2023  CEA, EDF, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -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
@@ -1884,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 );
       }
@@ -1955,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 )
@@ -1970,10 +2019,15 @@ namespace SMESH_MeshAlgos
         {
           if ( nodes[i] == nodes[i+1] || nodes[i] == nodes[i+2] || nodes[i+1] == nodes[i+2] )
           {
-#ifdef _DEBUG_
-            std::cerr << "BAD tria" << std::endl;
-            cf.Dump();
-#endif
+            if (SALOME::VerbosityActivated())
+            {
+              std::cerr << "BAD tria" << std::endl;
+              cf.Dump();
+            }
+            else
+            {
+              if ( i < 0 ) cf.Dump(); // avoid "CutFace::Dump() unused in release mode"
+            }
             continue;
           }
           if (!( tria = myMesh->FindFace( nodes[i], nodes[i+1], nodes[i+2] )))
@@ -2043,7 +2097,7 @@ namespace SMESH_MeshAlgos
           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;
@@ -2075,7 +2129,7 @@ namespace SMESH_MeshAlgos
       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
 
@@ -2448,7 +2502,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 )
@@ -2707,7 +2761,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 )
       {
@@ -2869,7 +2923,7 @@ 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;
@@ -2932,8 +2986,8 @@ namespace
             bool isAdded = checkedCoplanar.insert( myLinks[iE].myFace ).second;
             if ( !isAdded )
               continue;
-            toErase = SMESH_MeshAlgos::GetCommonNodes( myLinks[i ].myFace,
-                                                       myLinks[iE].myFace ).size() < 1;
+            toErase = ( SMESH_MeshAlgos::NbCommonNodes( myLinks[i ].myFace,
+                                                        myLinks[iE].myFace ) < 1 );
           }
         }
       }
@@ -3026,7 +3080,7 @@ namespace
       //check if the faces are connected
       int nbCommonNodes = 0;
       if ( e.myFace && myFace )
-        nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e.myFace, myFace ).size();
+        nbCommonNodes = SMESH_MeshAlgos::NbCommonNodes( e.myFace, myFace );
       bool toReplace = (( myIndex == _INTERNAL && nbCommonNodes > 1 ) ||
                         ( myIndex == _COPLANAR && nbCommonNodes < 2 ));
       if ( toReplace )
@@ -3175,7 +3229,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() )
@@ -3224,7 +3278,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 )
@@ -3251,6 +3305,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
   // ==================================================
@@ -3261,8 +3317,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;
@@ -3277,7 +3334,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 );
@@ -3293,7 +3350,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 )
       {
@@ -3312,7 +3369,7 @@ SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
       intersector.Cut( newFace, closeFace, nbCommonNodes );
     }
   }
-  intersector.MakeNewFaces( theNew2OldFaces, theNew2OldNodes, sign );
+  intersector.MakeNewFaces( theNew2OldFaces, theNew2OldNodes, sign, /*optimize=*/true );
 
   return newMesh;
 }