# patch with leaving old behavior, and adding new behavior by option

src/BRepMesh/BRepMesh_DiscretFactory.hxx

@@ -89,9 +89,11 @@ public:
   //! @param theShape shape to be meshed.
   //! @param theLinDeflection linear deflection to be used for meshing.
   //! @param theAngDeflection angular deflection to be used for meshing.
+  //! @param theTDOldBehavior tangential deflection behavior for default algorithm.
   Standard_EXPORT Handle(BRepMesh_DiscretRoot) Discret(const TopoDS_Shape& theShape,
                                                        const Standard_Real theLinDeflection,
-                                                       const Standard_Real theAngDeflection);
+                                                       const Standard_Real theAngDeflection,
+                                                       const Standard_Boolean theTDOldBehavior = Standard_True);
 
 protected:
   

src/BRepMesh/BRepMesh_FastDiscret.hxx

@@ -68,7 +68,8 @@ public:
      Relative(Standard_False),
      AdaptiveMin(Standard_False),
      InternalVerticesMode(Standard_True),
-     ControlSurfaceDeflection(Standard_True)
+     ControlSurfaceDeflection(Standard_True),
+     TDOldBehavior(Standard_True)
      {
      }
     
@@ -103,6 +104,10 @@ public:
     //! Prameter to check the deviation of triangulation and interior of
     //! the face
     Standard_Boolean                                 ControlSurfaceDeflection;
+
+    //! Prameter to set tangential deflection behavior
+    //! (old - for backward compatibility, may produce incorrect results in several cases)
+    Standard_Boolean                                 TDOldBehavior;
   };
 
 public:

src/BRepMesh/BRepMesh_IncrementalMesh.cxx

@@ -84,7 +84,8 @@ BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh( const TopoDS_Shape&    theSh
                                                     const Standard_Boolean isRelative,
                                                     const Standard_Real    theAngDeflection,
                                                     const Standard_Boolean isInParallel,
-                                                    const Standard_Boolean adaptiveMin)
+                                                    const Standard_Boolean adaptiveMin,
+                                                    const Standard_Boolean theTdOldBehavior)
 : myMaxShapeSize(0.),
   myModified(Standard_False),
   myStatus(0)
@@ -94,6 +95,7 @@ BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh( const TopoDS_Shape&    theSh
   myParameters.Angle = theAngDeflection;
   myParameters.InParallel = isInParallel;
   myParameters.AdaptiveMin = adaptiveMin;
+  myParameters.TDOldBehavior = theTdOldBehavior;
 
   myShape = theShape;
   Perform();
@@ -262,7 +264,7 @@ void BRepMesh_IncrementalMesh::discretizeFreeEdges()
     BRepAdaptor_Curve aCurve(aEdge);
     GCPnts_TangentialDeflection aDiscret(aCurve, aCurve.FirstParameter(),
       aCurve.LastParameter(), myParameters.Angle, aEdgeDeflection, 2,
-      Precision::PConfusion(), myParameters.MinSize);
+      Precision::PConfusion(), myParameters.MinSize, myParameters.TDOldBehavior);
 
     Standard_Integer aNodesNb = aDiscret.NbPoints();
     TColgp_Array1OfPnt   aNodes  (1, aNodesNb);
@@ -572,6 +574,23 @@ Standard_Integer BRepMesh_IncrementalMesh::Discret(
   return 0; // no error
 }
 
+Standard_Integer BRepMesh_IncrementalMesh::Discret(
+  const TopoDS_Shape&    theShape,
+  const Standard_Real    theDeflection,
+  const Standard_Real    theAngle,
+  const Standard_Boolean theTdOldBehavior,
+  BRepMesh_DiscretRoot* &theAlgo)
+{
+    BRepMesh_IncrementalMesh* anAlgo = new BRepMesh_IncrementalMesh();
+    anAlgo->ChangeParameters().Deflection = theDeflection;
+    anAlgo->ChangeParameters().Angle = theAngle;
+    anAlgo->ChangeParameters().InParallel = IS_IN_PARALLEL;
+    anAlgo->ChangeParameters().TDOldBehavior = theTdOldBehavior;
+    anAlgo->SetShape(theShape);
+    theAlgo = anAlgo;
+    return 0; // no error
+}
+
 //=======================================================================
 //function : IsParallelDefault
 //purpose  :

