[bos #29472] [EDF] (2022-T1) Advanced geometry features: curvature vector on a point of a face

idl/GEOM_Gen.idl

@@ -4716,6 +4716,24 @@ module GEOM
      */
     double MinSurfaceCurvatureByPoint (in GEOM_Object theShape, in GEOM_Object thePoint);
 
+    /*!
+     *  \brief Get vector of curvature of surface in the given point along the given direction.
+     *  \param theShape - face.
+     *  \param thePoint - point.
+     *  \param theDirection - direction.
+     *  \note  Before the calculation of curvature, the point and the direction
+     *         are projected to the face, if the point does not lay on it or
+     *         the direction is not tangent to it initially.
+     *  \return Vector of curvature. The returned vector is codirectional with
+     *          the normal to the face in the given point in case of positive
+     *          curvature value and opposite to the normal in case of negative
+     *          curvature. The normal of the returned vector is equal to the
+     *          absolute value of the curvature.
+     */
+    GEOM_Object SurfaceCurvatureByPointAndDirection (in GEOM_Object theShape,
+                                                     in GEOM_Object thePoint,
+                                                     in GEOM_Object theDirection);
+
   };
 
  // # GEOM_IGroupOperations:

src/GEOMImpl/GEOMImpl_IMeasure.hxx

@@ -30,7 +30,8 @@ class GEOMImpl_IMeasure
     MEASURE_ARG_BASE  = 1,
     MEASURE_ARG_POINT = 2,
     MEASURE_INDEX = 3,
-    MEASURE_USE_ORI = 4
+    MEASURE_USE_ORI = 4,
+    MEASURE_ARG_DIR = 5
   };
  public:
 
@@ -56,6 +57,11 @@ class GEOMImpl_IMeasure
                               !_func->IsDone() ); // old behavior was to useOri
   }
 
+  void SetDirection(Handle(GEOM_Function) theDir)
+  { _func->SetReference(MEASURE_ARG_DIR, theDir); }
+
+  Handle(GEOM_Function) GetDirection() { return _func->GetReference(MEASURE_ARG_DIR); }
+
  private:
 
   Handle(GEOM_Function) _func;

src/GEOMImpl/GEOMImpl_IMeasureOperations.cxx

@@ -2658,6 +2658,63 @@ Standard_Real GEOMImpl_IMeasureOperations::MinSurfaceCurvatureByPoint
   return getSurfaceCurvatures(aSurf, UV.X(), UV.Y(), false);
 }
 
+//=============================================================================
+/*!
+ *  SurfaceCurvatureByPointAndDirection
+ */
+//=============================================================================
+Handle(GEOM_Object) GEOMImpl_IMeasureOperations::SurfaceCurvatureByPointAndDirection
+                                                 (Handle(GEOM_Object) theSurf,
+                                                  Handle(GEOM_Object) thePoint,
+                                                  Handle(GEOM_Object) theDirection)
+{
+  SetErrorCode(KO);
+
+  if (theSurf.IsNull() || thePoint.IsNull() || theDirection.IsNull()) return NULL;
+
+  Handle(GEOM_Function) aSurf = theSurf->GetLastFunction();
+  Handle(GEOM_Function) aPoint = thePoint->GetLastFunction();
+  Handle(GEOM_Function) aDirection = theDirection->GetLastFunction();
+  if (aSurf.IsNull() || aPoint.IsNull() || aDirection.IsNull()) return NULL;
+
+  //Add a new CurvatureVector object
+  //Handle(GEOM_Object) aCV = GetEngine()->AddObject(GEOM_CURVATURE_VEC);
+  Handle(GEOM_Object) aCV = GetEngine()->AddObject(GEOM_VECTOR);
+
+  //Add a new CurvatureVector function
+  Handle(GEOM_Function) aFunction =
+    aCV->AddFunction(GEOMImpl_MeasureDriver::GetID(), CURVATURE_VEC_MEASURE);
+  if (aFunction.IsNull()) return NULL;
+
+  //Check if the function is set correctly
+  if (aFunction->GetDriverGUID() != GEOMImpl_MeasureDriver::GetID()) return NULL;
+
+  GEOMImpl_IMeasure aCI (aFunction);
+  aCI.SetBase(aSurf);
+  aCI.SetPoint(aPoint);
+  aCI.SetDirection(aDirection);
+
+  //Compute the CurvatureVector
+  try {
+    OCC_CATCH_SIGNALS;
+    if (!GetSolver()->ComputeFunction(aFunction)) {
+      SetErrorCode("Measure driver failed to compute a surface curvature");
+      return NULL;
+    }
+  }
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
+    return NULL;
+  }
+
+  //Make a Python command
+  GEOM::TPythonDump(aFunction) << aCV << " = geompy.CurvatureOnFace(" << theSurf
+                               << ", " << thePoint << ", " << theDirection << ")";
+
+  SetErrorCode(OK);
+  return aCV;
+}
+
 //=======================================================================
 //function : FillErrorsSub
 //purpose  : Fill the errors list of subshapes on shape.

