OpenCASCADE 參數曲面面積

来源:http://www.cnblogs.com/opencascade/archive/2017/12/09/OpenCASCADE_GaussIntegration.html
-Advertisement-
Play Games

Abstract. 本文介紹了參數曲面的第一基本公式,並應用曲面的第一基本公式,結合OpenCASCADE中計算多重積分的類,對任意參數曲面的面積進行計算。 Key Words. Parametric Curve, Parametric Surface, Gauss Integration, Gl... ...


OpenCASCADE 參數曲面面積

[email protected]

Abstract. 本文介紹了參數曲面的第一基本公式,並應用曲面的第一基本公式,結合OpenCASCADE中計算多重積分的類,對任意參數曲面的面積進行計算。

Key Words. Parametric Curve, Parametric Surface, Gauss Integration, Global Properties

1.Introduction

我們知道一元函數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來介紹參數曲線曲面積分分別計算曲線弧長和曲面的面積。結合《微分幾何》來更好地理解曲線曲面相關知識。

2.Curve Natural Parametric Equations

設曲線C的參數方程是r=r(t),命:

wps_clip_image-10887

則s是該曲線的一個不變數,即它與空間中的坐標系的選擇無關,也與該曲線的參數變換無關。前者是因為在笛卡爾直角坐標變換下,切向量的長度|r’(t)|是不變的,故s不變。關於後者可以通過積分的變數的替換來證明,設參數變換是:

wps_clip_image-21189

並且

wps_clip_image-16009

因此

wps_clip_image-31864

根據積分的變數替換公式有:

wps_clip_image-29413

不變數s的幾何意義就是曲線段的弧長。這說明曲線參數t可以是任意的,但選擇不同的參數得到的參數方程會有不同,但是曲線段的弧長是不變的。以曲線弧長作為曲線方程的參數,這樣的方程稱為曲線的自然參數方程Natural Parametric Equations。

由曲線的參數方程可知,曲線弧長的計算公式為:

wps_clip_image-8988

幾何意義就是在每個微元處的切向量的長度求和。

3.Gauss Integration for Arc Length

曲線弧長的計算就是一元函數的積分。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;

};

 

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表示曲面上曲線的弧長,則由複合函數求導公式可得弧長微分公式:

wps_clip_image-11118

在古典微分幾何中,上式稱為曲面的第一基本公式,E、F、G稱為第一基本量。在曲面上,每一點的第一基本量與參數化無關。

利用曲面第一基本公式可以用於計算曲面的面積。參數曲面上與u,v參數空間的元素dudv對應的面積元為:

wps_clip_image-3399

由參數曲面法向的計算可知,曲面的面積元素即為u,v方向上的偏導數的乘積的模。

wps_clip_image-26176

其幾何意義可以理解為參數曲面的面積微元是由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);

}

}

 

通過theLower和theUpper指定定義域,由於採用了Gauss-Legendre演算法計算二重積分,所以需要指定階數,且階數越高積分結果精度越高,這裡使用了OpenCASCADE中最高的階數。

下麵通過對基本曲面的面積計算來驗證結果的正確性,並將計算結果和OpenCASCADE中計算面積的類BRepGProp::SurfaceProperties()結果進行對比。

6.Elementary Surface Area Test

下麵通過對OpenCASCADE中幾個初等曲面的面積進行計算,代碼如下所示:

