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-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=a54e77dbcd8159448aa99a0d8b0f9a1b39e44cf3;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..327a70750 --- /dev/null +++ b/doc/salome/examples/curvature_face.py @@ -0,0 +1,130 @@ +# 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 + +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') +# 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) + +# 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 = geompy.VectorCoordinates(curvature_Y) +curz = geompy.VectorCoordinates(curvature_Z) + +# Vectors should be opposite, scalar product should be negative +assert(cury[0]*curz[0] + cury[1]*curz[1] + cury[2]*curz[2] < -1e-07) + +# 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) 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 f280885b6..d0fa16858 100644 --- a/idl/GEOM_Gen.idl +++ b/idl/GEOM_Gen.idl @@ -4726,6 +4726,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 bbf4c4a8f..aa141c288 100644 --- a/src/GEOMImpl/GEOMImpl_IMeasureOperations.cxx +++ b/src/GEOMImpl/GEOMImpl_IMeasureOperations.cxx @@ -2730,6 +2730,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 bb178d3d1..827f695ff 100644 --- a/src/GEOMImpl/GEOMImpl_IMeasureOperations.hxx +++ b/src/GEOMImpl/GEOMImpl_IMeasureOperations.hxx @@ -217,6 +217,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..4a5780999 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,93 @@ 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 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_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; + } + 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; + } + + // 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 +401,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 +461,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 6ea33686a..59b284ec2 100644 --- a/src/GEOMImpl/GEOMImpl_Types.hxx +++ b/src/GEOMImpl/GEOMImpl_Types.hxx @@ -123,6 +123,8 @@ #define GEOM_PATCH_FACE 59 +#define GEOM_CURVATURE_VEC 60 + //GEOM_Function types #define COPY_WITH_REF 1 @@ -359,6 +361,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 4a6bace40..0c09e09c4 100644 --- a/src/GEOM_I/GEOM_IMeasureOperations_i.cc +++ b/src/GEOM_I/GEOM_IMeasureOperations_i.cc @@ -1214,3 +1214,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 fef4b7d6e..6a5aa39ee 100644 --- a/src/GEOM_I/GEOM_IMeasureOperations_i.hh +++ b/src/GEOM_I/GEOM_IMeasureOperations_i.hh @@ -166,6 +166,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 5c9fc66b2..5ad918d86 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