src/BRepMesh/BRepMesh_IncrementalMesh.hxx

@@ -57,7 +57,8 @@ public: //! @name mesher API
     const Standard_Boolean isRelative = Standard_False,
     const Standard_Real    theAngDeflection = 0.5,
     const Standard_Boolean isInParallel = Standard_False,
-    const Standard_Boolean adaptiveMin = Standard_False);  
+    const Standard_Boolean adaptiveMin = Standard_False,
+    const Standard_Boolean theTdOldBehavior = Standard_True);
 
   //! Constructor.
   //! Automatically calls method Perform.
@@ -109,6 +110,19 @@ public: //! @name plugin API
                                                   const Standard_Real    theAngDeflection,
                                                   BRepMesh_DiscretRoot* &theAlgo);
 
+  //! Plugin interface for the Mesh Factories.
+  //! Initializes meshing algorithm with the given parameters.
+  //! @param theShape shape to be meshed.
+  //! @param theLinDeflection linear deflection.
+  //! @param theAngDeflection angular deflection.
+  //! @param theTdOldBehavior angle comparison mode for tangential deflection
+  //! @param[out] theAlgo pointer to initialized algorithm.
+  Standard_EXPORT static Standard_Integer Discret(const TopoDS_Shape&    theShape,
+                                                  const Standard_Real    theLinDeflection,
+                                                  const Standard_Real    theAngDeflection,
+                                                  const Standard_Boolean theTdOldBehavior,
+                                                  BRepMesh_DiscretRoot* &theAlgo);
+  
   //! Returns multi-threading usage flag set by default in 
   //! Discret() static method (thus applied only to Mesh Factories).
   Standard_EXPORT static Standard_Boolean IsParallelDefault();

src/GCPnts/GCPnts_TangentialDeflection.cxx

@@ -71,7 +71,22 @@ static Standard_Real EstimAngl(const gp_Pnt& P1, const gp_Pnt& Pm, const gp_Pnt&
   //
   if(L > gp::Resolution())
   {
-    return V1.Dot(V2)/L;
+    return V1.Dot(V2) / L;
+  }
+  else
+  {
+    return 0.;
+  }
+}
+
+static Standard_Real EstimAnglOld(const gp_Pnt& P1, const gp_Pnt& Pm, const gp_Pnt& P2)
+{
+    gp_Vec V1(P1, Pm), V2(Pm, P2);
+    Standard_Real L = Sqrt(V1.SquareMagnitude() * V2.SquareMagnitude());
+    //
+    if(L > gp::Resolution())
+    {
+        return V1.CrossMagnitude(V2) / L;
   }
   else
   {

src/GCPnts/GCPnts_TangentialDeflection.hxx

@@ -75,21 +75,21 @@ public:
   
   Standard_EXPORT GCPnts_TangentialDeflection();
   
-  Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor3d_Curve& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7);
+  Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor3d_Curve& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7, const Standard_Boolean ACompOldBehavior = Standard_False);
   
-  Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor3d_Curve& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7);
+  Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor3d_Curve& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7, const Standard_Boolean ACompOldBehavior = Standard_False);
   
-  Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor2d_Curve2d& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7);
+  Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor2d_Curve2d& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7, const Standard_Boolean ACompOldBehavior = Standard_False);
   
-  Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor2d_Curve2d& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7);
+  Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor2d_Curve2d& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7, const Standard_Boolean ACompOldBehavior = Standard_False);
   
-  Standard_EXPORT void Initialize (const Adaptor3d_Curve& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7);
+  Standard_EXPORT void Initialize (const Adaptor3d_Curve& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7, const Standard_Boolean ACompOldBehavior = Standard_False);
   
