0026407: BRepMesh_Delaun should not take into account frontier edges on first pass of algorithm

Correction of new vertices insertion: remove all triangles shot by point even if they contain frontier edges in order to prevent cases of free edge glued with frontier;
findNextPolygonLink: choose link with opposite direction to previous one as the last resort in case if the is no another option;
Enlarge bounding boxes by Precision::PConfusion() in order to not to miss possible intersections;

Test-case for issue #26407

Return NoIntersection in case of end point touch and isConsiderEndPointTouch flag is not set

Update of test-cases according to the new behavior

src/BRepMesh/BRepMesh_Delaun.cxx

@@ -69,6 +69,13 @@ namespace {
   private:
     Handle(BRepMesh_DataStructureOfDelaun) myStructure;
   };
+
+  inline void UpdateBndBox(const gp_XY& thePnt1, const gp_XY& thePnt2, Bnd_B2d& theBox)
+  {
+    theBox.Add( thePnt1 );
+    theBox.Add( thePnt2 );
+    theBox.Enlarge(Precision);
+  }
 } // anonymous namespace
 
 //=======================================================================
@@ -424,6 +431,10 @@ void BRepMesh_Delaun::createTrianglesOnNewVertices(
 
   BRepMesh::MapOfIntegerInteger aLoopEdges(10, aAllocator);
 
+  Standard_Real aTolU, aTolV;
+  myMeshData->Data()->GetTolerance(aTolU, aTolV);
+  const Standard_Real aSqTol = aTolU * aTolU + aTolV * aTolV;
+
   // Insertion of nodes :
   Standard_Boolean isModify = Standard_True;
   
@@ -447,20 +458,23 @@ void BRepMesh_Delaun::createTrianglesOnNewVertices(
       // To add a node in the mesh it is necessary to check conditions: 
       // - the node should be within the boundaries of the mesh and so in an existing triangle
       // - all adjacent triangles should belong to a component connected with this triangle
-      if ( Contains( aCircleIt.Value(), aVertex, onEgdeId ) )
+      if ( Contains( aCircleIt.Value(), aVertex, aSqTol, onEgdeId ) )
       {
-        if ( onEgdeId == 0 )
+        if (onEgdeId != 0 && GetEdge(onEgdeId).Movability() != BRepMesh_Free)
         {
-          aTriangleId = aCircleIt.Value();
+          // We can skip free vertex too close to the frontier edge.
-          aCirclesList.Remove( aCircleIt );
+          if (aVertex.Movability() == BRepMesh_Free)
-          break;
+            continue;
-        }
+
-        else if ( GetEdge( onEgdeId ).Movability() == BRepMesh_Free )
+          // However, we should add vertex that have neighboring frontier edges.
-        {
-          aTriangleId = aCircleIt.Value();
-          aCirclesList.Remove( aCircleIt );
-          break;
         }
+
+        // Remove triangle even if it contains frontier edge in order 
+        // to prevent appearance of incorrect configurations like free 
+        // edge glued with frontier #26407
+        aTriangleId = aCircleIt.Value();
+        aCirclesList.Remove( aCircleIt );
+        break;
       }
     }
 
@@ -779,8 +793,7 @@ void BRepMesh_Delaun::fillBndBox(BRepMesh::SequenceOfBndB2d& theBoxes,
                                  const BRepMesh_Vertex&      theV2)
 {
   Bnd_B2d aBox;
-  aBox.Add( theV1.Coord() );
+  UpdateBndBox(theV1.Coord(), theV2.Coord(), aBox);
-  aBox.Add( theV2.Coord() );
   theBoxes.Append( aBox );
 }
 
@@ -1035,8 +1048,8 @@ Standard_Boolean BRepMesh_Delaun::checkIntersection(
   const Standard_Boolean             isSkipLastEdge,
   Bnd_B2d&                           theLinkBndBox ) const
 {
-  theLinkBndBox.Add( GetVertex( theLink.FirstNode() ).Coord() );
+  UpdateBndBox(GetVertex(theLink.FirstNode()).Coord(),
-  theLinkBndBox.Add( GetVertex( theLink.LastNode()  ).Coord() );
+    GetVertex(theLink.LastNode()).Coord(), theLinkBndBox);
 
   Standard_Integer aPolyLen = thePolygon.Length();
   // Don't check intersection with the last link
@@ -1457,8 +1470,7 @@ Standard_Integer BRepMesh_Delaun::createAndReplacePolygonLink(
     theNodes[0], theNodes[1], BRepMesh_Free ) );
 
   Bnd_B2d aNewBox;
-  aNewBox.Add( thePnts[0] );
+  UpdateBndBox(thePnts[0].Coord(), thePnts[1].Coord(), aNewBox);
-  aNewBox.Add( thePnts[1] );
 
   switch ( theReplaceFlag )
   {
@@ -1858,8 +1870,7 @@ void BRepMesh_Delaun::decomposeSimplePolygon(
       const gp_Pnt2d&         aLinkFirstVertex = aRefVertices[aRefLinkNodeIt];
 
       Bnd_B2d aBox;
-      aBox.Add( aLinkFirstVertex );
+      UpdateBndBox(aLinkFirstVertex.Coord(), aPivotVertex.Coord(), aBox);
-      aBox.Add( aPivotVertex );
 
       BRepMesh_Edge aCheckLink( aLinkFirstNode, aPivotNode, BRepMesh_Free );
 
@@ -1943,8 +1954,7 @@ void BRepMesh_Delaun::decomposeSimplePolygon(
     thePolyBoxes.Split(aUsedLinkId, thePolyBoxesCut);
 
     Bnd_B2d aBox;
-    aBox.Add( aRefVertices[0] );
+    UpdateBndBox(aRefVertices[0].Coord(), aRefVertices[2].Coord(), aBox);
-    aBox.Add( aRefVertices[2] );
     thePolyBoxesCut.Prepend( aBox );
   }
   else
@@ -1958,9 +1968,7 @@ void BRepMesh_Delaun::decomposeSimplePolygon(
     thePolygon.SetValue( 1, -aNewEdgesInfo[1] );
 
     Bnd_B2d aBox;
-    aBox.Add( aRefVertices[1] );
+    UpdateBndBox(aRefVertices[1].Coord(), aRefVertices[2].Coord(), aBox);
-    aBox.Add( aRefVertices[2] );
-
     thePolyBoxes.SetValue( 1, aBox );
   }
 }
@@ -2189,17 +2197,12 @@ BRepMesh::HMapOfInteger BRepMesh_Delaun::getEdgesByType(
 //=======================================================================
 Standard_Real BRepMesh_Delaun::calculateDist( const gp_XY            theVEdges[3],
                                               const gp_XY            thePoints[3],
-                                              const Standard_Integer theEdgesId[3],
                                               const BRepMesh_Vertex& theVertex,
                                               Standard_Real          theDistance[3],
                                               Standard_Real          theSqModulus[3],
                                               Standard_Integer&      theEdgeOn ) const
 {
-  Standard_Real aMinDist = -1;
+  Standard_Real aMinDist = RealLast();
-  if ( !theVEdges || !thePoints || !theEdgesId || 
-       !theDistance || !theSqModulus )
-    return aMinDist;
-    
   for( Standard_Integer i = 0; i < 3; ++i )
   {
     theSqModulus[i] = theVEdges[i].SquareModulus();
@@ -2211,9 +2214,9 @@ Standard_Real BRepMesh_Delaun::calculateDist( const gp_XY            theVEdges[3
     Standard_Real aDist = theDistance[i] * theDistance[i];
     aDist /= theSqModulus[i];
     
-    if ( aMinDist < 0 || aDist < aMinDist )
+    if ( aDist < aMinDist )
     {
-      theEdgeOn = theEdgesId[i];
+      theEdgeOn = i;
       aMinDist  = aDist;
     }
   }
@@ -2229,7 +2232,8 @@ Standard_Real BRepMesh_Delaun::calculateDist( const gp_XY            theVEdges[3
 //=======================================================================
 Standard_Boolean BRepMesh_Delaun::Contains( const Standard_Integer theTriangleId,
                                             const BRepMesh_Vertex& theVertex,
-                                            Standard_Integer&      theEdgeOn ) const
+                                            const Standard_Real    theSqTolerance,
+                                            Standard_Integer&      theEdgeOn) const
 {
   theEdgeOn = 0;
   
@@ -2264,34 +2268,24 @@ Standard_Boolean BRepMesh_Delaun::Contains( const Standard_Integer theTriangleId
   Standard_Real aDistance[3];
   Standard_Real aSqModulus[3];
 
-  Standard_Real aMinDist;  
+  Standard_Real aSqMinDist;
-  aMinDist = calculateDist( aVEdges, aPoints, e, theVertex, aDistance, aSqModulus, theEdgeOn );
+  Standard_Integer aEdgeOnId;
-  if ( aMinDist < 0 )
+  aSqMinDist = calculateDist( aVEdges, aPoints, theVertex, aDistance, aSqModulus, aEdgeOnId );
+  if ( aSqMinDist < 0 )
     return Standard_False;
 
-  if ( aMinDist > Precision2 )
+  const Standard_Boolean isNotFree = (anEdges[aEdgeOnId]->Movability() != BRepMesh_Free);
+  if ( aSqMinDist > theSqTolerance )
   {
-    Standard_Integer anEdgeId = theEdgeOn;
+    if (isNotFree && aDistance[aEdgeOnId] < ( aSqModulus[aEdgeOnId] / 5. ))
-    theEdgeOn = 0;
+      theEdgeOn = e[aEdgeOnId];
-    
-    if ( anEdgeId != 0 ) 
-    {
-      Standard_Integer i = 0;
-      for ( ; i < 3; ++i )
-      {
-        if( e[i] == anEdgeId )
-          break;
-      }
-      
-      if( anEdges[i]->Movability() != BRepMesh_Free )
-        if ( aDistance[i] < ( aSqModulus[i] / 5. ) )
-          theEdgeOn = e[i];
-    }
   }
+  else if (isNotFree)
+    return Standard_False;
+  else
+    theEdgeOn = e[aEdgeOnId];
 
-  return ( aDistance[0] + aDistance[1] + aDistance[2] != 0. &&
+  return (aDistance[0] >= 0. && aDistance[1] >= 0. && aDistance[2] >= 0.);
-            ( ( aDistance[0] >= 0. && aDistance[1] >= 0. && aDistance[2] >= 0. ) ||
-              ( aDistance[0] <= 0. && aDistance[1] <= 0. && aDistance[2] <= 0. ) ) );
 }
 
 //=============================================================================
@@ -2355,3 +2349,71 @@ Standard_Real BRepMesh_Delaun::polyArea(const BRepMesh::SequenceOfInteger& thePo
 
   return aArea / 2.;
 }
+
+#ifdef DEB
+//=======================================================================
+//function : BRepMesh_DumpPoly
+//purpose  : 
+//=======================================================================
+#include <TopoDS_Compound.hxx>
+#include <BRep_Builder.hxx>
+#include <Standard_ErrorHandler.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepTools.hxx>
+Standard_CString BRepMesh_DumpPoly(void*            thePolygon,
+                                   void*            theMeshHandlePtr,
+                                   Standard_CString theFileNameStr)
+{
+  if (thePolygon == 0 || theFileNameStr == 0)
+  {
+    return "Error: file name or polygon data is null";
+  }
+
+  BRepMesh::SequenceOfInteger& aPolygon = *(BRepMesh::SequenceOfInteger*)thePolygon;
+
+  Handle(BRepMesh_DataStructureOfDelaun) aMeshData = 
+    *(Handle(BRepMesh_DataStructureOfDelaun)*)theMeshHandlePtr;
+
+  if (aMeshData.IsNull())
+    return "Error: mesh data is empty";
+
+  TopoDS_Compound aMesh;
+  BRep_Builder aBuilder;
+  aBuilder.MakeCompound(aMesh);
+
+  try
+  {
+    OCC_CATCH_SIGNALS
+
+    BRepMesh::SequenceOfInteger::Iterator aLinksIt(aPolygon);
+    for (; aLinksIt.More(); aLinksIt.Next())
+    {
+      const BRepMesh_Edge& aLink = aMeshData->GetLink(Abs(aLinksIt.Value()));
+
+      gp_Pnt aPnt[2];
+      for (Standard_Integer i = 0; i < 2; ++i)
+      {
+        const Standard_Integer aNodeId = 
+          (i == 0) ? aLink.FirstNode() : aLink.LastNode();
+
+        const gp_XY& aNode = aMeshData->GetNode(aNodeId).Coord();
+        aPnt[i] = gp_Pnt(aNode.X(), aNode.Y(), 0.);
+      }
+
+      if (aPnt[0].SquareDistance(aPnt[1]) < Precision::SquareConfusion())
+        continue;
+
+      aBuilder.Add(aMesh, BRepBuilderAPI_MakeEdge(aPnt[0], aPnt[1]));
+    }
+
+    if (!BRepTools::Write(aMesh, theFileNameStr))
+      return "Error: write failed";
+  }
+  catch (Standard_Failure)
+  {
+    return Standard_Failure::Caught()->GetMessageString();
+  }
+
+  return theFileNameStr;
+}
+#endif

src/BRepMesh/BRepMesh_Delaun.hxx

@@ -114,10 +114,12 @@ public:
   }
 
   //! Test is the given triangle contains the given vertex.
-  //! If theEdgeOn != 0 the vertex lies onto the edge index
+  //! @param theSqTolerance square tolerance to check closeness to some edge
-  //! returned through this parameter.
+  //! @param theEdgeOn If it is != 0 the vertex lies onto the edge index
+  //!        returned through this parameter.
   Standard_EXPORT Standard_Boolean Contains (const Standard_Integer theTriangleId,
                                              const BRepMesh_Vertex& theVertex,
+                                             const Standard_Real    theSqTolerance,
                                              Standard_Integer&      theEdgeOn) const;
 
 private:
@@ -300,7 +302,6 @@ private:
   //! Calculates distances between the given point and edges of triangle.
   Standard_Real calculateDist (const gp_XY            theVEdges[3],
                                const gp_XY            thePoints[3],
-                               const Standard_Integer theEdgesId[3],
                                const BRepMesh_Vertex& theVertex,
                                Standard_Real          theDistance[3],
                                Standard_Real          theSqModulus[3],

src/BRepMesh/BRepMesh_FastDiscret.cxx

@@ -444,9 +444,13 @@ Standard_Integer BRepMesh_FastDiscret::Add(const TopoDS_Face& theFace)
       Standard_Real deltaY = 1.0;
 
       {
+        Standard_Real aTolU, aTolV;
+        myAttribute->ChangeStructure()->Data()->GetTolerance(aTolU, aTolV);
+        const Standard_Real aTol = Sqrt(aTolU * aTolU + aTolV * aTolV);
+
         BRepMesh::HClassifier& aClassifier = myAttribute->ChangeClassifier();
-        BRepMesh_WireChecker aDFaceChecker(aFace, Precision::PConfusion(),
+        BRepMesh_WireChecker aDFaceChecker(aFace, aTol, aInternalEdges, 
-          aInternalEdges, aVertexEdgeMap, myAttribute->ChangeStructure(),
+          aVertexEdgeMap, myAttribute->ChangeStructure(),
           myumin, myumax, myvmin, myvmax, myInParallel );
 
         aDFaceChecker.ReCompute(aClassifier);

src/BRepMesh/BRepMesh_GeomTool.cxx

@@ -223,10 +223,12 @@ BRepMesh_GeomTool::IntFlag BRepMesh_GeomTool::IntSegSeg(
   };
 
   // Consider case when edges have shared vertex
-  if ( isConsiderEndPointTouch )
+  if ( aPointHash[0] < 0 || aPointHash[1] < 0 )
   {
-    if ( aPointHash[0] < 0 || aPointHash[1] < 0 )
+    if ( isConsiderEndPointTouch )
       return BRepMesh_GeomTool::EndPointTouch;
+
+    return BRepMesh_GeomTool::NoIntersection;
   }
 
   Standard_Integer aPosHash = 

src/BRepMesh/BRepMesh_VertexTool.hxx

@@ -80,6 +80,17 @@ public:
     myTolerance[1] = theToleranceY;
   }
 
+  //! Gets the tolerance to be used for identification of 
+  //! coincident vertices.
+  //! @param theToleranceX tolerance for X dimension.
+  //! @param theToleranceY tolerance for Y dimension.
+  Standard_EXPORT void GetTolerance(Standard_Real& theToleranceX,
+                                    Standard_Real& theToleranceY)
+  {
+    theToleranceX = myTolerance[0];
+    theToleranceY = myTolerance[1];
+  }
+
   //! Adds vertex with empty data to the tool.
   //! @param theVertex node to be added to the mesh.
   //! @param isForceAdd adds the given node to structure without 

src/QABugs/QABugs_19.cxx

@@ -3679,6 +3679,74 @@ static Standard_Integer OCC26448 (Draw_Interpretor& theDI, Standard_Integer, con
   return 0;
 }
 
+//=======================================================================
+//function : OCC26407
+//purpose  :
+//=======================================================================
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepMesh_IncrementalMesh.hxx>
+#include <TCollection_AsciiString.hxx>
+static Standard_Integer OCC26407 (Draw_Interpretor& theDI, Standard_Integer theArgNb, const char** theArgVec)
+{
+  if (theArgNb != 2)
+  {
+    std::cerr << "Error: wrong number of arguments! See usage:\n";
+    theDI.PrintHelp (theArgVec[0]);
+    return 1;
+  }
+
+  // Construct vertices.
+  std::vector<TopoDS_Vertex> wire_vertices;
+  wire_vertices.push_back(BRepBuilderAPI_MakeVertex(gp_Pnt(587.90000000000009094947, 40.6758179230516248026106, 88.5)));
+  wire_vertices.push_back(BRepBuilderAPI_MakeVertex(gp_Pnt(807.824182076948432040808, 260.599999999999965893949, 88.5)));
+  wire_vertices.push_back(BRepBuilderAPI_MakeVertex(gp_Pnt(644.174182076948454778176, 424.249999999999943156581, 88.5000000000000142108547)));
+  wire_vertices.push_back(BRepBuilderAPI_MakeVertex(gp_Pnt(629.978025792618950617907, 424.25, 88.5)));
+  wire_vertices.push_back(BRepBuilderAPI_MakeVertex(gp_Pnt(793.628025792618700506864, 260.599999999999852207111, 88.5)));
+  wire_vertices.push_back(BRepBuilderAPI_MakeVertex(gp_Pnt(587.900000000000204636308, 54.8719742073813492311274, 88.5)));
+  wire_vertices.push_back(BRepBuilderAPI_MakeVertex(gp_Pnt(218.521974207381418864315, 424.250000000000056843419, 88.5)));
+  wire_vertices.push_back(BRepBuilderAPI_MakeVertex(gp_Pnt(204.325817923051886282337, 424.249999999999943156581, 88.5)));
+
+  // Construct wire.
+  BRepBuilderAPI_MakeWire wire_builder;
+  for (size_t i = 0; i < wire_vertices.size(); i++)
+  {
+    const TopoDS_Vertex &v = wire_vertices[i];
+    const TopoDS_Vertex &w = wire_vertices[(i+1) % wire_vertices.size()];
+
+    wire_builder.Add(BRepBuilderAPI_MakeEdge(v, w));
+  }
+
+  // Create face and triangulate it.
+  // Construct face.
+  gp_Pnt v0 = BRep_Tool::Pnt(wire_vertices[0]);
+  gp_Pnt v1 = BRep_Tool::Pnt(wire_vertices[1]);
+  gp_Pnt v2 = BRep_Tool::Pnt(wire_vertices[wire_vertices.size() - 1]);
+
+  gp_Vec face_normal = gp_Vec(v0, v1).Crossed(gp_Vec(v0, v2));
+
+  TopoDS_Face face = BRepBuilderAPI_MakeFace(gp_Pln(v0, face_normal), wire_builder);
+  BRepMesh_IncrementalMesh m(face, 1e-7);
+
+  if (m.GetStatusFlags() != 0)
+  {
+    theDI << "Failed. Status for face constructed from vertices: " << m.GetStatusFlags() << "\n";
+    return 1;
+  }
+  DBRep::Set(theArgVec[1], face);
+  char buf[256];
+  sprintf(buf, "isos %s 0", theArgVec[1]);
+  theDI.Eval(buf);
+
+  sprintf(buf, "triangles %s", theArgVec[1]);
+  theDI.Eval(buf);
+
+  theDI.Eval("smallview; fit");
+
+  theDI << "Test completed\n";
+  return 0;
+}
+
 void QABugs::Commands_19(Draw_Interpretor& theCommands) {
   const char *group = "QABugs";
 
@@ -3751,5 +3819,6 @@ void QABugs::Commands_19(Draw_Interpretor& theCommands) {
   theCommands.Add ("OCC26139", "OCC26139 [-boxsize value] [-boxgrid value] [-compgrid value]", __FILE__, OCC26139, group);
   theCommands.Add ("OCC26284", "OCC26284", __FILE__, OCC26284, group);
   theCommands.Add ("OCC26448", "OCC26448: check method Prepend() of sequence", __FILE__, OCC26448, group);
+  theCommands.Add ("OCC26407", "OCC26407 result_name", __FILE__, OCC26407, group);
   return;
 }

tests/bugs/mesh/bug26407

@@ -0,0 +1,11 @@
+puts "========"
+puts "OCC26407"
+puts "========"
+puts ""
+##########################################################################################
+# BRepMesh_Delaun should not take into account frontier edges on first pass of algorithm
+##########################################################################################
+
+pload QAcommands
+
+OCC26407 result

tests/mesh/data/standard/U7

@@ -7,16 +7,18 @@ set bug_freenodes "OCC22687"
 
 if { [string compare $command "shading"] == 0 } {
 ##set nbt 14
-  set nbt 9
+  set nbt 8
+  set nbl 8
   set nbn 83
   set nbwithouttri($env(os_type)) $nbt
+  set nbfree($env(os_type)) $nbl
   set nbfreenodes($env(os_type)) $nbn
 } else {
   set bug_withouttri "OCC23105"
 ##set nbt 14
   set nbt 8
   set nbn 60
-  set nbl 3
+  set nbl 12
   set nbwithouttri($env(os_type)) $nbt
   set nbfree($env(os_type)) $nbl
   set nbfreenodes($env(os_type)) $nbn

tests/mesh/data/standard/V4

@@ -2,6 +2,6 @@ set TheFileName shading_wrongshape_015.brep
 set bug_withouttri "OCC22687"
 if { [string compare $command "shading"] != 0 } {
    set bug_freenodes "OCC23105"
-   set nbfreenodes(ALL) 4
+#   set nbfreenodes(ALL) 4
 }
 set nbwithouttri(All) 7