0028131: BRepOffset_MakeOffset can't create offset with a face which created by filli...

Handling of degenerated points (with all derivatives zero) in GeomEvaluator_OffsetSurface is improved by iterative movement towards middle of the surface (extension of one-step approach implemented earlier in #28112).

diff --git a/src/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx b/src/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx
index 47b1fac738aa46a65d68ec6424eae742595be91a..c41d9a7efb63ce4aaa64f491f4993a25d2a64b30 100644 (file)
--- a/src/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx
+++ b/src/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx
@@ -24,9 +24,78 @@
 #include <gp_Vec.hxx>
 #include <Standard_RangeError.hxx>
 
-
 IMPLEMENT_STANDARD_RTTIEXT(GeomEvaluator_OffsetSurface,GeomEvaluator_Surface)
 
+namespace {
+
+// tolerance for considering derivative to be null
+const Standard_Real the_D1MagTol = 1.e-9;
+
+// If calculation of normal fails, try shifting the point towards the center 
+// of the parametric space of the surface, in the hope that derivatives 
+// are better defined there.
+//
+// This shift is iterative, starting with Precision::PConfusion()
+// and increasing by multiple of 2 on each step.
+//
+// NB: temporarily this is made as static function and not class method, 
+// hence code duplications
+static Standard_Boolean shiftPoint (const Standard_Real theUStart, const Standard_Real theVStart,
+                                    Standard_Real& theU, Standard_Real& theV, 
+                                    const Handle(Geom_Surface)& theSurf,
+                                    const Handle(GeomAdaptor_HSurface)& theAdaptor,
+                                    const gp_Vec& theD1U, const gp_Vec& theD1V)
+{
+  // Get parametric bounds and closure status
+  Standard_Real aUMin, aUMax, aVMin, aVMax;
+  Standard_Boolean isUPeriodic, isVPeriodic;
+  if (! theSurf.IsNull())
+  {
+    theSurf->Bounds (aUMin, aUMax, aVMin, aVMax);
+    isUPeriodic = theSurf->IsUPeriodic();
+    isVPeriodic = theSurf->IsVPeriodic();
+  }
+  else
+  {
+    aUMin = theAdaptor->FirstUParameter();
+    aUMax = theAdaptor->LastUParameter();
+    aVMin = theAdaptor->FirstVParameter();
+    aVMax = theAdaptor->LastVParameter();
+    isUPeriodic = theAdaptor->IsUPeriodic();
+    isVPeriodic = theAdaptor->IsVPeriodic();
+  }
+
+  // check if either U or V is singular (normally one of them is)
+  Standard_Boolean isUSingular = (theD1U.SquareMagnitude() < the_D1MagTol * the_D1MagTol);
+  Standard_Boolean isVSingular = (theD1V.SquareMagnitude() < the_D1MagTol * the_D1MagTol);
+
+  // compute vector to shift from start point to center of the surface;
+  // if surface is periodic or singular in some direction, take shift in that direction zero
+  Standard_Real aDirU = (isUPeriodic || (isUSingular && ! isVSingular) ? 0. : 0.5 * (aUMin + aUMax) - theUStart);
+  Standard_Real aDirV = (isVPeriodic || (isVSingular && ! isUSingular) ? 0. : 0.5 * (aVMin + aVMax) - theVStart);
+  Standard_Real aDist = Sqrt (aDirU * aDirU + aDirV * aDirV);
+
+  // shift current point from its current position towards center, by value of twice
+  // current distance from it to start (but not less than Precision::PConfusion());
+  // fail if center is overpassed.
+  Standard_Real aDU = theU - theUStart;
+  Standard_Real aDV = theV - theVStart;
+  Standard_Real aStep = Max (2. * Sqrt (aDU * aDU + aDV * aDV), Precision::PConfusion());
+  if (aStep >= aDist)
+  {
+    return Standard_False;
+  }
+
+  aStep /= aDist;
+  theU += aDirU * aStep;
+  theV += aDirV * aStep;
+
+//  std::cout << "Normal calculation failed at (" << theUStart << ", " << theVStart << "), shifting to (" << theU << ", " << theV << ")" << std::endl;
+
+  return Standard_True;
+}
+
+// auxiliary function
 template<class SurfOrAdapt>
 static void derivatives(Standard_Integer theMaxOrder,
                         Standard_Integer theMinOrder,
@@ -39,9 +108,103 @@ static void derivatives(Standard_Integer theMaxOrder,
                         const Standard_Boolean theAlongV,
                         const Handle(Geom_BSplineSurface)& theL,
                         TColgp_Array2OfVec& theDerNUV,
-                        TColgp_Array2OfVec& theDerSurf);
+                        TColgp_Array2OfVec& theDerSurf)
+{
+  Standard_Integer i, j;
+  gp_Pnt P;
+  gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V;
 
+  if (theAlongU || theAlongV)
+  {
+    theMaxOrder = 0;
+    TColgp_Array2OfVec DerSurfL(0, theMaxOrder + theNU + 1, 0, theMaxOrder + theNV + 1);
+    switch (theMinOrder)
+    {
+    case 1:
+      theL->D1(theU, theV, P, DL1U, DL1V);
+      DerSurfL.SetValue(1, 0, DL1U);
+      DerSurfL.SetValue(0, 1, DL1V);
+      break;
+    case 2:
+      theL->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV);
+      DerSurfL.SetValue(1, 0, DL1U);
+      DerSurfL.SetValue(0, 1, DL1V);
+      DerSurfL.SetValue(1, 1, DL2UV);
+      DerSurfL.SetValue(2, 0, DL2U);
+      DerSurfL.SetValue(0, 2, DL2V);
+      break;
+    case 3:
+      theL->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV);
+      DerSurfL.SetValue(1, 0, DL1U);
+      DerSurfL.SetValue(0, 1, DL1V);
+      DerSurfL.SetValue(1, 1, DL2UV);
+      DerSurfL.SetValue(2, 0, DL2U);
+      DerSurfL.SetValue(0, 2, DL2V);
+      DerSurfL.SetValue(3, 0, DL3U);
+      DerSurfL.SetValue(2, 1, DL3UUV);
+      DerSurfL.SetValue(1, 2, DL3UVV);
+      DerSurfL.SetValue(0, 3, DL3V);
+      break;
+    default:
+      break;
+    }
 
+    if (theNU <= theNV)
+    {
+      for (i = 0; i <= theMaxOrder + 1 + theNU; i++)
+        for (j = i; j <= theMaxOrder + theNV + 1; j++)
+          if (i + j > theMinOrder)
+          {
+            DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
+            theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
+            if (i != j && j <= theNU + 1)
+            {
+              theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
+              DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
+            }
+          }
+    }
+    else
+    {
+      for (j = 0; j <= theMaxOrder + 1 + theNV; j++)
+        for (i = j; i <= theMaxOrder + theNU + 1; i++)
+          if (i + j > theMinOrder)
+          {
+            DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
+            theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
+            if (i != j && i <= theNV + 1)
+            {
+              theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
+              DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
+            }
+          }
+    }
+    for (i = 0; i <= theMaxOrder + theNU; i++)
+      for (j = 0; j <= theMaxOrder + theNV; j++)
+      {
+        if (theAlongU)
+          theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf));
+        if (theAlongV)
+          theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL));
+      }
+  }
+  else
+  {
+    for (i = 0; i <= theMaxOrder + theNU+ 1; i++)
+      for (j = i; j <= theMaxOrder + theNV + 1; j++)
+        if (i + j > theMinOrder)
+        {
+          theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
+          if (i != j)
+            theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
+        }
+    for (i = 0; i <= theMaxOrder + theNU; i++)
+      for (j = 0; j <= theMaxOrder + theNV; j++)
+        theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf));
+  }
+}
+
+} // end of anonymous namespace
 
 GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface(
         const Handle(Geom_Surface)& theBase,
@@ -72,89 +235,27 @@ GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface(
 {
 }
 
-// If point is on parametric boundary, and calculation of normal fails,
-// try shifting it towards the inside in the hope that derivatives 
-// are better defined there.
-//
-// NB: temporarily this is made as static function and not class method, 
-// hence code duplications
-static Standard_Boolean shiftPoint (Standard_Real& theU, Standard_Real& theV, 
-                                    const Handle(Geom_Surface)& theSurf,
-                                    const Handle(GeomAdaptor_HSurface)& theAdaptor)
-{
-  // Get parametric bounds and closure status
-  Standard_Real aUMin, aUMax, aVMin, aVMax;
-  Standard_Boolean isUPeriodic, isVPeriodic;
-  if (! theSurf.IsNull())
-  {
-    theSurf->Bounds (aUMin, aUMax, aVMin, aVMax);
-    isUPeriodic = theSurf->IsUPeriodic();
-    isVPeriodic = theSurf->IsVPeriodic();
-  }
-  else
-  {
-    aUMin = theAdaptor->FirstUParameter();
-    aUMax = theAdaptor->LastUParameter();
-    aVMin = theAdaptor->FirstVParameter();
-    aVMax = theAdaptor->LastVParameter();
-    isUPeriodic = theAdaptor->IsUPeriodic();
-    isVPeriodic = theAdaptor->IsVPeriodic();
-  }
-
-  Standard_Boolean isShifted = Standard_False;
-
-  // shift by U
-  if (! isUPeriodic && aUMax - aUMin > 2 * Precision::PConfusion())
-  {
-    if (Abs (theU - aUMin) < Precision::PConfusion())
-    {
-      theU += Precision::PConfusion();
-      isShifted = Standard_True;
-    }
-    else if (Abs (theU - aUMax) < Precision::PConfusion())
-    {
-      theU -= Precision::PConfusion();
-      isShifted = Standard_True;
-    }
-  }
-
-  // shift by V
-  if (! isVPeriodic && aVMax - aVMin > 2 * Precision::PConfusion())
-  {
-    if (Abs (theV - aVMin) < Precision::PConfusion())
-    {
-      theV += Precision::PConfusion();
-      isShifted = Standard_True;
-    }
-    else if (Abs (theV - aVMax) < Precision::PConfusion())
-    {
-      theV -= Precision::PConfusion();
-      isShifted = Standard_True;
-    }
-  }
-
-  return isShifted;
-}
-
 void GeomEvaluator_OffsetSurface::D0(
     const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue) const
 {
-  gp_Vec aD1U, aD1V;
-  BaseD1(theU, theV, theValue, aD1U, aD1V);
-  try 
+  Standard_Real aU = theU, aV = theV;
+  for (;;)
   {
-    CalculateD0(theU, theV, theValue, aD1U, aD1V);
-  } 
-  catch (Geom_UndefinedValue&)
-  {
-    // if failed at parametric boundary, try taking derivative at shifted point
-    Standard_Real aU = theU, aV = theV;
-    if (! shiftPoint (aU, aV, myBaseSurf, myBaseAdaptor))
+    gp_Vec aD1U, aD1V;
+    BaseD1 (aU, aV, theValue, aD1U, aD1V);
+    try 
     {
-      throw;
+      CalculateD0(aU, aV, theValue, aD1U, aD1V);
+      break;
+    } 
+    catch (Geom_UndefinedValue&)
+    {
+      // if failed at parametric boundary, try taking derivative at shifted point
+      if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, aD1U, aD1V))
+      {
+        throw;
+      }
     }
-    BaseD1(aU, aV, theValue, aD1U, aD1V);
-    CalculateD0(theU, theV, theValue, aD1U, aD1V);
   }
 }
 