-  Standard_EXPORT void Initialize (const Adaptor3d_Curve& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7);
+  Standard_EXPORT void Initialize (const Adaptor3d_Curve& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7, const Standard_Boolean ACompOldBehavior = Standard_False);
   
-  Standard_EXPORT void Initialize (const Adaptor2d_Curve2d& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7);
+  Standard_EXPORT void Initialize (const Adaptor2d_Curve2d& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7, const Standard_Boolean ACompOldBehavior = Standard_False);
   
-  Standard_EXPORT void Initialize (const Adaptor2d_Curve2d& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7);
+  Standard_EXPORT void Initialize (const Adaptor2d_Curve2d& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7, const Standard_Boolean ACompOldBehavior = Standard_False);
   
   //! Add point to already calculated points (or replace existing)
   //! Returns index of new added point
@@ -111,6 +111,12 @@ public:
     return points.Value (I);
   }
 
+  //! Returns if the old behavior for angle comparison (based on sin) or the new one (based on cos) should be used.
+  Standard_Boolean ACompOldBehavior()
+  {
+      return aCompOldBehavior;
+  }
+  
   //! Computes angular step for the arc using the given parameters.
   Standard_EXPORT static Standard_Real ArcAngularStep (const Standard_Real theRadius, const Standard_Real theLinearDeflection, const Standard_Real theAngularDeflection, const Standard_Real theMinLength);
 
@@ -148,6 +154,7 @@ private:
   Standard_Real firstu;
   TColgp_SequenceOfPnt points;
   TColStd_SequenceOfReal parameters;
+  Standard_Boolean aCompOldBehavior;
 };
 
 #endif // _GCPnts_TangentialDeflection_HeaderFile

src/GCPnts/GCPnts_TangentialDeflection.pxx

@@ -42,7 +42,7 @@ void GCPnts_TangentialDeflection::EvaluateDu (
     Standard_Real Lc = N.CrossMagnitude (T);
     Standard_Real Ln = Lc/Lt;
     if (Ln > LTol) {
-      Du = sqrt (8.0 * Max(curvatureDeflection, myMinLen) / Ln);
+      Du = Sqrt (8.0 * Max(curvatureDeflection, myMinLen) / Ln);
       NotDone = Standard_False;
     }
   }
@@ -60,10 +60,11 @@ GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
  const Standard_Real    CurvatureDeflection,
  const Standard_Integer MinimumOfPoints,
  const Standard_Real    UTol,
- const Standard_Real    theMinLen)
+ const Standard_Real    theMinLen, 
+ const Standard_Boolean ACompOldBehavior)
 
 { 
-  Initialize (C,AngularDeflection,CurvatureDeflection,MinimumOfPoints,UTol,theMinLen); 
+  Initialize (C,AngularDeflection,CurvatureDeflection,MinimumOfPoints,UTol,theMinLen, ACompOldBehavior);
 }
 
 
@@ -80,7 +81,8 @@ GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
  const Standard_Real    CurvatureDeflection,
  const Standard_Integer MinimumOfPoints,
  const Standard_Real    UTol,
- const Standard_Real    theMinLen)
+ const Standard_Real    theMinLen,
+ const Standard_Boolean ACompOldBehavior)
 
 { 
   Initialize (C, 
@@ -89,7 +91,8 @@ GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
         AngularDeflection, 
         CurvatureDeflection, 
         MinimumOfPoints,
-        UTol, theMinLen);
+        UTol, theMinLen,
+        ACompOldBehavior);
 }
 
 
@@ -105,7 +108,8 @@ void GCPnts_TangentialDeflection::Initialize (
  const Standard_Real    CurvatureDeflection,
  const Standard_Integer MinimumOfPoints,
  const Standard_Real    UTol,
- const Standard_Real    theMinLen)
+ const Standard_Real    theMinLen,
+ const Standard_Boolean ACompOldBehavior)
 
 {
   Initialize (C, 
@@ -114,7 +118,8 @@ void GCPnts_TangentialDeflection::Initialize (
               AngularDeflection, 
               CurvatureDeflection,
               MinimumOfPoints,
-              UTol, theMinLen);
+              UTol, theMinLen,
+              ACompOldBehavior);
 }
 
 
