上述代码的意思是直接对曲线的一阶导数的长度求积分,即是弧长。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;
};
4.First Fundamental Form of a Surface
曲面参数方程是个二元向量函数。根据《微分几何》中曲面的第一基本公式(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方向的偏导数的向量围成的一个四边形的面积,则整个曲面的面积即是对面积元素求积分。由于参数曲面有两个参数,所以若要计算曲面的面积,只需要对面积元素计算二重积分即可。
5.Gauss Integration for Area
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);
}
}