From: jfa Date: Wed, 8 Jun 2022 07:24:48 +0000 (+0300) Subject: [bos #29472] [EDF] (2022-T1) Advanced geometry features: curvature vector on a point... X-Git-Tag: V9_9_1b1~1^2 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=335b028279b21ca634bba57228daf66895d591de;p=modules%2Fgeom.git [bos #29472] [EDF] (2022-T1) Advanced geometry features: curvature vector on a point of a face --- 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