Add GeomAPI_PointsToBSpline as an optional alternative algorithm in CSTCurveBuilder

ChangeLog.md

@@ -3,14 +3,15 @@ Changelog
 
 Changes since last release
 -------------
-18/12/2024
+10/01/2025
 
 - Fixes
   - #936 A particular defined positioning (of the C++-type CCPACSPositioning) was not available via Python bindings since the std::vector<std::unique_ptr<**Type**>> is not exposed to swig. New getter functions have been implemented in CCPACSPositioning.h to make these elements accesible via index, similar to the implementation of for several other classes. For more information see https://github.com/RISCSoftware/cpacs_tigl_gen/issues/59.
 
 
 - General changes:
     - Update the C++ standard to C++17 (#1045).
+    - Added the possibility to switch between two algorithms in the `CCSTCurveBuilder`. The default algorithm based on a Chebychev approximation results in a small number of control points, but introduces small discontinuities in the curvature. The new option is to use OCCT's `GeomAPI_PointsToBSpline`, which results in a C3 continuous B-Spline.
 
 
 13/11/2024

src/geometry/CCSTCurveBuilder.cpp

@@ -20,6 +20,9 @@
 
 #include "tiglmathfunctions.h"
 #include "CFunctionToBspline.h"
+#include "CTiglError.h"
+
+#include "GeomAPI_PointsToBSpline.hxx"
 
 #include <cassert>
 
@@ -83,9 +86,15 @@
 namespace tigl
 {
 
-CCSTCurveBuilder::CCSTCurveBuilder(double N1, double N2, const std::vector<double>& B, double T)
+CCSTCurveBuilder::CCSTCurveBuilder(double N1, double N2, const std::vector<double>& B, double T, Algorithm method)
+    : _n1(N1)
+    , _n2(N2)
+    , _t(T)
+    , _b(B)
+    , _degree(4)
+    , _tol(1e-5)
+    , _algo(method)
 {
-    _n1 = N1; _n2 = N2; _b = B; _t = T;
 }
 
 double CCSTCurveBuilder::N1() const
@@ -111,8 +120,28 @@
 Handle(Geom_BSplineCurve) CCSTCurveBuilder::Curve()
 {
     CSTFunction function(this);
-    CFunctionToBspline approximator(function, 0., 1., 4, 1e-5, 10);
-    return approximator.Curve();
+    if (_algo == Algorithm::Piecewise_Chebychev_Approximation)
+    {
+        CFunctionToBspline approximator(function, 0., 1., _degree, _tol, 10);
+        return approximator.Curve();
+    }
+    else if (_algo == Algorithm::GeomAPI_PointsToBSpline) {
+
+        // sample the CST curve
+        int nsamples = 100;
+        TColgp_HArray1OfPnt points(1, nsamples);
+        for (int i = 0; i < nsamples; ++i) {
+            double t = (double)i*1/((double)nsamples-1);
+            points.SetValue(i+1, gp_Pnt(function.valueX(t), function.valueY(t), function.valueZ(t)));
+        }
+
+        // approximate sampled points using OCCT's internal algorithm GeomAPI_PointsToBSpline
+        GeomAPI_PointsToBSpline approximator(points, _degree, _degree, GeomAbs_C3, _tol);
+        return approximator.Curve();
+    }
+    else {
+        throw CTiglError("Unknown algorithm enum value passed to CCSTCurveBuilder", TIGL_ERROR);
+    }
 }
 
 }

src/geometry/CCSTCurveBuilder.h

@@ -35,7 +35,23 @@ namespace tigl
 class CCSTCurveBuilder
 {
 public:
-    TIGL_EXPORT CCSTCurveBuilder(double N1, double N2, const std::vector<double>& B, double T);
+
+    /**
+     * @brief The Algorithm enum gives a choice of the method.
+     *
+     * Piecewise_Chebychev_Approximation subdivides the curve into
+     * segments, where each segment is approximated using a Chebychev polynomials.
+     * The final result is the C1 concatenation of these polynomials to a B-Spline
+     *
+     * GeomAPI_PointsToBSpline uses OCCT's internal approximation algorithm to create
+     * a B-Spline that approximates a CST Curve.
+     */
+    enum class Algorithm {
+        Piecewise_Chebychev_Approximation = 0,
+        GeomAPI_PointsToBSpline
+    };
+
+    TIGL_EXPORT CCSTCurveBuilder(double N1, double N2, const std::vector<double>& B, double T, Algorithm method=Algorithm::Piecewise_Chebychev_Approximation);
 
     // returns parameters of cst curve
     TIGL_EXPORT double N1() const;
@@ -48,6 +64,9 @@ class CCSTCurveBuilder
 private:
     double _n1, _n2, _t;
     std::vector<double> _b;
+    int _degree;
+    double _tol;
+    Algorithm _algo;
 };
 
 } // namespace tigl

tests/unittests/tiglWingCSTProfile.cpp

@@ -152,7 +152,9 @@ TEST_F(WingCSTProfile, tiglWingCSTProfile_VTK_export)
     ASSERT_TRUE(tiglExportMeshedWingVTKSimpleByUID(tiglHandle, "CSTExample_W1", vtkWingFilename, 0.01) == TIGL_SUCCESS);
 }
 
