diff --git a/doc/salome/examples/curvature_face.py b/doc/salome/examples/curvature_face.py
new file mode 100644
index 000000000..3c22b34c0
--- /dev/null
+++ b/doc/salome/examples/curvature_face.py
@@ -0,0 +1,203 @@
+# Curvature of a Face along given direction
+
+import salome
+salome.salome_init_without_session()
+import GEOM
+from salome.geom import geomBuilder
+geompy = geomBuilder.New()
+import math
+import numpy as np
+
+def test_acceptance():
+ """
+ Acceptance test [tuleap29472]
+ """
+ Vector = [0,100,100]
+ O = geompy.MakeVertex(0, 0, 0)
+ OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
+ OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
+ OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
+ Cylinder_1 = geompy.MakeCylinderRH(100, 300)
+ Translation_1 = geompy.MakeTranslation(Cylinder_1, 0, 0, -150)
+ Vertex_1 = geompy.MakeVertex(100, 0, 0)
+ Vertex_2 = geompy.MakeVertex(100, -Vector[2], Vector[1])
+ Line_1 = geompy.MakeLineTwoPnt(Vertex_1, Vertex_2)
+ Plane_1 = geompy.MakePlane(Vertex_1, Line_1, 2000)
+ Rotation_1 = geompy.MakeRotation(Translation_1, OZ, 90*math.pi/180.0)# avoid to have degenerated edge across Vertex_1
+
+ [Face_1,Face_2,Face_3] = geompy.ExtractShapes(Rotation_1, geompy.ShapeType["FACE"], True)
+
+ curvature_29472 = np.array( geompy.VectorCoordinates( geompy.CurvatureOnFace(Face_2, Vertex_1, geompy.MakeVectorDXDYDZ(*Vector))) ).reshape(1,3)
+ expected_curvature = np.array( [-200.0,0.0,0.0] ).reshape(1,3)
+ assert( np.isclose( 0.0, np.linalg.norm( curvature_29472 - expected_curvature ) ,rtol=0,atol=1e-5 ) )
+
+ Intersection_1 = geompy.MakeSection(Face_2, Plane_1, True)
+ geompy.addToStudy( O, 'O' )
+ geompy.addToStudy( OX, 'OX' )
+ geompy.addToStudy( OY, 'OY' )
+ geompy.addToStudy( OZ, 'OZ' )
+ geompy.addToStudy( Vertex_1, 'Vertex_1' )
+ geompy.addToStudy( Cylinder_1, 'Cylinder_1' )
+ geompy.addToStudy( Translation_1, 'Translation_1' )
+ geompy.addToStudy( Vertex_2, 'Vertex_2' )
+ geompy.addToStudy( Line_1, 'Line_1' )
+ geompy.addToStudy( Plane_1, 'Plane_1' )
+ geompy.addToStudy( Rotation_1, 'Rotation_1' )
+ geompy.addToStudyInFather( Rotation_1, Face_1, 'Face_1' )
+ geompy.addToStudyInFather( Rotation_1, Face_2, 'Face_2' )
+ geompy.addToStudyInFather( Rotation_1, Face_3, 'Face_3' )
+ geompy.addToStudy( Intersection_1, 'Intersection_1' )
+ angle = math.asin(Vector[2]/math.sqrt(Vector[1]*Vector[1]+Vector[2]*Vector[2]))
+ tmp = geompy.MakeTranslation(Intersection_1,*[-elt for elt in geompy.PointCoordinates(Vertex_1)])
+ tmp = geompy.MakeRotation(tmp,OX,-angle)
+ Intersection_1_OXY = geompy.MakeTranslation(tmp,*geompy.PointCoordinates(Vertex_1))
+ geompy.addToStudy( Intersection_1_OXY, 'Intersection_1_OXY' )
+
+ eps = 0.01
+ offset = 0.75
+ p0 = np.array( geompy.PointCoordinates( geompy.MakeVertexOnCurve(Intersection_1_OXY,offset-eps) ) ).reshape(1,3)
+ p1 = np.array( geompy.PointCoordinates( geompy.MakeVertexOnCurve(Intersection_1_OXY,offset) ) ).reshape(1,3)
+ p2 = np.array( geompy.PointCoordinates( geompy.MakeVertexOnCurve(Intersection_1_OXY,offset+eps) ) ).reshape(1,3)
+ assert( np.isclose(0.0,np.linalg.norm(p1- np.array(geompy.PointCoordinates(Vertex_1)).reshape(1,3) ),rtol=0,atol=1e-8) )
+ p01=(p0+p1)/2
+ p12=(p1+p2)/2
+ v0 = (p1-p0)/np.linalg.norm(p1-p0)
+ v1 = (p2-p1)/np.linalg.norm(p2-p1)
+ computedRadius = 1/np.linalg.norm((v1-v0)/np.linalg.norm(p12-p01))
+ # manual detection of radius : https://fr.wikipedia.org/wiki/Courbure_d%27un_arc
+ circle = geompy.MakeCircle(O,OZ,computedRadius)
+ circle = geompy.MakeTranslation(circle,100-computedRadius,0,0)
+ geompy.addToStudy(circle, "expectedCircle")
+ print("Radius expected is {}".format(computedRadius))
+ print("Radius obtain by CurvatureOnFace is {}".format(np.linalg.norm(curvature_29472)))
+
+O = geompy.MakeVertex(0, 0, 0, 'O')
+OX = geompy.MakeVectorDXDYDZ(1, 0, 0, 'OX')
+OY = geompy.MakeVectorDXDYDZ(0, 1, 0, 'OY')
+OZ = geompy.MakeVectorDXDYDZ(0, 0, 1, 'OZ')
+
+pXYZ = geompy.MakeVertex(105, 105, 105, 'pXYZ')
+pY = geompy.MakeVertex(0, 105, 0, 'pY')
+pZ = geompy.MakeVertex(0, 0, 105, 'pZ')
+
+vZ_XY = geompy.MakeVectorDXDYDZ(-1, -1, 1, 'vZ-XY')
+vZ_XY2 = geompy.MakeVectorDXDYDZ(-1, -1, 10, 'vZ-XY')
+vZ_XY3 = geompy.MakeVectorDXDYDZ(-1, -1, 100, 'vZ-XY')
+
+R = 100.0
+
+# I. Curvature of a Sphere
+Sphere_1 = geompy.MakeSphereR(R, 'Sphere_1')
+[Sph] = geompy.ExtractShapes(Sphere_1, geompy.ShapeType["FACE"], True, "Sph")
+
+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')
+
+# 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)
+
+# Pole
+isExcept = False
+try:
+ geompy.CurvatureOnFace(Sph, pZ, OX)
+except:
+ isExcept = True
+assert(isExcept)
+
+# Normal direction
+isExcept = False
+try:
+ geompy.CurvatureOnFace(Sph, pY, OY)
+except:
+ isExcept = True
+assert(isExcept)
+
+# II. Curvature of a Cylinder
+Cylinder_1 = geompy.MakeCylinderRH(R, 300, 'Cylinder_1')
+[Face_1,Face_2,Face_3] = geompy.ExtractShapes(Cylinder_1, geompy.ShapeType["FACE"], True, "Face")
+
+# Curvature radius of a cylinder along any direction, orthogonal to its Z axis, equal to R
+curvature_4 = geompy.CurvatureOnFace(Face_2, pY, OX, 'curvature_cyl_pY_OX')
+assert(abs(geompy.BasicProperties(curvature_4)[0] - R) < 1e-07)
+
+# Curvature radius of a cylinder along its Z direction is infinite
+curvature_zero = geompy.CurvatureOnFace(Face_2, pY, OZ)
+assert(geompy.MeasuOp.GetErrorCode() == "ZERO_CURVATURE")
+assert(not curvature_zero)
+
+# Curvature radius of a cylinder along some direction, different from two above
+curvature_5 = geompy.CurvatureOnFace(Face_2, pY, vZ_XY, 'curvature_cyl_pY_vZ_XY')
+curvature_6 = geompy.CurvatureOnFace(Face_2, pY, vZ_XY2, 'curvature_cyl_pY_vZ_XY2')
+curvature_7 = geompy.CurvatureOnFace(Face_2, pY, vZ_XY3, 'curvature_cyl_pY_vZ_XY3')
+
+# R < r5 < r6 < r7
+# r5 = 100.01, r6 = 101.0, r7 = 200
+r5 = geompy.BasicProperties(curvature_5)[0]
+r6 = geompy.BasicProperties(curvature_6)[0]
+r7 = geompy.BasicProperties(curvature_7)[0]
+
+assert(R + 1e-07 < r5)
+assert(r5 + 1e-07 < r6)
+assert(r6 + 1e-07 < r7)
+
+# Projection aborted. Point is out of the face boundaries.
+isExcept = False
+try:
+ pXY_Z = geompy.MakeVertex(105, 105, -105, 'pXY_Z')
+ geompy.CurvatureOnFace(Face_2, pXY_Z, OX, 'curvature_cyl_pXY_Z')
+except:
+ isExcept = True
+assert(isExcept)
+
+# Projection aborted (point on axis). Equal distances to many points.
+isExcept = False
+try:
+ geompy.CurvatureOnFace(Face_2, O, vZ_XY, 'curvature_cyl_O')
+except:
+ isExcept = True
+assert(isExcept)
+
+# Curvature radius of a planar face is infinite
+curvature_zero_2 = geompy.CurvatureOnFace(Face_1, pZ, OX)
+assert(geompy.MeasuOp.GetErrorCode() == "ZERO_CURVATURE")
+assert(not curvature_zero_2)
+
+# III. Curvature of a "Horse saddle"
+[Edge_1,Edge_2,Edge_3] = geompy.ExtractShapes(Sphere_1, geompy.ShapeType["EDGE"], True)
+geompy.addToStudyInFather( Sphere_1, Edge_1, 'Edge_1' )
+geompy.addToStudyInFather( Sphere_1, Edge_2, 'Edge_2' )
+geompy.addToStudyInFather( Sphere_1, Edge_3, 'Edge_3' )
+
+Rotation_1 = geompy.MakeRotation(Edge_3, OX, 90*math.pi/180.0, 'Rotation_1')
+Rotation_2 = geompy.MakeRotation(Rotation_1, OY, 180*math.pi/180.0, 'Rotation_2')
+Translation_1 = geompy.MakeTranslation(Rotation_2, 200, 0, 0, 'Translation_1')
+Translation_2 = geompy.MakeTranslation(Edge_3, 100, 100, 0, 'Translation_2')
+Translation_3 = geompy.MakeTranslation(Translation_2, 0, -200, 0, 'Translation_3')
+Filling_1 = geompy.MakeFilling([Translation_2, Edge_3, Translation_3])
+geompy.addToStudy(Filling_1, 'Filling_1')
+Vertex_2 = geompy.MakeVertex(100, 0, 0, 'Vertex_2')
+
+curvature_Y = geompy.CurvatureOnFace(Filling_1, Vertex_2, OY, 'curvature_Y')
+curvature_Z = geompy.CurvatureOnFace(Filling_1, Vertex_2, OZ, 'curvature_Z')
+
+cury = np.array( geompy.VectorCoordinates(curvature_Y) ).reshape(1,3)
+curz = np.array( geompy.VectorCoordinates(curvature_Z) ).reshape(1,3)
+cury_expected = np.array( [50,0,0] ).reshape(1,3)
+curz_expected = np.array( [-100,0,0] ).reshape(1,3)
+assert( np.isclose( 0.0, np.linalg.norm( cury - cury_expected ) ,rtol=0,atol=1e-5 ) )
+assert( np.isclose( 0.0, np.linalg.norm( curz - curz_expected ) ,rtol=0,atol=1e-5 ) )
+
+# Normal direction
+norm_1 = geompy.GetNormal(Filling_1, Vertex_2, "Normal_1")
+isExcept = False
+try:
+ geompy.CurvatureOnFace(Filling_1, Vertex_2, norm_1)
+except:
+ isExcept = True
+assert(isExcept)
+
+# acceptance case
+test_acceptance()
\ No newline at end of file
diff --git a/doc/salome/examples/tests.set b/doc/salome/examples/tests.set
index 47d9c909c..327d64070 100644
--- a/doc/salome/examples/tests.set
+++ b/doc/salome/examples/tests.set
@@ -77,6 +77,7 @@ SET(GOOD_TESTS
import_export.py
inertia.py
min_distance.py
+ curvature_face.py
normal_face.py
notebook_geom.py
polyline.py
diff --git a/doc/salome/gui/GEOM/input/tui_test_all.doc b/doc/salome/gui/GEOM/input/tui_test_all.doc
index 7c2ddd136..af06ea0fa 100644
--- a/doc/salome/gui/GEOM/input/tui_test_all.doc
+++ b/doc/salome/gui/GEOM/input/tui_test_all.doc
@@ -108,6 +108,9 @@
\until geompy.GetSubShapesWithTolerance(Box, GEOM.FACE, GEOM.CC_LE, 1.e-7, "le")
\anchor swig_MakeExtraction
-\until print "DONE"
+\until geompy.MakeExtraction(Box, [16], "Ext_no_vertex")
+
+\anchor swig_CurvatureOnFace
+\until print("DONE")
*/
diff --git a/idl/GEOM_Gen.idl b/idl/GEOM_Gen.idl
index 6af219e02..36c8059b2 100644
--- a/idl/GEOM_Gen.idl
+++ b/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:
diff --git a/src/GEOMGUI/GEOM_msg_en.ts b/src/GEOMGUI/GEOM_msg_en.ts
index abc4146e7..ef6042244 100644
--- a/src/GEOMGUI/GEOM_msg_en.ts
+++ b/src/GEOMGUI/GEOM_msg_en.ts
@@ -4980,6 +4980,10 @@ Please, select face, shell or solid and try again
STB_NORMALE
Compute normal to the face
+
+ MEN_CURVATURE
+ Vector of curvature
+
TOP_MEASURE_ANGLE
Angle
diff --git a/src/GEOMGUI/GEOM_msg_fr.ts b/src/GEOMGUI/GEOM_msg_fr.ts
index f47f8e987..b58f0550b 100644
--- a/src/GEOMGUI/GEOM_msg_fr.ts
+++ b/src/GEOMGUI/GEOM_msg_fr.ts
@@ -4972,6 +4972,10 @@ Choisissez une face, une coque ou un solide et essayez de nouveau
STB_NORMALE
Vecteur normal à une face
+
+ MEN_CURVATURE
+ Vecteur de courbure
+
TOP_MEASURE_ANGLE
Angle
diff --git a/src/GEOMGUI/GEOM_msg_ja.ts b/src/GEOMGUI/GEOM_msg_ja.ts
index c6708f292..6fbd8f6ef 100644
--- a/src/GEOMGUI/GEOM_msg_ja.ts
+++ b/src/GEOMGUI/GEOM_msg_ja.ts
@@ -4975,6 +4975,10 @@
STB_NORMALE
フェースに垂直
+
+ MEN_CURVATURE
+ Vector_of_curvature
+
TOP_MEASURE_ANGLE
角度
diff --git a/src/GEOMImpl/GEOMImpl_IMeasure.hxx b/src/GEOMImpl/GEOMImpl_IMeasure.hxx
index 1e3489440..fa94ffc46 100644
--- a/src/GEOMImpl/GEOMImpl_IMeasure.hxx
+++ b/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;
diff --git a/src/GEOMImpl/GEOMImpl_IMeasureOperations.cxx b/src/GEOMImpl/GEOMImpl_IMeasureOperations.cxx
index 69c6e409e..ff539a6c9 100644
--- a/src/GEOMImpl/GEOMImpl_IMeasureOperations.cxx
+++ b/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.
diff --git a/src/GEOMImpl/GEOMImpl_IMeasureOperations.hxx b/src/GEOMImpl/GEOMImpl_IMeasureOperations.hxx
index 58b27bbda..2d2002613 100644
--- a/src/GEOMImpl/GEOMImpl_IMeasureOperations.hxx
+++ b/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
diff --git a/src/GEOMImpl/GEOMImpl_MeasureDriver.cxx b/src/GEOMImpl/GEOMImpl_MeasureDriver.cxx
index 9d430e4c7..122a4e9a5 100644
--- a/src/GEOMImpl/GEOMImpl_MeasureDriver.cxx
+++ b/src/GEOMImpl/GEOMImpl_MeasureDriver.cxx
@@ -33,10 +33,12 @@
#include
#include
#include
+#include
#include
#include
#include
#include
+#include
#include
#include
@@ -49,7 +51,9 @@
#include
#include
+#include
#include
+#include
#include
#include
@@ -80,6 +84,91 @@ 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)
+{
+ // Point
+ if (thePoint.IsNull())
+ Standard_NullObject::Raise("Point for curvature measurement is null");
+ if (thePoint.ShapeType() != TopAbs_VERTEX)
+ Standard_TypeMismatch::Raise("Point for curvature calculation is not a vertex");
+ gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(thePoint));
+
+ // Point projection on the face
+ Standard_Real U, V;
+ aPnt = GEOMUtils::ProjectPointOnFace(aPnt, theFace, U, V);
+ gp_Pnt2d UV (U, V);
+
+ // Face and Vector
+ TopoDS_Face aFace = TopoDS::Face(theFace);
+ gp_Vec aV = GEOMUtils::GetVector(theDir, Standard_False);
+
+ // Calculate differential properties
+ BRepAdaptor_Surface aSurfAdapt (aFace);
+ BRepLProp_SLProps Props (aSurfAdapt, UV.X(), UV.Y(), 2, 1e-7);
+ if (!Props.IsNormalDefined())
+ Standard_ConstructionError::Raise
+ ("Curvature calculation failed: normal direction is not defined");
+
+ // Get differential properties
+ gp_Vec n = Props.Normal();
+ if (aV.IsParallel(n, Precision::Angular()))
+ Standard_ConstructionError::Raise
+ ("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 {
+ 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);
+ BRepBuilderAPI_MakeEdge aBuilder (aPnt, aPntEnd);
+ if (!aBuilder.IsDone())
+ Standard_ConstructionError::Raise("Curvature calculation failed: edge is not built");
+
+ return aBuilder.Shape();
+}
+
//=======================================================================
//function : Execute
//purpose :
@@ -310,6 +399,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 {
}
@@ -359,6 +459,12 @@ GetCreationInformation(std::string& theOperationName,
AddParam( theParams, "Face", aCI.GetBase() );
AddParam( theParams, "Point", aCI.GetPoint(), "face center" );
break;
+ case CURVATURE_VEC_MEASURE:
+ theOperationName = "CURVATURE";
+ AddParam( theParams, "Face", aCI.GetBase() );
+ AddParam( theParams, "Point", aCI.GetPoint(), "point of interest" );
+ AddParam( theParams, "Vector", aCI.GetDirection(), "direction of interest" );
+ break;
default:
return false;
}
diff --git a/src/GEOMImpl/GEOMImpl_ProjectionDriver.cxx b/src/GEOMImpl/GEOMImpl_ProjectionDriver.cxx
index d131b0702..0929093d9 100644
--- a/src/GEOMImpl/GEOMImpl_ProjectionDriver.cxx
+++ b/src/GEOMImpl/GEOMImpl_ProjectionDriver.cxx
@@ -127,70 +127,11 @@ Standard_Integer GEOMImpl_ProjectionDriver::Execute(Handle(TFunction_Logbook)& l
Handle(GEOM_Function) aTargetFunction = TI.GetPlane();
if (aTargetFunction.IsNull()) return 0;
TopoDS_Shape aFaceShape = aTargetFunction->GetValue();
- //if (aFaceShape.IsNull() || aFaceShape.ShapeType() != TopAbs_FACE) {
- // Standard_ConstructionError::Raise
- // ("Projection aborted : the target shape is not a face");
- //}
-
- Standard_Real tol = 1.e-4;
if (anOriginal.ShapeType() == TopAbs_VERTEX) {
- if (aFaceShape.IsNull() || aFaceShape.ShapeType() != TopAbs_FACE) {
- Standard_ConstructionError::Raise
- ("Projection aborted : the target shape is not a face");
- }
- TopoDS_Face aFace = TopoDS::Face(aFaceShape);
- Handle(Geom_Surface) surface = BRep_Tool::Surface(aFace);
- double U1, U2, V1, V2;
- //surface->Bounds(U1, U2, V1, V2);
- BRepTools::UVBounds(aFace, U1, U2, V1, V2);
-
- // projector
- GeomAPI_ProjectPointOnSurf proj;
- proj.Init(surface, U1, U2, V1, V2, tol);
-
- gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(anOriginal));
- proj.Perform(aPnt);
- if (!proj.IsDone()) {
- Standard_ConstructionError::Raise
- ("Projection aborted : the algorithm failed");
- }
- int nbPoints = proj.NbPoints();
- if (nbPoints < 1) {
- Standard_ConstructionError::Raise("No solution found");
- }
-
Standard_Real U, V;
- proj.LowerDistanceParameters(U, V);
- gp_Pnt2d aProjPnt (U, V);
-
- // classifier
- BRepClass_FaceClassifier aClsf (aFace, aProjPnt, tol);
- if (aClsf.State() != TopAbs_IN && aClsf.State() != TopAbs_ON) {
- bool isSol = false;
- double minDist = RealLast();
- for (int i = 1; i <= nbPoints; i++) {
- Standard_Real Ui, Vi;
- proj.Parameters(i, Ui, Vi);
- aProjPnt = gp_Pnt2d(Ui, Vi);
- aClsf.Perform(aFace, aProjPnt, tol);
- if (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON) {
- isSol = true;
- double dist = proj.Distance(i);
- if (dist < minDist) {
- minDist = dist;
- U = Ui;
- V = Vi;
- }
- }
- }
- if (!isSol) {
- Standard_ConstructionError::Raise("No solution found");
- }
- }
-
- gp_Pnt surfPnt = surface->Value(U, V);
-
+ gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(anOriginal));
+ gp_Pnt surfPnt = GEOMUtils::ProjectPointOnFace(aPnt, aFaceShape, U, V);
aShape = BRepBuilderAPI_MakeVertex(surfPnt).Shape();
}
else {
diff --git a/src/GEOMImpl/GEOMImpl_Types.hxx b/src/GEOMImpl/GEOMImpl_Types.hxx
index 1f64191fa..faa8f9953 100644
--- a/src/GEOMImpl/GEOMImpl_Types.hxx
+++ b/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
diff --git a/src/GEOMUtils/GEOMUtils.cxx b/src/GEOMUtils/GEOMUtils.cxx
index 3bd842996..4c69a9a05 100644
--- a/src/GEOMUtils/GEOMUtils.cxx
+++ b/src/GEOMUtils/GEOMUtils.cxx
@@ -37,6 +37,7 @@
#include
#include
+#include
#include
#include
@@ -64,6 +65,8 @@
#include
#include
+#include
+
#include
#include
#include
@@ -1021,6 +1024,65 @@ Standard_Real GEOMUtils::GetMinDistance
return aResult;
}
+//=======================================================================
+// function : ProjectPointOnFace()
+// purpose : Returns the projection (3d point) if found, throws an exception otherwise
+//=======================================================================
+gp_Pnt GEOMUtils::ProjectPointOnFace(const gp_Pnt& thePoint,
+ const TopoDS_Shape& theFace,
+ double& theU, double& theV)
+{
+ if (theFace.IsNull() || theFace.ShapeType() != TopAbs_FACE)
+ Standard_TypeMismatch::Raise
+ ("Projection aborted : the target shape is not a face");
+
+ TopoDS_Face aFace = TopoDS::Face(theFace);
+ Handle(Geom_Surface) surface = BRep_Tool::Surface(aFace);
+ double U1, U2, V1, V2;
+ BRepTools::UVBounds(aFace, U1, U2, V1, V2);
+
+ // projector
+ Standard_Real tol = 1.e-4;
+ GeomAPI_ProjectPointOnSurf proj;
+ proj.Init(surface, U1, U2, V1, V2, tol);
+ proj.Perform(thePoint);
+ if (!proj.IsDone())
+ StdFail_NotDone::Raise("Projection aborted : the algorithm failed");
+ int nbPoints = proj.NbPoints();
+ if (nbPoints < 1)
+ Standard_ConstructionError::Raise("Projection aborted : No solution found");
+ proj.LowerDistanceParameters(theU, theV);
+ gp_Pnt2d aProjPnt (theU, theV);
+
+ // classifier
+ BRepClass_FaceClassifier aClsf (aFace, aProjPnt, tol);
+ if (aClsf.State() != TopAbs_IN && aClsf.State() != TopAbs_ON) {
+ bool isSol = false;
+ double minDist = RealLast();
+ for (int i = 1; i <= nbPoints; i++) {
+ Standard_Real Ui, Vi;
+ proj.Parameters(i, Ui, Vi);
+ aProjPnt = gp_Pnt2d(Ui, Vi);
+ aClsf.Perform(aFace, aProjPnt, tol);
+ if (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON) {
+ isSol = true;
+ double dist = proj.Distance(i);
+ if (dist < minDist) {
+ minDist = dist;
+ theU = Ui;
+ theV = Vi;
+ }
+ }
+ }
+ if (!isSol) {
+ Standard_ConstructionError::Raise("Projection aborted : No solution found");
+ }
+ }
+
+ gp_Pnt surfPnt = surface->Value(theU, theV);
+ return surfPnt;
+}
+
//=======================================================================
// function : ConvertClickToPoint()
// purpose : Returns the point clicked in 3D view
diff --git a/src/GEOMUtils/GEOMUtils.hxx b/src/GEOMUtils/GEOMUtils.hxx
index f41b21650..5546211fa 100644
--- a/src/GEOMUtils/GEOMUtils.hxx
+++ b/src/GEOMUtils/GEOMUtils.hxx
@@ -215,6 +215,19 @@ namespace GEOMUtils
const TopoDS_Shape& theShape2,
gp_Pnt& thePnt1, gp_Pnt& thePnt2);
+ /*!
+ * \brief Computes normal projection of \a thePoint to \a theFace.
+ *
+ * \param thePoint the 3d point
+ * \param theFace the face shape
+ * \param theU the output U parameter of the point on the face
+ * \param theV the output V parameter of the point on the face
+ * \retval the projection (3d point) if found, throws an exception otherwise
+ */
+ Standard_EXPORT gp_Pnt ProjectPointOnFace(const gp_Pnt& thePoint,
+ const TopoDS_Shape& theFace,
+ double& theU, double& theV);
+
/*!
* \brief Returns the point clicked in 3D view.
*
diff --git a/src/GEOM_I/GEOM_IMeasureOperations_i.cc b/src/GEOM_I/GEOM_IMeasureOperations_i.cc
index c2acb2158..ccc759c3e 100644
--- a/src/GEOM_I/GEOM_IMeasureOperations_i.cc
+++ b/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);
+}
diff --git a/src/GEOM_I/GEOM_IMeasureOperations_i.hh b/src/GEOM_I/GEOM_IMeasureOperations_i.hh
index bebaedeec..30b89faeb 100644
--- a/src/GEOM_I/GEOM_IMeasureOperations_i.hh
+++ b/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(); }
};
diff --git a/src/GEOM_SWIG/GEOM_TestAll.py b/src/GEOM_SWIG/GEOM_TestAll.py
index b3dfabb72..a182836db 100644
--- a/src/GEOM_SWIG/GEOM_TestAll.py
+++ b/src/GEOM_SWIG/GEOM_TestAll.py
@@ -595,5 +595,20 @@ def TestAll (geompy, math):
geompy.MakeExtraction(Box, [18], "Ext_no_edge")
geompy.MakeExtraction(Box, [16], "Ext_no_vertex")
+ # CurvatureOnFace
+ Cylinder_1 = geompy.MakeCylinderRH(100, 50, 'Cylinder_r100_h150')
+ [Face_1,Face_2,Face_3] = geompy.ExtractShapes(Cylinder_1, geompy.ShapeType["FACE"], True, "Face")
+ curvature_1 = geompy.CurvatureOnFace(Face_2, px, vy, 'curvature_cyl_px_vy')
+ assert(abs(geompy.BasicProperties(curvature_1)[0] - 100) < 1e-07)
+ curvature_zero = geompy.CurvatureOnFace(Face_2, px, vz)
+ assert(geompy.MeasuOp.GetErrorCode() == "ZERO_CURVATURE")
+ assert(not curvature_zero)
+ isExcept = False
+ try:
+ # p0 is on cylinder axis, projection should fail
+ geompy.CurvatureOnFace(Face_2, p0, vy)
+ except:
+ isExcept = True
+ assert(isExcept)
print("DONE")
diff --git a/src/GEOM_SWIG/geomBuilder.py b/src/GEOM_SWIG/geomBuilder.py
index e7d0364b6..d10e9118a 100644
--- a/src/GEOM_SWIG/geomBuilder.py
+++ b/src/GEOM_SWIG/geomBuilder.py
@@ -11178,6 +11178,56 @@ class geomBuilder(GEOM._objref_GEOM_Gen):
return aSurf
## @}
+ ## Measure curvature radius 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 radius.
+ # Null shape is returned in case of infinite radius
+ # (zero curvature), for example, in case of flat face.
+ #
+ ## @ref swig_CurvatureOnFace "Example"
+ @ManageTransactions("MeasuOp")
+ def CurvatureOnFace(self, theSurf, thePoint, theDirection, theName=None):
+ """
+ Measure curvature radius 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 radius.
+ Null shape is returned in case of infinite radius
+ (zero curvature), for example, in case of flat face.
+
+ Example of usage:
+ curvature_1 = geompy.CurvatureOnFace(Face_1, Vertex_1, OX)
+ """
+ aVec = self.MeasuOp.SurfaceCurvatureByPointAndDirection(theSurf,thePoint,theDirection)
+ if self.MeasuOp.GetErrorCode() != "ZERO_CURVATURE":
+ 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