src/GEOMImpl/GEOMImpl_IMeasureOperations.hxx

@@ -216,6 +216,11 @@ class GEOMImpl_IMeasureOperations : public GEOM_IOperations {
   Standard_EXPORT Standard_Real MinSurfaceCurvatureByPoint (Handle(GEOM_Object) theSurf,
                                                             Handle(GEOM_Object) thePoint);
 
+  Standard_EXPORT Handle(GEOM_Object) SurfaceCurvatureByPointAndDirection
+                                      (Handle(GEOM_Object) theSurf,
+                                       Handle(GEOM_Object) thePoint,
+                                       Handle(GEOM_Object) theDirection);
+
  private:
 
    void FillErrorsSub

src/GEOMImpl/GEOMImpl_MeasureDriver.cxx

@@ -33,6 +33,7 @@
 #include <BRep_Tool.hxx>
 #include <BRepTools.hxx>
 #include <BRepGProp.hxx>
+#include <BRepLProp_SLProps.hxx>
 #include <BRepBuilderAPI_Copy.hxx>
 #include <BRepBuilderAPI_MakeEdge.hxx>
 #include <BRepBuilderAPI_MakeVertex.hxx>
@@ -80,6 +81,101 @@ GEOMImpl_MeasureDriver::GEOMImpl_MeasureDriver()
 {
 }
 
+//! This function is designed to evaluate normal curvature of the surface
+//! in the given point along the given direction.
+//! param[in] theFace face of interest.
+//! param[in] thePoint point of interest.
+//! param[in] theDir edge, giving the direction of interest.
+//! return Edge, representing the curvature vector
+TopoDS_Shape EvaluateAlongCurvature(const TopoDS_Shape& theFace,
+                                    const TopoDS_Shape& thePoint,
+                                    const TopoDS_Shape& theDir)
+{
+  if (theFace.IsNull())
+    Standard_NullObject::Raise("Face for curvature calculation is null");
+  if (theFace.ShapeType() != TopAbs_FACE)
+    Standard_NullObject::Raise("Shape for curvature calculation is not a face");
+  TopoDS_Face aFace = TopoDS::Face(theFace);
+  Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
+  if (aSurf.IsNull())
+    Standard_NullObject::Raise("Surface for curvature calculation is null");
+
+  if (thePoint.IsNull())
+    Standard_NullObject::Raise("Point for curvature measurement is null");
+  if (thePoint.ShapeType() != TopAbs_VERTEX)
+    Standard_NullObject::Raise("Point for curvature calculation is not a vertex");
+  gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(thePoint));
+
+  gp_Vec aV = GEOMUtils::GetVector(theDir, Standard_False);
+
+  // Point projection parameters on surface
+  ShapeAnalysis_Surface aSAS (aSurf);
+  gp_Pnt2d UV = aSAS.ValueOfUV(aPnt, Precision::Confusion());
+  aPnt = aSAS.Value(UV);
+
+  // Calculate differential properties
+  BRepAdaptor_Surface aSurfAdapt (aFace);
+  BRepLProp_SLProps Props (aSurfAdapt, UV.X(), UV.Y(), 2, 1e-7);
+  if (!Props.IsCurvatureDefined())
+    Standard_NullObject::Raise("Curvature calculation failed");
+
+  // 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())
+    Standard_NullObject::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 curvature using the coefficients of both fundamental forms
+  double k = 0.;
+  if (Abs(T.X()) < Precision::Confusion()) {
+    //if (Abs(G) < Precision::Confusion())
+    //  Standard_NullObject::Raise("Curvature calculation failed: G = 0");
+    //k = N / G; // curvature
+    if (Abs(N) < Precision::Confusion())
+      Standard_NullObject::Raise("Curvature calculation failed: N = 0");
+    k = G / N; // radius of curvature
+  }
+  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(detE) < Precision::Confusion())
+    if (Abs(detL) < Precision::Confusion())
+      Standard_NullObject::Raise("Curvature calculation failed: det = 0");
+    //k = detL / detE; // curvature
+    k = detE / detL; // radius of curvature
+  }
+
+  // Result
+  gp_Dir aNormal (n);
+  gp_Pnt aPntEnd (aPnt.XYZ() + aNormal.XYZ() * k);
+  BRepBuilderAPI_MakeEdge aBuilder (aPnt, aPntEnd);
+  if (!aBuilder.IsDone())
+    Standard_NullObject::Raise("Curvature calculation failed: edge is not built");
+
+  return aBuilder.Shape();
+}
+
 //=======================================================================
 //function : Execute
 //purpose  :