-TEST(CSTApprox, simpleAirfoil)
+class CSTApprox : public testing::TestWithParam<tigl::CCSTCurveBuilder::Algorithm>{};
+
+TEST_P(CSTApprox, simpleAirfoil1)
 {
     double N1 = 0.5;
     double N2 = 1.0;
@@ -161,7 +163,7 @@ TEST(CSTApprox, simpleAirfoil)
     B.push_back(1.0);
     B.push_back(1.0);
 
-    tigl::CCSTCurveBuilder builder(N1, N2, B, 0.);
+    tigl::CCSTCurveBuilder builder(N1, N2, B, 0., GetParam());
     Handle(Geom_BSplineCurve) curve = builder.Curve();
 
     double devmax=0.0;
@@ -183,7 +185,51 @@ TEST(CSTApprox, simpleAirfoil)
     ASSERT_NEAR(0.0, devmax, 1e-5);
 }
 
-TEST(CSTApprox, ellipticBody)
+TEST_P(CSTApprox, simpleAirfoil2)
+{
+    // see https://github.com/DLR-SC/tigl/issues/1048
+    double N1 = 0.5;
+    double N2 = 1.0;
+    std::vector<double> B;
+    B.push_back(0.11809019);
+    B.push_back(0.18951797);
+    B.push_back(0.20255648);
+    double T = 0.0025;
+
+    auto algo = GetParam();
+    tigl::CCSTCurveBuilder builder(N1, N2, B, T, algo);
+    Handle(Geom_BSplineCurve) curve = builder.Curve();
+
+    double devmax=0.0;
+    // project sample points on curve and calculate distance
+    std::vector<double> x;
+    for (int i = 0; i < 100; ++i) {
+        x.push_back(double(i)/100.0);
+    }
+
+    for (unsigned int i = 0; i < x.size(); ++i) {
+        gp_Pnt samplePoint(x[i], Standard_Real(tigl::cstcurve(N1, N2, B, T, x[i])), 0);
+        GeomAPI_ProjectPointOnCurve projection(samplePoint, curve);
+        gp_Pnt projectedPoint=projection.NearestPoint();
+        double deviation=samplePoint.Distance(projectedPoint);
+        if (deviation >= devmax) {
+            devmax=deviation;
+        }
+    }
+
+    if (algo == tigl::CCSTCurveBuilder::Algorithm::Piecewise_Chebychev_Approximation) {
+        // I am not sure why this algorithm fails to satisfy the tolerance 1e-5 for this profile
+        ASSERT_NEAR(0.0, devmax, 1e-4);
+    }
+
+    if (algo == tigl::CCSTCurveBuilder::Algorithm::GeomAPI_PointsToBSpline) {
+        auto edge = BRepBuilderAPI_MakeEdge(curve);
+        BRepTools::Write(edge, "cst_edge.brep");
+        ASSERT_NEAR(0.0, devmax, 1e-5);
+    }
+}
+
+TEST_P(CSTApprox, ellipticBody)
 {
     double N1 = 0.5;
     double N2 = 0.5;
@@ -192,7 +238,7 @@ TEST(CSTApprox, ellipticBody)
     B.push_back(1.0);
     B.push_back(1.0);
 
-    tigl::CCSTCurveBuilder builder(N1, N2, B, 0.);
+    tigl::CCSTCurveBuilder builder(N1, N2, B, 0., GetParam());
     Handle(Geom_BSplineCurve) curve = builder.Curve();
 
     double devmax=0.0;
@@ -214,16 +260,17 @@ TEST(CSTApprox, ellipticBody)
     ASSERT_NEAR(0.0, devmax, 1e-5);
 }
 
-TEST(CSTApprox, hypersonicAirfoil)
+TEST_P(CSTApprox, hypersonicAirfoil)
 {
     double N1 = 1.0;
     double N2 = 1.0;
     std::vector<double> B;
     B.push_back(1.0);
     B.push_back(1.0);
     B.push_back(1.0);
-
-    tigl::CCSTCurveBuilder builder(N1, N2, B, 0.);
+
+    auto algo = GetParam();
+    tigl::CCSTCurveBuilder builder(N1, N2, B, 0.,  algo);
     Handle(Geom_BSplineCurve) curve = builder.Curve();
 
     double devmax=0.0;
@@ -242,6 +289,23 @@ TEST(CSTApprox, hypersonicAirfoil)
             devmax=deviation;
         }
     }
-    // approximation should be exact since we require only degree 2 spline
-    ASSERT_NEAR(0.0, devmax, 1e-12);
+
+    if (algo == tigl::CCSTCurveBuilder::Algorithm::Piecewise_Chebychev_Approximation) {
+        // approximation should be exact since we require only degree 2 spline
+        ASSERT_NEAR(0.0, devmax, 1e-12);
+    }
+
+    if (algo == tigl::CCSTCurveBuilder::Algorithm::GeomAPI_PointsToBSpline) {
+        // I doubt the optimzation algo used by GeomAPI_PointsToBSpline guarantees the exact solution.
+        ASSERT_NEAR(0.0, devmax, 1e-5);
+    }
 }
+
+INSTANTIATE_TEST_SUITE_P(
+    CSTApprox_DifferentAlgos,
+    CSTApprox,
+    testing::Values(
+        tigl::CCSTCurveBuilder::Algorithm::Piecewise_Chebychev_Approximation,
+        tigl::CCSTCurveBuilder::Algorithm::GeomAPI_PointsToBSpline
+    )
+);