@@ -131,12 +136,14 @@ void GCPnts_TangentialDeflection::Initialize (
  const Standard_Real    CurvatureDeflection,
  const Standard_Integer MinimumOfPoints,
  const Standard_Real    UTol,
- const Standard_Real    theMinLen)
+ const Standard_Real    theMinLen,
+ const Standard_Boolean ACompOldBehavior)
 
 {
   
   Standard_ConstructionError_Raise_if (CurvatureDeflection <= Precision::Confusion () || AngularDeflection   <= Precision::Angular (), "GCPnts_TangentialDeflection::Initialize - Zero Deflection")
 
+ aCompOldBehavior = ACompOldBehavior;
  parameters.Clear ();
  points    .Clear ();
  if (FirstParameter < LastParameter) {
@@ -462,7 +469,7 @@ void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
       L2 = V2.SquareModulus ();
       if (L1 > MinLen2 && L2 > MinLen2)
       {
-        Standard_Real angg = V1.CrossMagnitude(V2) / (Sqrt(L1 * L2));
+        Standard_Real angg = V1.CrossMagnitude(V2) / Sqrt(L1 * L2);
         ACoef = angg / AngleMax;
       }
       else
@@ -607,7 +614,7 @@ void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
   //Additional check for intervals
   Standard_Integer MaxNbp = 10 * Nbp;
 
-  Standard_Real CosAngleMax = Sqrt(1.0 - AngleMax * AngleMax);
+  const Standard_Real AngleMaxFixed = aCompOldBehavior ? AngleMax : Sqrt(1.0 - AngleMax * AngleMax);
 
   for(i = 1; i < Nbp; ++i)
   {
@@ -616,14 +623,16 @@ void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
     // Check maximal deflection on interval;
     Standard_Real dmax = 0.;
     Standard_Real umax = 0.;
-    Standard_Real amax = 0.;
     EstimDefl(C, U1, U2, dmax, umax);
     const gp_Pnt& P1 = points(i);
     const gp_Pnt& P2 = points(i+1);
     D0(C, umax, MiddlePoint);
-    amax = EstimAngl(P1, MiddlePoint, P2);
+
-    if(dmax > curvatureDeflection || amax < CosAngleMax)
+    Standard_Real amax = aCompOldBehavior ? EstimAnglOld(P1, MiddlePoint, P2) : EstimAngl(P1, MiddlePoint, P2);
+
+    if(dmax > curvatureDeflection || (aCompOldBehavior ? amax - AngleMaxFixed : AngleMaxFixed - amax) > ATol)
     {
+      
       if(umax - U1 > uTol && U2 - umax > uTol)
       {
         if (P1.SquareDistance(MiddlePoint) > MinLen2 && P2.SquareDistance(MiddlePoint) > MinLen2)
@@ -660,7 +669,7 @@ void GCPnts_TangentialDeflection::EstimDefl (const TheCurve& C,
   Standard_Real reltol = Max(1.e-3, 2.*uTol/((Abs(U1) + Abs(U2))));
   //
   math_BrentMinimum anOptLoc(reltol, aNbIter, uTol);
-  anOptLoc.Perform(aFunc, U1, (U1+U2)/2., U2);
+  anOptLoc.Perform(aFunc, U1, (U1+U2)*0.5, U2);
   if(anOptLoc.IsDone())
   {
     MaxDefl = Sqrt(-anOptLoc.Minimum());

src/IVtkOCC/IVtkOCC_ShapeMesher.cxx

@@ -128,7 +128,8 @@ void IVtkOCC_ShapeMesher::meshShape()
     Handle(BRepMesh_DiscretRoot) anAlgo;
     anAlgo = BRepMesh_DiscretFactory::Get().Discret (anOcctShape,
                                                      aDeflection,
-                                                     GetDeviationAngle());
+                                                     GetDeviationAngle(),
+                                                     GetTDOldBehavior());
     if (!anAlgo.IsNull())
     {
       anAlgo->Perform();
@@ -458,6 +459,7 @@ static void FindLimits (const Adaptor3d_Curve& theCurve,
 static void DrawCurve (Adaptor3d_Curve&       theCurve,
                        const Standard_Real    theDeflection,
                        const Standard_Real    theAngle,
+                       const Standard_Boolean theTDOldBehavior,
                        const Standard_Real    theU1,
                        const Standard_Real    theU2,
                        IVtk_Polyline&         thePoints)
@@ -493,7 +495,7 @@ static void DrawCurve (Adaptor3d_Curve&    theCurve,
           anU1 = Max(anU1, anU1);
           anU2 = Min(anU2, anU2);
 
-          GCPnts_TangentialDeflection anAlgo (theCurve, anU1, anU2, theAngle, theDeflection);
+          GCPnts_TangentialDeflection anAlgo (theCurve, anU1, anU2, theAngle, theDeflection, theTDOldBehavior);
           NumberOfPoints = anAlgo.NbPoints();
 
           if (NumberOfPoints > 0)
@@ -766,7 +768,7 @@ void IVtkOCC_ShapeMesher::buildIsoLines (const Handle(BRepAdaptor_HSurface)& the
         if (aB2 - aB1 > Precision::Confusion())
         {
           IVtk_Polyline aPoints;
-          DrawCurve (aGeomCurve, aDeflection, anAngle, aB1, aB2, aPoints);
+          DrawCurve (aGeomCurve, aDeflection, anAngle, GetTDOldBehavior(), aB1, aB2, aPoints);
           thePolylines.Append (aPoints);
         }
       }
@@ -784,7 +786,7 @@ void IVtkOCC_ShapeMesher::buildIsoLines (const Handle(BRepAdaptor_HSurface)& the
         if (aB2 - aB1>Precision::Confusion())
         {
           IVtk_Polyline aPoints;
-          DrawCurve (anIso, aDeflection, anAngle, aB1, aB2, aPoints);
+          DrawCurve (anIso, aDeflection, anAngle, GetTDOldBehavior(), aB1, aB2, aPoints);
           thePolylines.Append (aPoints);
         }
       }

src/IVtkOCC/IVtkOCC_ShapeMesher.hxx

@@ -50,9 +50,11 @@ public:
   IVtkOCC_ShapeMesher (const Standard_Real& theDevCoeff = 0.0001,
                        const Standard_Real& theDevAngle = 12.0 * M_PI / 180.0,
                        const Standard_Integer theNbUIsos = 1,
-                       const Standard_Integer theNbVIsos = 1)
+                       const Standard_Integer theNbVIsos = 1,
+                       const Standard_Boolean theTDOldBehavior = Standard_True)
  : myDevCoeff (theDevCoeff),
    myDevAngle (theDevAngle),
+   myTDOldBehavior(theTDOldBehavior),
    myDeflection (0.0)
   {
     myNbIsos[0] = theNbUIsos;
@@ -85,6 +87,13 @@ public:
     return myDevAngle;
   }
 
+  //! Returns tangential deflection behavior used in the algorithm
+  //! @return tangential deflection behavior
+  Standard_Boolean GetTDOldBehavior() const
+  {
+      return myTDOldBehavior;
+  }
+
 protected:
   //! Executes the mesh generation algorithms. To be defined in implementation class.
   Standard_EXPORT virtual void internalBuild() Standard_OVERRIDE;
@@ -186,6 +195,7 @@ private:
   Standard_Real         myDevAngle;
   mutable Standard_Real myDeflection;
   Standard_Integer      myNbIsos[2];
+  Standard_Boolean      myTDOldBehavior;
 };
 
 #endif //  __IVTKOCC_SHAPEMESHER_H__