@@ -310,6 +406,17 @@ Standard_Integer GEOMImpl_MeasureDriver::Execute(Handle(TFunction_Logbook)& log)
       Standard_NullObject::Raise("Vector construction failed");
     aShape = aBuilder.Shape();
   }
+  else if (aType == CURVATURE_VEC_MEASURE) {
+    Handle(GEOM_Function) aSrfFunc = aCI.GetBase();
+    Handle(GEOM_Function) aPntFunc = aCI.GetPoint();
+    Handle(GEOM_Function) aDirFunc = aCI.GetDirection();
+
+    TopoDS_Shape aFace   = aSrfFunc->GetValue();
+    TopoDS_Shape aVertex = aPntFunc->GetValue();
+    TopoDS_Shape anEdge  = aDirFunc->GetValue();
+
+    aShape = EvaluateAlongCurvature(aFace, aVertex, anEdge);
+  }
   else {
   }
 

src/GEOMImpl/GEOMImpl_Types.hxx

@@ -121,6 +121,8 @@
 
 #define GEOM_EXTRACTION 58
 
+#define GEOM_CURVATURE_VEC 59
+
 //GEOM_Function types
 
 #define COPY_WITH_REF    1
@@ -357,6 +359,7 @@
 #define BND_BOX_MEASURE_PRECISE 3
 #define VECTOR_FACE_NORMALE 4
 #define VERTEX_BY_INDEX 5
+#define CURVATURE_VEC_MEASURE 6
 
 #define GROUP_FUNCTION 1
 

src/GEOM_I/GEOM_IMeasureOperations_i.cc

@@ -1188,3 +1188,32 @@ CORBA::Double GEOM_IMeasureOperations_i::MinSurfaceCurvatureByPoint
 
   return GetOperations()->MinSurfaceCurvatureByPoint(aShape,aPoint);
 }
