Salome HOME
[bos #40653][CEA] New mesh import export formats with meshio.
[modules/smesh.git] / src / SMESHUtils / SMESH_Delaunay.cxx
index bc5576797f0ba94bb15c793d9bc249a7c8fa88dc..14ae9faf57748f6b32c46b823ff12265ab7d0557 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2021  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2024  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
@@ -191,6 +191,11 @@ const BRepMesh_Triangle* SMESH_Delaunay::FindTriangle( const gp_XY&
 
   gp_XY uv = UV.Multiplied( _scale );
 
+  // prevent infinite loop in case of numerical instability
+  // test case NRT_GRIDS_GEOM_BUGS_15_R7
+  const BRepMesh_Triangle* tria1 = 0;
+  const BRepMesh_Triangle* tria2 = tria;
+
   while ( tria )
   {
     // check if the uv is in tria
@@ -203,7 +208,9 @@ const BRepMesh_Triangle* SMESH_Delaunay::FindTriangle( const gp_XY&
     SMESH_MeshAlgos::GetBarycentricCoords( uv,
                                            nodeUVs[0], nodeUVs[1], nodeUVs[2],
                                            bc[0], bc[1] );
-    if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 )
+    if ( (bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1) ||
+         (tria == tria1 &&
+          bc[0] >= -1e-14 && bc[1] >= -1e-14 && bc[0] + bc[1] <= 1 + 1e-14) )
     {
       if ( _triaDS->GetNode( nodeIDs[0] ).Movability() != BRepMesh_Frontier ||
            _triaDS->GetNode( nodeIDs[1] ).Movability() != BRepMesh_Frontier ||
@@ -218,6 +225,10 @@ const BRepMesh_Triangle* SMESH_Delaunay::FindTriangle( const gp_XY&
       return tria;
     }
 
+    if (tria == tria1) return 0;
+    tria1 = tria2;
+    tria2 = tria;
+
     // look for a neighbor triangle, which is adjacent to a link intersected
     // by a segment( triangle center -> uv )
 
@@ -273,7 +284,7 @@ const BRepMesh_Triangle* SMESH_Delaunay::GetTriangleNear( int iBndNode )
   if ( iBndNode >= _triaDS->NbNodes() )
     return 0;
   int nodeIDs[3];
-  int nbNbNodes = _bndNodes.size();
+  int nbBndNodes = _bndNodes.size();
 #if OCC_VERSION_LARGE <= 0x07030000
   typedef BRepMesh::ListOfInteger TLinkList;
 #else
@@ -289,9 +300,9 @@ const BRepMesh_Triangle* SMESH_Delaunay::GetTriangleNear( int iBndNode )
       if ( tria.Movability() != BRepMesh_Deleted )
       {
         _triaDS->ElementNodes( tria, nodeIDs );
-        if ( nodeIDs[0]-1 < nbNbNodes &&
-             nodeIDs[1]-1 < nbNbNodes &&
-             nodeIDs[2]-1 < nbNbNodes )
+        if ( nodeIDs[0]-1 < nbBndNodes &&
+             nodeIDs[1]-1 < nbBndNodes &&
+             nodeIDs[2]-1 < nbBndNodes )
           return &tria;
       }
     }
@@ -301,9 +312,9 @@ const BRepMesh_Triangle* SMESH_Delaunay::GetTriangleNear( int iBndNode )
       if ( tria.Movability() != BRepMesh_Deleted )
       {
         _triaDS->ElementNodes( tria, nodeIDs );
-        if ( nodeIDs[0]-1 < nbNbNodes &&
-             nodeIDs[1]-1 < nbNbNodes &&
-             nodeIDs[2]-1 < nbNbNodes )
+        if ( nodeIDs[0]-1 < nbBndNodes &&
+             nodeIDs[1]-1 < nbBndNodes &&
+             nodeIDs[2]-1 < nbBndNodes )
           return &tria;
       }
     }
@@ -373,14 +384,29 @@ void SMESH_Delaunay::ToPython() const
   }
 
   int nodeIDs[3];
+  const char* dofName[] = { "Free",
+                            "InVolume",
+                            "OnSurface",
+                            "OnCurve",
+                            "Fixed",
+                            "Frontier",
+                            "Deleted" };
+  text << "# nb elements = " << _triaDS->NbElements() << endl;
+  std::vector< int > deletedElems;
   for ( int i = 0; i < _triaDS->NbElements(); ++i )
   {
     const BRepMesh_Triangle& t = _triaDS->GetElement( i+1 );
     if ( t.Movability() == BRepMesh_Deleted )
-      continue;
+      deletedElems.push_back( i+1 );
+    //   continue;
     _triaDS->ElementNodes( t, nodeIDs );
-    text << "mesh.AddFace([ " << nodeIDs[0] << ", " << nodeIDs[1] << ", " << nodeIDs[2] << " ])" << endl;
+    text << "mesh.AddFace([ " << nodeIDs[0] << ", " << nodeIDs[1] << ", " << nodeIDs[2] << " ]) # "
+         <<  dofName[ t.Movability() ] << endl;
   }
+  text << "mesh.MakeGroupByIds( 'deleted elements', SMESH.FACE, [";
+  for ( int id : deletedElems )
+    text << id << ",";
+  text << "])" << endl;
 
   const char* fileName = "/tmp/Delaunay.py";
   SMESH_File file( fileName, false );