eryar@163.com
Abstract. 本文介绍了参数曲面的第一基本公式,并应用曲面的第一基本公式,结合OpenCASCADE中计算多重积分的类,对任意参数曲面的面积进行计算。
Key Words. Parametric Curve, Parametric Surface, Gauss Integration, Global Properties
我们知道一元函数y=f(x)的图像是一条曲线,二元函数z=f(x,y)的图像是一张曲面。但是,把曲线曲面表示成参数方程则更加便利于研究,这种表示方法首先是由欧洲瑞士数学家Euler引进的。例如,在空间中的一条曲线可以表示为三个一元函数:
X=x(t), Y=y(t), Z=z(t)
在向量的概念出现后,空间中的一条曲线可以自然地表示为一个一元向量函数:
r=r(t)=(x(t), y(t), z(t))
用向量函数来表示曲线和曲面后,使曲线曲面一些量的计算方式比较统一。如曲线可以表示为一元向量函数,曲面可以表示为二元向量函数。
本文结合OpenCASCADE来介绍参数曲线曲面积分分别计算曲线弧长和曲面的面积。结合《微分几何》来更好地理解曲线曲面相关知识。
设曲线C的参数方程是r=r(t),命:
则s是该曲线的一个不变量,即它与空间中的坐标系的选择无关,也与该曲线的参数变换无关。前者是因为在笛卡尔直角坐标变换下,切向量的长度|r’(t)|是不变的,故s不变。关于后者可以通过积分的变量的替换来证明,设参数变换是:
并且
因此
根据积分的变量替换公式有:
不变量s的几何意义就是曲线段的弧长。这说明曲线参数t可以是任意的,但选择不同的参数得到的参数方程会有不同,但是曲线段的弧长是不变的。以曲线弧长作为曲线方程的参数,这样的方程称为曲线的自然参数方程Natural Parametric Equations。
由曲线的参数方程可知,曲线弧长的计算公式为:
几何意义就是在每个微元处的切向量的长度求和。
曲线弧长的计算就是一元函数的积分。OpenCASCADE中是如何计算任意曲线弧长的呢?直接找到相关的源码列举如下:(在类CPnts_AbscissaPoint中)
- // auxiliary functions to compute the length of the derivative
- static Standard_Real f3d(const Standard_Real X, const Standard_Address C)
- {
- gp_Pnt P;
- gp_Vec V;
- ((Adaptor3d_Curve*)C)->D1(X,P,V);
- return V.Magnitude();
- }
- static Standard_Real f2d(const Standard_Real X, const Standard_Address C)
- {
- gp_Pnt2d P;
- gp_Vec2d V;
- ((Adaptor2d_Curve2d*)C)->D1(X,P,V);
- return V.Magnitude();
- }
- //==================================================================
- //function : Length
- //purpose : 3d with parameters
- //==================================================================
- Standard_Real CPnts_AbscissaPoint::Length(const Adaptor3d_Curve& C,
- const Standard_Real U1,
- const Standard_Real U2)
- {
- CPnts_MyGaussFunction FG;
- //POP pout WNT
- CPnts_RealFunction rf = f3d;
- FG.Init(rf,(Standard_Address)&C);
- // FG.Init(f3d,(Standard_Address)&C);
- math_GaussSingleIntegration TheLength(FG, U1, U2, order(C));
- if (!TheLength.IsDone()) {
- throw Standard_ConstructionError();
- }
- return Abs(TheLength.Value());
- }
- //==================================================================
- //function : Length
- //purpose : 2d with parameters
- //==================================================================
- Standard_Real CPnts_AbscissaPoint::Length(const Adaptor2d_Curve2d& C,
- const Standard_Real U1,
- const Standard_Real U2)
- {
- CPnts_MyGaussFunction FG;
- //POP pout WNT
- CPnts_RealFunction rf = f2d;
- FG.Init(rf,(Standard_Address)&C);
- // FG.Init(f2d,(Standard_Address)&C);
- math_GaussSingleIntegration TheLength(FG, U1, U2, order(C));
- if (!TheLength.IsDone()) {
- throw Standard_ConstructionError();
- }
- return Abs(TheLength.Value());
- }
上述代码的意思是直接对曲线的一阶导数的长度求积分,即是弧长。OpenCASCADE的代码写得有点难懂,根据意思把对三维曲线求弧长的代码改写下,更便于理解:
- //! Function for curve length evaluation.
- class math_LengthFunction: public math_Function
- {
- public:
- math_LengthFunction(const Handle(Geom_Curve) & theCurve)
- : myCurve(theCurve)
- {
- }
- virtual Standard_Boolean Value(const Standard_Real X, Standard_Real & F)
- {
- gp_Pnt aP;
- gp_Vec aV;
- myCurve - >D1(X, aP, aV);
- F = aV.Magnitude();
- return Standard_True;
- }
- private:
- Handle(Geom_Curve) myCurve;
- };
曲面参数方程是个二元向量函数。根据《微分几何》中曲面的第一基本公式(First Fundamental Form of a Surface)可知,曲面上曲线的表达式为:
r=r(u(t), v(t)) = (x(t), y(t), z(t))
若以s表示曲面上曲线的弧长,则由复合函数求导公式可得弧长微分公式:
在古典微分几何中,上式称为曲面的第一基本公式,E、F、G称为第一基本量。在曲面上,每一点的第一基本量与参数化无关。
利用曲面第一基本公式可以用于计算曲面的面积。参数曲面上与u,v参数空间的元素dudv对应的面积元为:
由参数曲面法向的计算可知,曲面的面积元素即为u,v方向上的偏导数的乘积的模。
其几何意义可以理解为参数曲面的面积微元是由u,v方向的偏导数的向量围成的一个四边形的面积,则整个曲面的面积即是对面积元素求积分。由于参数曲面有两个参数,所以若要计算曲面的面积,只需要对面积元素计算二重积分即可。
OpenCASCADE的math包中提供了多重积分的计算类math_GaussMultipleIntegration,由类名可知积分算法采用了Gauss积分算法。下面通过具体代码来说明OpenCASCADE中计算曲面积分的过程。要计算积分,先要定义被积函数。因为参数曲面与参数曲线不同,参数曲线只有一个参数,而参数曲面有两个参数,所以是一个多元函数。
- //! 2D variable function for surface area evaluation.
- class math_AreaFunction: public math_MultipleVarFunction
- {
- public:
- math_AreaFunction(const Handle(Geom_Surface) & theSurface)
- : mySurface(theSurface)
- {
- }
- virtual Standard_Integer NbVariables() const
- {
- return 2;
- }
- virtual Standard_Boolean Value(const math_Vector & X, Standard_Real & Y)
- {
- gp_Pnt aP;
- gp_Vec aDu;
- gp_Vec aDv;
- mySurface - >D1(X(1), X(2), aP, aDu, aDv);
- Y = aDu.Crossed(aDv).Magnitude();
- return Standard_True;
- }
- private:
- Handle(Geom_Surface) mySurface;
- };
由于参数曲面是多元函数,所以从类math_MultipleVarFunction派生,并在虚函数NbVariables()中说明有两个变量。在虚函数Value()中计算面积元素的值,即根据曲面第一基本公式中面积元素的定义,对参数曲面求一阶导数,计算两个偏导数向量的叉乘的模。
有了被积函数,只需要在定义域对其计算二重积分,相应代码如下所示:
- void evalArea(const Handle(Geom_Surface) & theSurface, const math_Vector & theLower, const math_Vector & theUpper)
- {
- math_IntegerVector aOrder(1, 2, math: :GaussPointsMax());
- math_AreaFunction aFunction(theSurface);
- math_GaussMultipleIntegration anIntegral(aFunction, theLower, theUpper, aOrder);
- if (anIntegral.IsDone())
- {
- anIntegral.Dump(std: :cout);
- }
- }
通过theLower和theUpper指定定义域,由于采用了Gauss-Legendre算法计算二重积分,所以需要指定阶数,且阶数越高积分结果精度越高,这里使用了OpenCASCADE中最高的阶数。
下面通过对基本曲面的面积计算来验证结果的正确性,并将计算结果和OpenCASCADE中计算面积的类BRepGProp::SurfaceProperties()结果进行对比。
下面通过对OpenCASCADE中几个初等曲面的面积进行计算,代码如下所示:
- /*
- Copyright(C) 2017 Shing Liu(eryar@163.com)
- Permission is hereby granted, free of charge, to any person obtaining a copy
- of this software and associated documentation files(the "Software"), to deal
- in the Software without restriction, including without limitation the rights
- to use, copy, modify, merge, publish, distribute, sublicense, and / or sell
- copies of the Software, and to permit persons to whom the Software is
- furnished to do so, subject to the following conditions :
- The above copyright notice and this permission notice shall be included in all
- copies or substantial portions of the Software.
- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
- IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
- FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
- AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
- LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
- OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
- SOFTWARE.
- */
- #include <gp_Pnt.hxx>
- #include <gp_Vec.hxx>
- #include <math.hxx>
- #include <math_Function.hxx>
- #include <math_MultipleVarFunction.hxx>
- #include <math_GaussMultipleIntegration.hxx>
- #include <Geom_Plane.hxx>
- #include <Geom_ConicalSurface.hxx>
- #include <Geom_CylindricalSurface.hxx>
- #include <Geom_SphericalSurface.hxx>
- #include <Geom_ToroidalSurface.hxx>
- #include <Geom_BSplineSurface.hxx>
- #include <Geom_RectangularTrimmedSurface.hxx>
- #include <GeomConvert.hxx>
- #include <GProp_GProps.hxx>
- #include <TopoDS_Face.hxx>
- #include <BRepGProp.hxx>
- #include <BRepBuilderAPI_MakeFace.hxx>
- #pragma comment(lib, "TKernel.lib")
- #pragma comment(lib, "TKMath.lib")
- #pragma comment(lib, "TKG2d.lib")
- #pragma comment(lib, "TKG3d.lib")
- #pragma comment(lib, "TKGeomBase.lib")
- #pragma comment(lib, "TKGeomAlgo.lib")
- #pragma comment(lib, "TKBRep.lib")
- #pragma comment(lib, "TKTopAlgo.lib")
- //! 2D variable function for surface area evaluation.
- class math_AreaFunction : public math_MultipleVarFunction
- {
- public:
- math_AreaFunction(const Handle(Geom_Surface)& theSurface)
- : mySurface(theSurface)
- {
- }
- virtual Standard_Integer NbVariables() const
- {
- return 2;
- }
- virtual Standard_Boolean Value(const math_Vector& X, Standard_Real& Y)
- {
- gp_Pnt aP;
- gp_Vec aDu;
- gp_Vec aDv;
- Standard_Real E = 0.0;
- Standard_Real F = 0.0;
- Standard_Real G = 0.0;
- mySurface->D1(X(1), X(2), aP, aDu, aDv);
- E = aDu.Dot(aDu);
- F = aDu.Dot(aDv);
- G = aDv.Dot(aDv);
- Y = Sqrt(E * G - F * F);
- //Y = aDu.Crossed(aDv).Magnitude();
- return Standard_True;
- }
- private:
- Handle(Geom_Surface) mySurface;
- };
- void evalArea(const Handle(Geom_Surface)& theSurface, const math_Vector& theLower, const math_Vector& theUpper)
- {
- math_IntegerVector aOrder(1, 2, math::GaussPointsMax());
- math_AreaFunction aFunction(theSurface);
- math_GaussMultipleIntegration anIntegral(aFunction, theLower, theUpper, aOrder);
- if (anIntegral.IsDone())
- {
- anIntegral.Dump(std::cout);
- }
- }
- void evalArea(const Handle(Geom_BoundedSurface)& theSurface)
- {
- math_IntegerVector aOrder(1, 2, math::GaussPointsMax());
- math_Vector aLower(1, 2, 0.0);
- math_Vector aUpper(1, 2, 0.0);
- theSurface->Bounds(aLower(1), aUpper(1), aLower(2), aUpper(2));
- math_AreaFunction aFunction(theSurface);
- math_GaussMultipleIntegration anIntegral(aFunction, aLower, aUpper, aOrder);
- if (anIntegral.IsDone())
- {
- anIntegral.Dump(std::cout);
- }
- }
- void testFace(const TopoDS_Shape& theFace)
- {
- GProp_GProps aSurfaceProps;
- BRepGProp::SurfaceProperties(theFace, aSurfaceProps);
- std::cout << "Face area: " << aSurfaceProps.Mass() << std::endl;
- }
- void testPlane()
- {
- std::cout << "====== Test Plane Area =====" << std::endl;
- Handle(Geom_Plane) aPlaneSurface = new Geom_Plane(gp::XOY());
- math_Vector aLower(1, 2);
- math_Vector aUpper(1, 2);
- // Parameter U range.
- aLower(1) = 0.0;
- aUpper(1) = 2.0;
- // Parameter V range.
- aLower(2) = 0.0;
- aUpper(2) = 3.0;
- evalArea(aPlaneSurface, aLower, aUpper);
- // Convert to BSpline Surface.
- Handle(Geom_RectangularTrimmedSurface) aTrimmedSurface =
- new Geom_RectangularTrimmedSurface(aPlaneSurface, aLower(1), aUpper(1), aLower(2), aUpper(2));
- Handle(Geom_BSplineSurface) aBSplineSurface = GeomConvert::SurfaceToBSplineSurface(aTrimmedSurface);
- evalArea(aBSplineSurface);
- // Test Face.
- TopoDS_Face aFace = BRepBuilderAPI_MakeFace(aTrimmedSurface, Precision::Confusion()).Face();
- testFace(aFace);
- aFace = BRepBuilderAPI_MakeFace(aBSplineSurface, Precision::Confusion()).Face();
- testFace(aFace);
- }
- void testCylinder()
- {
- std::cout << "====== Test Cylinder Area =====" << std::endl;
- Handle(Geom_CylindricalSurface) aCylindrialSurface = new Geom_CylindricalSurface(gp::XOY(), 1.0);
- math_Vector aLower(1, 2);
- math_Vector aUpper(1, 2);
- aLower(1) = 0.0;
- aUpper(1) = M_PI * 2.0;
- aLower(2) = 0.0;
- aUpper(2) = 3.0;
- evalArea(aCylindrialSurface, aLower, aUpper);
- // Convert to BSpline Surface.
- Handle(Geom_RectangularTrimmedSurface) aTrimmedSurface =
- new Geom_RectangularTrimmedSurface(aCylindrialSurface, aLower(1), aUpper(1), aLower(2), aUpper(2));
- Handle(Geom_BSplineSurface) aBSplineSurface = GeomConvert::SurfaceToBSplineSurface(aTrimmedSurface);
- evalArea(aBSplineSurface);
- // Test Face.
- TopoDS_Face aFace = BRepBuilderAPI_MakeFace(aTrimmedSurface, Precision::Confusion()).Face();
- testFace(aFace);
- aFace = BRepBuilderAPI_MakeFace(aBSplineSurface, Precision::Confusion()).Face();
- testFace(aFace);
- }
- void testSphere()
- {
- std::cout << "====== Test Sphere Area =====" << std::endl;
- Handle(Geom_SphericalSurface) aSphericalSurface = new Geom_SphericalSurface(gp::XOY(), 1.0);
- math_Vector aLower(1, 2);
- math_Vector aUpper(1, 2);
- aSphericalSurface->Bounds(aLower(1), aUpper(1), aLower(2), aUpper(2));
- evalArea(aSphericalSurface, aLower, aUpper);
- // Convert to BSpline Surface.
- Handle(Geom_BSplineSurface) aBSplineSurface = GeomConvert::SurfaceToBSplineSurface(aSphericalSurface);
- evalArea(aBSplineSurface);
- // Test Face.
- TopoDS_Face aFace = BRepBuilderAPI_MakeFace(aSphericalSurface, Precision::Confusion()).Face();
- testFace(aFace);
- aFace = BRepBuilderAPI_MakeFace(aBSplineSurface, Precision::Confusion()).Face();
- testFace(aFace);
- }
- void test()
- {
- testPlane();
- testSphere();
- testCylinder();
- }
- int main(int argc, char* argv[])
- {
- test();
- return 0;
- }
计算结果如下图所示:
上述代码计算了曲面的面积,再将曲面转换成B样条曲面,再使用算法计算面积。再将曲面和转换的B样条曲面生成拓朴面,利用OpenCASCADE中计算曲面面积功能进行对比。使用自定义函数math_AreaFunction利用多重积分类计算的结果与OpenCASCADE中计算曲面面积的值是一致的。当把曲面转换成B样条曲面后,OpenCASCADE计算的曲面面积偏大。
在学习《高等数学》的积分时,其主要的一个应用就是计算弧长、面积和体积等。学习高数抽象概念时,总会问学了高数有什么用?就从计算机图形方面来看,可以利用数学工具对任意曲线求弧长,对任意曲面计算面积等,更具一般性。
通过自定义被积函数再利用积分算法来计算任意曲面的面积,将理论与实践结合起来了。即将曲面的第一基本公式与具体的代码甚至可以利用OpenCASCADE生成对应的图形,这样抽象的理论就直观了,更便于理解相应的概念。
1.朱心雄. 自由曲线曲面造型技术. 科学出版社. 2000
2.陈维桓. 微分几何. 北京大学出版社. 2006
3.同济大学数学教研室. 高等数学. 高等教育出版社. 1996
来源: http://www.cnblogs.com/opencascade/p/OpenCASCADE_GaussIntegration.html