@@ -162,22 +263,24 @@ void GeomEvaluator_OffsetSurface::D1(
     const Standard_Real theU, const Standard_Real theV,
     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V) const
 {
-  gp_Vec aD2U, aD2V, aD2UV;
-  BaseD2(theU, theV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
-  try 
+  Standard_Real aU = theU, aV = theV;
+  for (;;)
   {
-    CalculateD1(theU, theV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
-  } 
-  catch (Geom_UndefinedValue&)
-  {
-    // if failed at parametric boundary, try taking derivative at shifted point
-    Standard_Real aU = theU, aV = theV;
-    if (! shiftPoint (aU, aV, myBaseSurf, myBaseAdaptor))
+    gp_Vec aD2U, aD2V, aD2UV;
+    BaseD2 (aU, aV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
+    try 
     {
-      throw;
+      CalculateD1(aU, aV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
+      break;
+    } 
+    catch (Geom_UndefinedValue&)
+    {
+      // if failed at parametric boundary, try taking derivative at shifted point
+      if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V))
+      {
+        throw;
+      }
     }
-    BaseD2 (aU, aV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
-    CalculateD1(theU, theV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
   }
 }
 
@@ -186,26 +289,26 @@ void GeomEvaluator_OffsetSurface::D2(
     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV) const
 {
-  gp_Vec aD3U, aD3V, aD3UUV, aD3UVV;
-  BaseD3(theU, theV, theValue, theD1U, theD1V,
-         theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
-  try
+  Standard_Real aU = theU, aV = theV;
+  for (;;)
   {
-    CalculateD2(theU, theV, theValue, theD1U, theD1V,
-                theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
-  } 
-  catch (Geom_UndefinedValue&)
-  {
-    // if failed at parametric boundary, try taking derivative at shifted point
-    Standard_Real aU = theU, aV = theV;
-    if (! shiftPoint (aU, aV, myBaseSurf, myBaseAdaptor))
+    gp_Vec aD3U, aD3V, aD3UUV, aD3UVV;
+    BaseD3 (aU, aV, theValue, theD1U, theD1V,
+            theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
+    try
+    {
+      CalculateD2 (aU, aV, theValue, theD1U, theD1V,
+                   theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
+      break;
+    } 
+    catch (Geom_UndefinedValue&)
     {
-      throw;
+      // if failed at parametric boundary, try taking derivative at shifted point
+      if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V))
+      {
+        throw;
+      }
     }
-    BaseD3(theU, theV, theValue, theD1U, theD1V,
-         theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
-    CalculateD2(theU, theV, theValue, theD1U, theD1V,
-                theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
   }
 }
 
@@ -215,25 +318,25 @@ void GeomEvaluator_OffsetSurface::D3(
     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
 {
-  BaseD3(theU, theV, theValue, theD1U, theD1V,
-         theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
-  try
+  Standard_Real aU = theU, aV = theV;
+  for (;;)
   {
-    CalculateD3(theU, theV, theValue, theD1U, theD1V,
-                theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
-  } 
-  catch (Geom_UndefinedValue&)
-  {
-    // if failed at parametric boundary, try taking derivative at shifted point
-    Standard_Real aU = theU, aV = theV;
-    if (! shiftPoint (aU, aV, myBaseSurf, myBaseAdaptor))
+    BaseD3 (aU, aV, theValue, theD1U, theD1V,
+            theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
+    try
+    {
+      CalculateD3 (aU, aV, theValue, theD1U, theD1V,
+                   theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
+      break;
+    } 
+    catch (Geom_UndefinedValue&)
     {
-      throw;
+      // if failed at parametric boundary, try taking derivative at shifted point
+      if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V))
+      {
+        throw;
+      }
     }
-    BaseD3(aU, aV, theValue, theD1U, theD1V,
-         theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
-    CalculateD3(theU, theV, theValue, theD1U, theD1V,
-                theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
   }
 }
 
@@ -246,23 +349,24 @@ gp_Vec GeomEvaluator_OffsetSurface::DN(
   Standard_RangeError_Raise_if(theDerU + theDerV < 1,
       "GeomEvaluator_OffsetSurface::DN(): theDerU + theDerV < 1");
 
-  gp_Pnt aP;
-  gp_Vec aD1U, aD1V;
-  BaseD1(theU, theV, aP, aD1U, aD1V);
-  try
-  {
-    return CalculateDN(theU, theV, theDerU, theDerV, aD1U, aD1V);
-  } 
-  catch (Geom_UndefinedValue&)
+  Standard_Real aU = theU, aV = theV;
+  for (;;)
   {
-    // if failed at parametric boundary, try taking derivative at shifted point
-    Standard_Real aU = theU, aV = theV;
-    if (! shiftPoint (aU, aV, myBaseSurf, myBaseAdaptor))
+    gp_Pnt aP;
+    gp_Vec aD1U, aD1V;
+    BaseD1 (aU, aV, aP, aD1U, aD1V);
+    try
+    {
+      return CalculateDN (aU, aV, theDerU, theDerV, aD1U, aD1V);
+    } 
+    catch (Geom_UndefinedValue&)
     {
-      throw;
+      // if failed at parametric boundary, try taking derivative at shifted point
+      if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, aD1U, aD1V))
+      {
+        throw;
+      }
     }
-    BaseD1 (aU, aV, aP, aD1U, aD1V);
-    return CalculateDN (theU, theV, theDerU, theDerV, aD1U, aD1V);
   }
 }
 
@@ -315,8 +419,6 @@ void GeomEvaluator_OffsetSurface::CalculateD0(
     gp_Pnt& theValue,
     const gp_Vec& theD1U, const gp_Vec& theD1V) const
 {
-  const Standard_Real MagTol = 1.e-9;
-
   // Normalize derivatives before normal calculation because it gives more stable result.
   // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
   gp_Vec aD1U(theD1U);
@@ -329,7 +431,7 @@ void GeomEvaluator_OffsetSurface::CalculateD0(
     aD1V /= Sqrt(aD1VNorm2);
 
   gp_Vec aNorm = aD1U.Crossed(aD1V);
-  if (aNorm.SquareMagnitude() > MagTol * MagTol)
+  if (aNorm.SquareMagnitude() > the_D1MagTol * the_D1MagTol)
   {
     // Non singular case. Simple computations.
     aNorm.Normalize();
@@ -365,15 +467,15 @@ void GeomEvaluator_OffsetSurface::CalculateD0(
 
     gp_Dir Normal;
     CSLib_NormalStatus NStatus;
-    CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
+    CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
                   NStatus, Normal, OrderU, OrderV);
     if (NStatus == CSLib_InfinityOfSolutions)
     {
       // Replace zero derivative and try to calculate normal
       gp_Vec aNewDU = theD1U;
       gp_Vec aNewDV = theD1V;
-      if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, MagTol * MagTol))
-        CSLib::Normal(aNewDU, aNewDV, MagTol, NStatus, Normal);
+      if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, the_D1MagTol * the_D1MagTol))
+        CSLib::Normal(aNewDU, aNewDV, the_D1MagTol, NStatus, Normal);
     }
 
     if (NStatus != CSLib_Defined)
@@ -389,8 +491,6 @@ void GeomEvaluator_OffsetSurface::CalculateD1(
     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
     const gp_Vec& theD2U, const gp_Vec& theD2V, const gp_Vec& theD2UV) const
 {
-  const Standard_Real MagTol = 1.e-9;
-
   // Check offset side.
   Handle(Geom_BSplineSurface) L;
   Standard_Boolean isOpposite = Standard_False;
@@ -411,7 +511,7 @@ void GeomEvaluator_OffsetSurface::CalculateD1(
   Standard_Boolean isSingular = Standard_False;
   Standard_Integer MaxOrder = 0;
   gp_Vec aNorm = aD1U.Crossed(aD1V);
-  if (aNorm.SquareMagnitude() <= MagTol * MagTol)
+  if (aNorm.SquareMagnitude() <= the_D1MagTol * the_D1MagTol)
   {
     if (!myOscSurf.IsNull())
     {
@@ -478,13 +578,13 @@ void GeomEvaluator_OffsetSurface::CalculateD1(
 
   gp_Dir Normal;
   CSLib_NormalStatus NStatus;
-  CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax, NStatus, Normal, OrderU, OrderV);
+  CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax, NStatus, Normal, OrderU, OrderV);
   if (NStatus == CSLib_InfinityOfSolutions)
   {
     gp_Vec aNewDU = theD1U;
     gp_Vec aNewDV = theD1V;
     // Replace zero derivative and try to calculate normal
-    if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, MagTol * MagTol))
+    if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, the_D1MagTol * the_D1MagTol))
     {
       DerSurf.SetValue(1, 0, aNewDU);
       DerSurf.SetValue(0, 1, aNewDV);
@@ -492,7 +592,7 @@ void GeomEvaluator_OffsetSurface::CalculateD1(
         derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
       else
         derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
-      CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
+      CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
                     NStatus, Normal, OrderU, OrderV);
     }
   }
@@ -513,11 +613,9 @@ void GeomEvaluator_OffsetSurface::CalculateD2(
     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
     const gp_Vec& theD3U, const gp_Vec& theD3V, const gp_Vec& theD3UUV, const gp_Vec& theD3UVV) const
 {
-  const Standard_Real MagTol = 1.e-9;
-
   gp_Dir Normal;
   CSLib_NormalStatus NStatus;
-  CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
+  CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal);
 
   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
   Standard_Integer OrderU, OrderV;
@@ -554,7 +652,7 @@ void GeomEvaluator_OffsetSurface::CalculateD2(
   else
     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
 
-  CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
+  CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
                 NStatus, Normal, OrderU, OrderV);
   if (NStatus != CSLib_Defined)
     throw Geom_UndefinedValue(
@@ -589,11 +687,9 @@ void GeomEvaluator_OffsetSurface::CalculateD3(
     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
 {
-  const Standard_Real MagTol = 1.e-9;
-
   gp_Dir Normal;
   CSLib_NormalStatus NStatus;
-  CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
+  CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal);
   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
   Standard_Integer OrderU, OrderV;
   TColgp_Array2OfVec DerNUV(0, MaxOrder + 3, 0, MaxOrder + 3);
@@ -629,7 +725,7 @@ void GeomEvaluator_OffsetSurface::CalculateD3(
   else
     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
 
-  CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
+  CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
                 NStatus, Normal, OrderU, OrderV);
   if (NStatus != CSLib_Defined)
     throw Geom_UndefinedValue(
@@ -675,11 +771,9 @@ gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(
     const Standard_Integer theNu, const Standard_Integer theNv,
     const gp_Vec& theD1U, const gp_Vec& theD1V) const
 {
-  const Standard_Real MagTol = 1.e-9;
-
   gp_Dir Normal;
   CSLib_NormalStatus NStatus;
-  CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
+  CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal);
   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
   Standard_Integer OrderU, OrderV;
   TColgp_Array2OfVec DerNUV(0, MaxOrder + theNu, 0, MaxOrder + theNu);
@@ -709,7 +803,7 @@ gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(
   else
     derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
 
-  CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
+  CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
                 NStatus, Normal, OrderU, OrderV);
   if (NStatus != CSLib_Defined)
     throw Geom_UndefinedValue(
@@ -811,117 +905,3 @@ Standard_Boolean GeomEvaluator_OffsetSurface::ReplaceDerivative(
   }
   return isReplaced;
 }
-
-
-
-
-
-// =====================   Auxiliary functions   ===================================
-template<class SurfOrAdapt>
-void derivatives(Standard_Integer theMaxOrder,
-                 Standard_Integer theMinOrder,
-                 const Standard_Real theU,
-                 const Standard_Real theV,
-                 const SurfOrAdapt& theBasisSurf,
-                 const Standard_Integer theNU,
-                 const Standard_Integer theNV,
-                 const Standard_Boolean theAlongU,
-                 const Standard_Boolean theAlongV,
-                 const Handle(Geom_BSplineSurface)& theL,
-                 TColgp_Array2OfVec& theDerNUV,
-                 TColgp_Array2OfVec& theDerSurf)
-{
-  Standard_Integer i, j;
-  gp_Pnt P;
-  gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V;
-
-  if (theAlongU || theAlongV)
-  {
-    theMaxOrder = 0;
-    TColgp_Array2OfVec DerSurfL(0, theMaxOrder + theNU + 1, 0, theMaxOrder + theNV + 1);
-    switch (theMinOrder)
-    {
-    case 1:
-      theL->D1(theU, theV, P, DL1U, DL1V);
-      DerSurfL.SetValue(1, 0, DL1U);
-      DerSurfL.SetValue(0, 1, DL1V);
-      break;
-    case 2:
-      theL->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV);
-      DerSurfL.SetValue(1, 0, DL1U);
-      DerSurfL.SetValue(0, 1, DL1V);
-      DerSurfL.SetValue(1, 1, DL2UV);
-      DerSurfL.SetValue(2, 0, DL2U);
-      DerSurfL.SetValue(0, 2, DL2V);
-      break;
-    case 3:
-      theL->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV);
-      DerSurfL.SetValue(1, 0, DL1U);
-      DerSurfL.SetValue(0, 1, DL1V);
-      DerSurfL.SetValue(1, 1, DL2UV);
-      DerSurfL.SetValue(2, 0, DL2U);
-      DerSurfL.SetValue(0, 2, DL2V);
-      DerSurfL.SetValue(3, 0, DL3U);
-      DerSurfL.SetValue(2, 1, DL3UUV);
-      DerSurfL.SetValue(1, 2, DL3UVV);
-      DerSurfL.SetValue(0, 3, DL3V);
-      break;
-    default:
-      break;
-    }
-
-    if (theNU <= theNV)
-    {
-      for (i = 0; i <= theMaxOrder + 1 + theNU; i++)
-        for (j = i; j <= theMaxOrder + theNV + 1; j++)
-          if (i + j > theMinOrder)
-          {
-            DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
-            theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
-            if (i != j && j <= theNU + 1)
-            {
-              theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
-              DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
-            }
-          }
-    }
-    else
-    {
-      for (j = 0; j <= theMaxOrder + 1 + theNV; j++)
-        for (i = j; i <= theMaxOrder + theNU + 1; i++)
-          if (i + j > theMinOrder)
-          {
-            DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
-            theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
-            if (i != j && i <= theNV + 1)
-            {
-              theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
-              DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
-            }
-          }
-    }
-    for (i = 0; i <= theMaxOrder + theNU; i++)
-      for (j = 0; j <= theMaxOrder + theNV; j++)
-      {
-        if (theAlongU)
-          theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf));
-        if (theAlongV)
-          theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL));
-      }
-  }
-  else
-  {
-    for (i = 0; i <= theMaxOrder + theNU+ 1; i++)
-      for (j = i; j <= theMaxOrder + theNV + 1; j++)
-        if (i + j > theMinOrder)
-        {
-          theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
-          if (i != j)
-            theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
-        }
-    for (i = 0; i <= theMaxOrder + theNU; i++)
-      for (j = 0; j <= theMaxOrder + theNV; j++)
-        theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf));
-  }
-}
-

diff --git a/src/QABugs/QABugs_20.cxx b/src/QABugs/QABugs_20.cxx
index e6a6a960d55df8f3f8e9db23d4858f5cdd309a6a..4b7a1160d0bd42290bd3700ff6f980554a92e2e9 100644 (file)
--- a/src/QABugs/QABugs_20.cxx
+++ b/src/QABugs/QABugs_20.cxx
@@ -24,9 +24,12 @@
 #include <TColStd_Array2OfReal.hxx>
 #include <TColStd_Array1OfReal.hxx>
 #include <TColStd_Array1OfInteger.hxx>
+#include <Geom_BezierCurve.hxx>
 #include <Geom_BSplineSurface.hxx>
+#include <GeomConvert.hxx>
 #include <Geom2d_Curve.hxx>
 #include <Geom2d_Line.hxx>
+#include <GeomFill_BSplineCurves.hxx>
 #include <Draw.hxx>
 #include <DrawTrSurf.hxx>
 #include <ShapeConstruct_ProjectCurveOnSurface.hxx>
@@ -2435,6 +2438,99 @@ static Standard_Integer OCC28887 (Draw_Interpretor&, Standard_Integer theNbArgs,
   return 0;
 }
 
+static Standard_Integer OCC28131 (Draw_Interpretor&, Standard_Integer theNbArgs, const char** theArgVec)
+{
+  if (theNbArgs != 2)
+  {
+    std::cerr << "Error: wrong number of arguments" << std::endl;
+    return 1;
+  }
+
+  double height = 8.5;
+  gp_Pnt JiZhunXian2_v0 = gp_Pnt(-17.6, 0.0, 0.0);
+  gp_Pnt JiZhunXian2_v1 = gp_Pnt(0, 32.8, 0.0);
+
+  // Outline
+  TColgp_Array1OfPnt outer_e_bzr_geom_v(1, 4);
+  {
+    outer_e_bzr_geom_v(1) = JiZhunXian2_v0;
+    outer_e_bzr_geom_v(4) = JiZhunXian2_v1;
+
+    Standard_Real ratio1 = 5.4 / 13.2;
+    outer_e_bzr_geom_v(2) = gp_Pnt(outer_e_bzr_geom_v(1).X(), ratio1*outer_e_bzr_geom_v(4).Y(), 0);
+    Standard_Real ratio2 = 6.0 / 6.8;
+    outer_e_bzr_geom_v(3) = gp_Pnt(ratio2*outer_e_bzr_geom_v(1).X(), outer_e_bzr_geom_v(4).Y(), 0);
+  }
+
+  Handle(Geom_BezierCurve) outer_e_bzr_geom = new Geom_BezierCurve(outer_e_bzr_geom_v);
+  Handle(Geom_BSplineCurve) outer_e_bsp_geom = GeomConvert::CurveToBSplineCurve(outer_e_bzr_geom);
+  TopoDS_Edge outer_e = BRepBuilderAPI_MakeEdge(outer_e_bsp_geom);
+
+  Handle(Geom_BSplineCurve) curve1;
+  {
+    Handle(TColgp_HArray1OfPnt2d) harray = new TColgp_HArray1OfPnt2d(1, 2); // sizing harray
+    harray->SetValue(1, gp_Pnt2d(-JiZhunXian2_v1.Y(), 0));
+    harray->SetValue(2, gp_Pnt2d(0, height + height / 2));
+
+    Geom2dAPI_Interpolate anInterpolation(harray, Standard_False, 1e-6);
+
+    gp_Vec2d vtangent1(0, 1);
+    gp_Vec2d vtangent2(1, 0);
+    anInterpolation.Load(vtangent1, vtangent2);
+    anInterpolation.Perform();
+
+    Handle(Geom2d_BSplineCurve) c = anInterpolation.Curve();
+
+    gp_Pln pln(gp_Ax3(gp_Pnt(), gp_Dir(1, 0, 0), gp_Dir(0, -1, 0)));
+
+    Handle(Geom_BSplineCurve) c3d = Handle(Geom_BSplineCurve)::DownCast(GeomAPI::To3d(c, pln));
+    curve1 = c3d;
+  }
+
+  Handle(Geom_BSplineCurve) curve2;
+  {
+    Handle(TColgp_HArray1OfPnt2d) harray = new TColgp_HArray1OfPnt2d(1, 3); // sizing harray
+    harray->SetValue(1, gp_Pnt2d(-JiZhunXian2_v0.X(), 0));
+    harray->SetValue(2, gp_Pnt2d(-JiZhunXian2_v0.X() - 2.6, height));
+    harray->SetValue(3, gp_Pnt2d(0, height + height / 2));
+
+    Geom2dAPI_Interpolate anInterpolation(harray, Standard_False, 1e-6);
+    anInterpolation.Perform();
+
+    Handle(Geom2d_BSplineCurve) c = anInterpolation.Curve();
+    gp_Pln pln(gp_Ax3(gp_Pnt(), gp_Dir(0, -1, 0), gp_Dir(-1, 0, 0)));
+    Handle(Geom_BSplineCurve) c3d = Handle(Geom_BSplineCurve)::DownCast(GeomAPI::To3d(c, pln));
+    curve2 = c3d;
+  }
+
+  //////////////////////////////////////
+  GeomFill_BSplineCurves fill2;
+  fill2.Init(outer_e_bsp_geom, curve1, curve2, GeomFill_CoonsStyle);
+
+  const Handle(Geom_BSplineSurface)& surf_geom = fill2.Surface();
+
+  TopoDS_Shape filled_face = BRepBuilderAPI_MakeFace(surf_geom, 0);
+
+  DBRep::Set (theArgVec[1], filled_face);
+
+/*
+  ///////////////////////////////////////////////////////////////////////
+  TopoDS_Solid first_solid;
+  {
+    BRepOffset_MakeOffset myOffsetShape(filled_face, -offset_thick, 1e-4,
+      BRepOffset_Skin, //Mode
+      Standard_False, //Intersection
+      Standard_False, //SelfInter
+      GeomAbs_Intersection, //Join
+      Standard_True, //Thickening
+      Standard_False //RemoveIntEdges
+      ); //RemoveInvalidFaces
+    first_solid = TopoDS::Solid(myOffsetShape.Shape());
+  }
+*/
+  return 0;
+}
+
 void QABugs::Commands_20(Draw_Interpretor& theCommands) {
   const char *group = "QABugs";
 
@@ -2462,6 +2558,7 @@ void QABugs::Commands_20(Draw_Interpretor& theCommands) {
                   "OCC28887 filePath result"
                   "\n\t\t: Check interface for reading BRep from memory.",
                   __FILE__, OCC28887, group);
+  theCommands.Add("OCC28131", "OCC28131 name: creates face problematic for offset", __FILE__, OCC28131, group);
 
   return;
 }

diff --git a/tests/bugs/modalg_7/bug25730 b/tests/bugs/modalg_7/bug25730
index 18ad471f844bde1d78e3c93dfa3c0b9c7f0d150c..82b990debd41d5f3cdfe94250ed2580102c7fab8 100755 (executable)
--- a/tests/bugs/modalg_7/bug25730
+++ b/tests/bugs/modalg_7/bug25730
@@ -1,27 +1,25 @@
-puts "TODO OCC25730 ALL: result of MakeThickSolid aborts the BOPCheck in Geom_OffsetSurface::SetD0"
-
-puts "============"
-puts "OCC25730"
-puts "============"
+puts "# ===================================================================="
+puts "# OCC25730: result of MakeThickSolid aborts the BOPCheck in Geom_OffsetSurface::SetD0"
+puts "# ===================================================================="
 puts ""
-#############################################################################################
-## result of MakeThickSolid aborts the BOPCheck in Geom_OffsetSurface::SetD0
-#############################################################################################
 
+puts "# Load shape"
 restore [locate_data_file bug25730_thickness8-draw-fillet001.brep] Fillet001
 
 explode Fillet001 F
 
+puts "# Perform offset"
 offsetparameter 1e-7 p a
 offsetload Fillet001 -1 Fillet001_4
 offsetperform Thickness
 
-if { [regexp "There were errors during the operation, so the list may be incomplete" [bopcheck Fillet001]] == 1 } {
-    puts "Error : result of MakeThickSolid aborts the BOPCheck in Geom_OffsetSurface::SetD0"
+puts "# Check result"
+if { [regexp "errors" [bopcheck Fillet001]] == 1 } {
+    puts "Error : bopcheck fails on initial shape"
 }
 
-if { [regexp "There were errors during the operation, so the list may be incomplete" [bopcheck Thickness]] == 1 } {
-    puts "Error : result of MakeThickSolid aborts the BOPCheck in Geom_OffsetSurface::SetD0"
+if { [regexp "errors" [bopcheck Thickness]] == 1 } {
+    puts "Error : bopcheck fails on offsetted shape"
 }
 
 checkview -display Fillet001 -2d -path ${imagedir}/${test_image}-Fillet001-2d.png

diff --git a/tests/bugs/modalg_7/bug28131 b/tests/bugs/modalg_7/bug28131
new file mode 100644 (file)
index 0000000..8af534c
--- /dev/null
+++ b/tests/bugs/modalg_7/bug28131
@@ -0,0 +1,38 @@
+puts "# ==============================================================="
+puts "# 0028131: BRepOffset_MakeOffset can't create offset with a face which created by filling 3 bsplinecurve"
+puts "# ==============================================================="
+puts ""
+
+puts "# Create face to be offset, by dedicated command"
+pload QAcommands
+OCC28131 f
+
+puts "# Try simple offset"
+offsetshapesimple result_simple f 10.
+checkshape result_simple
+checkmaxtol result_simple -ref 0.205
+checkprops result_simple -s 1693.7
+
+puts "# Try standard offset"
+offsetshape result_std f 10.
+fixshape result_std result_std ;# need to fix it....
+checkshape result_std
+checkmaxtol result_std -ref 0.408
+checkprops result_std -s 1693.76
+
+puts "# Make snapshots (overall and zoom to degenerated point)"
+
+smallview -Y+Z
+fit
+checkview -2d -screenshot -path ${imagedir}/${test_image}.png
+
+smallview -Y+Z
+zoom 400
+pu; pu; pu
+pr; pr; pr
+
+donly result_simple
+checkview -2d -screenshot -path ${imagedir}/${test_image}_zoom_simple.png
+
+donly result_std
+checkview -2d -screenshot -path ${imagedir}/${test_image}_zoom_standard.png