+
+//=============================================================================
+/*!
+ *  SurfaceCurvatureByPointAndDirection
+ */
+//=============================================================================
+GEOM::GEOM_Object_ptr GEOM_IMeasureOperations_i::SurfaceCurvatureByPointAndDirection
+                                          (GEOM::GEOM_Object_ptr theSurf,
+                                           GEOM::GEOM_Object_ptr thePoint,
+                                           GEOM::GEOM_Object_ptr theDirection)
+{
+  GEOM::GEOM_Object_var aGEOMObject;
+
+  //Set a not done flag
+  GetOperations()->SetNotDone();
+
+  //Get the reference shape
+  Handle(::GEOM_Object) aShape = GetObjectImpl(theSurf);
+  Handle(::GEOM_Object) aPoint = GetObjectImpl(thePoint);
+  Handle(::GEOM_Object) aDirection = GetObjectImpl(theDirection);
+  if (aShape.IsNull() || aPoint.IsNull() || aDirection.IsNull()) return aGEOMObject._retn();
+
+  Handle(::GEOM_Object) anObject =
+    GetOperations()->SurfaceCurvatureByPointAndDirection(aShape,aPoint,aDirection);
+  if (!GetOperations()->IsDone() || anObject.IsNull())
+    return aGEOMObject._retn();
+
+  return GetObject(anObject);
+}

src/GEOM_I/GEOM_IMeasureOperations_i.hh

@@ -164,6 +164,10 @@ class GEOM_I_EXPORT GEOM_IMeasureOperations_i :
   CORBA::Double MinSurfaceCurvatureByPoint (GEOM::GEOM_Object_ptr theSurf,
                                             GEOM::GEOM_Object_ptr thePoint);
 
+  GEOM::GEOM_Object_ptr SurfaceCurvatureByPointAndDirection (GEOM::GEOM_Object_ptr theSurf,
+                                                             GEOM::GEOM_Object_ptr thePoint,
+                                                             GEOM::GEOM_Object_ptr theDirection);
+
   ::GEOMImpl_IMeasureOperations* GetOperations()
   { return (::GEOMImpl_IMeasureOperations*)GetImpl(); }
 };

src/GEOM_SWIG/geomBuilder.py

@@ -11172,6 +11172,51 @@ class geomBuilder(GEOM._objref_GEOM_Gen):
             return aSurf
         ## @}
 
+        ## Measure curvature of surface in the given point along the given direction.
+        #  @param theSurf the given face.
+        #  @param thePoint given point.
+        #  @param theDirection given direction.
+        #  @param theName Object name; when specified, this parameter is used
+        #         for result publication in the study. Otherwise, if automatic
+        #         publication is switched on, default value is used for result name.
+        #
+        #  @return New GEOM.GEOM_Object, containing vector of curvature of theSurf.
+        #          The returned vector is codirectional with the normal to the face
+        #          in the given point in case of positive curvature value
+        #          and opposite to the normal in case of negative curvature.
+        #          The normal of the returned vector is equal to the
+        #          absolute value of the curvature.
+        #
+        ## @ref swig_todo "Example"
+        @ManageTransactions("MeasuOp")
+        def CurvatureOnFace(self, theSurf, thePoint, theDirection, theName=None):
+            """
+            Measure curvature of surface in the given point along the given direction.
+
+            Parameters:
+                theSurf the given face.
+                thePoint given point.
+                theDirection given direction.
+                theName Object name; when specified, this parameter is used
+                        for result publication in the study. Otherwise, if automatic
+                        publication is switched on, default value is used for result name.
+
+            Returns:
+                New GEOM.GEOM_Object, containing vector of curvature of theSurf.
+                The returned vector is codirectional with the normal to the face
+                in the given point in case of positive curvature value
+                and opposite to the normal in case of negative curvature.
+                The normal of the returned vector is equal to the
+                absolute value of the curvature.
+
+            Example of usage:
+                curvature_1 = geompy.CurvatureOnFace(Face_1, Vertex_1, OX)
+            """
+            aVec = self.MeasuOp.SurfaceCurvatureByPointAndDirection(theSurf,thePoint,theDirection)
+            RaiseIfFailed("CurvatureOnFace", self.MeasuOp)
+            self._autoPublish(aVec, theName, "curvature")
+            return aVec
+
         ## Get min and max tolerances of sub-shapes of theShape
         #  @param theShape Shape, to get tolerances of.
         #  @return [FaceMin,FaceMax, EdgeMin,EdgeMax, VertMin,VertMax]\n