/*
Copyright(C) 2017 Shing Liu([email protected])

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;
}

 

計算結果如下圖所示:

wps_clip_image-21244

上述代碼計算了曲面的面積,再將曲面轉換成B樣條曲面,再使用演算法計算面積。再將曲面和轉換的B樣條曲面生成拓樸面,利用OpenCASCADE中計算曲面面積功能進行對比。使用自定義函數math_AreaFunction利用多重積分類計算的結果與OpenCASCADE中計算曲面面積的值是一致的。當把曲面轉換成B樣條曲面後,OpenCASCADE計算的曲面面積偏大。

7.Conclusion

在學習《高等數學》的積分時,其主要的一個應用就是計算弧長、面積和體積等。學習高數抽象概念時,總會問學了高數有什麼用?就從電腦圖形方面來看,可以利用數學工具對任意曲線求弧長,對任意曲面計算面積等,更具一般性。

通過自定義被積函數再利用積分演算法來計算任意曲面的面積,將理論與實踐結合起來了。即將曲面的第一基本公式與具體的代碼甚至可以利用OpenCASCADE生成對應的圖形,這樣抽象的理論就直觀了,更便於理解相應的概念。

8.References

1.朱心雄. 自由曲線曲面造型技術. 科學出版社. 2000

2.陳維桓. 微分幾何. 北京大學出版社. 2006

3.同濟大學數學教研室. 高等數學. 高等教育出版社. 1996


您的分享是我們最大的動力!

-Advertisement-
Play Games
更多相關文章
  • 使用RSA公鑰解密,用openssl命令就是openssl rsautl -verify -in cipher_text -inkey public.pem -pubin -out clear_text,但其python網上還真沒有找到有博文去寫,只有hash的rsa解簽名。 這裡使用rsa庫,如果 ...
  • 一、python特性概要 1. python是解釋性腳本語言。 2. python特性總結 2.1 位元組碼 2.2 動態語義 在賦值是確定數據類型 2.3 縮進(4個空格) 3. python定義編碼類型 # -*- coding: utf-8 -*- 或 # -*- coding= utf-8 - ...
  • 題目描述 請你找出M個和為N的正整數,他們的乘積要儘可能的大。 輸出字典序最小的一種方案。 輸入輸出格式 輸入格式: 一行,兩個正整數N,M 輸出格式: M個和為N的,乘積儘可能的大的正整數。 輸入輸出樣例 輸入樣例#1: 複製 6 3 輸出樣例#1: 複製 2 2 2 輸入樣例#1: 複製 6 3 ...
  • 系統環境:win10 軟體:sublime Text 3 安裝插件: markdown editing.markdownlivepreview 修改prference-markdownlivepreview,複製左邊到右邊第一項設置為true。 下麵記錄了一些基本使用語法: #這是標題##這個是二級 ...
  • 前言 在javaEE中,最老生常談的莫過於SSH框架了,雖然現在類似struts的框架已經有很多了,但是我們還是不能夠忽視這個框架蘊含的一些知識。所以接下來的博文,筆者會根據struts的相關知識進行一些相關的講解。 一、struts框架的概述 1.1什麼是struts2? 以上是百度百科的對於st ...
  • 題目描述 Farmer John suffered a terrible loss when giant Australian cockroaches ate the entirety of his hay inventory, leaving him with nothing to feed th ...
  • 今天我們聊聊 Java 線程的中斷機制。 線程中斷機制提供了一種方法,用於將線程從阻塞等待中喚醒,並作出相應的“受控中斷”處理。 這段代碼使用了 Java 提供的 wait/notify 機制,線程執行 lock.wait() 會阻塞,有三種情況使線程恢復運行。 超時 1000ms 結束,正常執行下 ...
  • 題目背景 本題為題目 普通平衡樹 的可持久化加強版。 數據已經經過強化 題目描述 您需要寫一種數據結構(可參考題目標題),來維護一些數,其中需要提供以下操作(對於各個以往的歷史版本): 插入x數 刪除x數(若有多個相同的數,因只刪除一個,如果沒有請忽略該操作) 查詢x數的排名(排名定義為比當前數小的 ...
一周排行
    -Advertisement-
    Play Games
  • 移動開發(一):使用.NET MAUI開發第一個安卓APP 對於工作多年的C#程式員來說,近來想嘗試開發一款安卓APP,考慮了很久最終選擇使用.NET MAUI這個微軟官方的框架來嘗試體驗開發安卓APP,畢竟是使用Visual Studio開發工具,使用起來也比較的順手,結合微軟官方的教程進行了安卓 ...
  • 前言 QuestPDF 是一個開源 .NET 庫,用於生成 PDF 文檔。使用了C# Fluent API方式可簡化開發、減少錯誤並提高工作效率。利用它可以輕鬆生成 PDF 報告、發票、導出文件等。 項目介紹 QuestPDF 是一個革命性的開源 .NET 庫,它徹底改變了我們生成 PDF 文檔的方 ...
  • 項目地址 項目後端地址: https://github.com/ZyPLJ/ZYTteeHole 項目前端頁面地址: ZyPLJ/TreeHoleVue (github.com) https://github.com/ZyPLJ/TreeHoleVue 目前項目測試訪問地址: http://tree ...
  • 話不多說,直接開乾 一.下載 1.官方鏈接下載: https://www.microsoft.com/zh-cn/sql-server/sql-server-downloads 2.在下載目錄中找到下麵這個小的安裝包 SQL2022-SSEI-Dev.exe,運行開始下載SQL server; 二. ...
  • 前言 隨著物聯網(IoT)技術的迅猛發展,MQTT(消息隊列遙測傳輸)協議憑藉其輕量級和高效性,已成為眾多物聯網應用的首選通信標準。 MQTTnet 作為一個高性能的 .NET 開源庫,為 .NET 平臺上的 MQTT 客戶端與伺服器開發提供了強大的支持。 本文將全面介紹 MQTTnet 的核心功能 ...
  • Serilog支持多種接收器用於日誌存儲,增強器用於添加屬性,LogContext管理動態屬性,支持多種輸出格式包括純文本、JSON及ExpressionTemplate。還提供了自定義格式化選項,適用於不同需求。 ...
  • 目錄簡介獲取 HTML 文檔解析 HTML 文檔測試參考文章 簡介 動態內容網站使用 JavaScript 腳本動態檢索和渲染數據,爬取信息時需要模擬瀏覽器行為,否則獲取到的源碼基本是空的。 本文使用的爬取步驟如下: 使用 Selenium 獲取渲染後的 HTML 文檔 使用 HtmlAgility ...
  • 1.前言 什麼是熱更新 游戲或者軟體更新時,無需重新下載客戶端進行安裝,而是在應用程式啟動的情況下,在內部進行資源或者代碼更新 Unity目前常用熱更新解決方案 HybridCLR,Xlua,ILRuntime等 Unity目前常用資源管理解決方案 AssetBundles,Addressable, ...
  • 本文章主要是在C# ASP.NET Core Web API框架實現向手機發送驗證碼簡訊功能。這裡我選擇是一個互億無線簡訊驗證碼平臺,其實像阿裡雲,騰訊雲上面也可以。 首先我們先去 互億無線 https://www.ihuyi.com/api/sms.html 去註冊一個賬號 註冊完成賬號後,它會送 ...
  • 通過以下方式可以高效,並保證數據同步的可靠性 1.API設計 使用RESTful設計,確保API端點明確,並使用適當的HTTP方法(如POST用於創建,PUT用於更新)。 設計清晰的請求和響應模型,以確保客戶端能夠理解預期格式。 2.數據驗證 在伺服器端進行嚴格的數據驗證,確保接收到的數據符合預期格 ...