From: jfa <jfa@opencascade.com>
Date: Wed, 3 Aug 2022 21:54:59 +0000 (+0300)
Subject: Fix a problem with wrong curvature values
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=55e7a4baa50daf34a491e22c34a7c406fdd3af7c;p=modules%2Fgeom.git

Fix a problem with wrong curvature values
---

diff --git a/doc/salome/examples/curvature_face.py b/doc/salome/examples/curvature_face.py
index 327a70750..262f51b3e 100644
--- a/doc/salome/examples/curvature_face.py
+++ b/doc/salome/examples/curvature_face.py
@@ -29,14 +29,19 @@ Sphere_1 = geompy.MakeSphereR(R, 'Sphere_1')
 curvature_1 = geompy.CurvatureOnFace(Sph, pXYZ, OX,    'curvature_sph_pXYZ_OX')
 curvature_2 = geompy.CurvatureOnFace(Sph, pXYZ, vZ_XY, 'curvature_sph_pXYZ_vt')
 curvature_3 = geompy.CurvatureOnFace(Sph, pY,   OX,    'curvature_sph_pY_OX')
-# Pole
-curvature_p = geompy.CurvatureOnFace(Sph, pZ,   OX,    'curvature_sph_pZ_OX')
 
 # All sphere curvature radiuces = R
 assert(abs(geompy.BasicProperties(curvature_1)[0] - R) < 1e-07)
 assert(abs(geompy.BasicProperties(curvature_2)[0] - R) < 1e-07)
 assert(abs(geompy.BasicProperties(curvature_3)[0] - R) < 1e-07)
-assert(abs(geompy.BasicProperties(curvature_p)[0] - R) < 1e-07)
+
+# Pole
+isExcept = False
+try:
+  geompy.CurvatureOnFace(Sph, pZ, OX)
+except:
+  isExcept = True
+assert(isExcept)
 
 # Normal direction
 isExcept = False
diff --git a/src/GEOMImpl/GEOMImpl_MeasureDriver.cxx b/src/GEOMImpl/GEOMImpl_MeasureDriver.cxx
index 4a5780999..122a4e9a5 100644
--- a/src/GEOMImpl/GEOMImpl_MeasureDriver.cxx
+++ b/src/GEOMImpl/GEOMImpl_MeasureDriver.cxx
@@ -118,49 +118,47 @@ TopoDS_Shape EvaluateAlongCurvature(const TopoDS_Shape& theFace,
       ("Curvature calculation failed: normal direction is not defined");
 
   // Get differential properties
-  gp_Vec Xu  = Props.D1U();
-  gp_Vec Xv  = Props.D1V();
-  gp_Vec Xuu = Props.D2U();
-  gp_Vec Xuv = Props.DUV();
-  gp_Vec Xvv = Props.D2V();
-  gp_Vec n   = Props.Normal();
-
-  // Direction in 2d
-  gp_Dir aDirU (Xu);
-  gp_Dir aDirV (Xv);
-  gp_Vec2d T (aV.Dot(aDirU), aV.Dot(aDirV));
-  if (Abs(T.X()) < Precision::Confusion() &&
-      Abs(T.Y()) < Precision::Confusion())
+  gp_Vec n = Props.Normal();
+  if (aV.IsParallel(n, Precision::Angular()))
     Standard_ConstructionError::Raise
-                   ("Curvature calculation failed: direction is normal to the face");
-
-  // Coefficients of the FFF
-  double E = Xu.Dot(Xu);
-  double F = Xu.Dot(Xv);
-  double G = Xv.Dot(Xv);
-
-  // Coefficients of the SFF
-  double L = n.Dot(Xuu);
-  double M = n.Dot(Xuv);
-  double N = n.Dot(Xvv);
-
-  // Calculate radius (or -radius) of curvature
-  // using the coefficients of both fundamental forms
-  double r = 0.;
-  if (Abs(T.X()) < Precision::Confusion()) {
-    if (Abs(N) < Precision::Confusion())
-      Standard_Failure::Raise("ZERO_CURVATURE");
-    r = G / N;
+      ("Curvature calculation failed: direction is normal to the face");
+
+  if (!Props.IsCurvatureDefined())
+    Standard_ConstructionError::Raise
+      ("Curvature calculation failed: curvature is not defined");
+
+  // Curvature
+  Standard_Real k = 0.;
+
+  if (Props.IsUmbilic()) {
+    k = Props.MaxCurvature();
   }
   else {
-    double lambda = T.Y() / T.X();
-    double detE = E + 2*F*lambda + G*lambda*lambda;
-    double detL = L + 2*M*lambda + N*lambda*lambda;
-    if (Abs(detL) < Precision::Confusion())
-      Standard_Failure::Raise("ZERO_CURVATURE");
-    r = detE / detL;
+    gp_Dir aMaxDir, aMinDir;
+    Props.CurvatureDirections(aMaxDir, aMinDir);
+
+    gp_Vec2d T (aV.Dot(aMinDir), aV.Dot(aMaxDir));
+    if (Abs(T.X()) < Precision::Confusion() &&
+        Abs(T.Y()) < Precision::Confusion())
+      Standard_ConstructionError::Raise
+                     ("Curvature calculation failed: direction is normal to the face");
+    gp_Dir2d aDirT (T);
+    Standard_Real cosAlpha = aDirT.X();
+
+    Standard_Real aMinCurv = Props.MinCurvature();
+    Standard_Real aMaxCurv = Props.MaxCurvature();
+
+    // Euler formula: k = k1 * cos^2(alpha) + k2 * sin^2(alpha)
+    //                k = k2 + (k1 - k2) * cos^2(alpha)
+    k = aMaxCurv + (aMinCurv - aMaxCurv) * cosAlpha * cosAlpha;
   }
 
+  if (Abs(k) < Precision::Confusion())
+    Standard_Failure::Raise("ZERO_CURVATURE");
+
+  // Radius or -radius of curvature
+  Standard_Real r = 1.0 / k;
+
   // Result
   gp_Dir aNormal (n);
   gp_Pnt aPntEnd (aPnt.XYZ() + aNormal.XYZ() * r);