diff --git a/src/geometry/CCSTCurveBuilder.cpp b/src/geometry/CCSTCurveBuilder.cpp index eec79a713..301c6506a 100644 --- a/src/geometry/CCSTCurveBuilder.cpp +++ b/src/geometry/CCSTCurveBuilder.cpp @@ -20,6 +20,9 @@ #include "tiglmathfunctions.h" #include "CFunctionToBspline.h" +#include "CTiglError.h" + +#include "GeomAPI_PointsToBSpline.hxx" #include @@ -83,9 +86,15 @@ namespace namespace tigl { -CCSTCurveBuilder::CCSTCurveBuilder(double N1, double N2, const std::vector& B, double T) +CCSTCurveBuilder::CCSTCurveBuilder(double N1, double N2, const std::vector& 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 @@ std::vector CCSTCurveBuilder::B() const 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); + } } } diff --git a/src/geometry/CCSTCurveBuilder.h b/src/geometry/CCSTCurveBuilder.h index 665756dfc..3b9155948 100644 --- a/src/geometry/CCSTCurveBuilder.h +++ b/src/geometry/CCSTCurveBuilder.h @@ -35,7 +35,23 @@ namespace tigl class CCSTCurveBuilder { public: - TIGL_EXPORT CCSTCurveBuilder(double N1, double N2, const std::vector& 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& 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 _b; + int _degree; + double _tol; + Algorithm _algo; }; } // namespace tigl diff --git a/tests/unittests/tiglWingCSTProfile.cpp b/tests/unittests/tiglWingCSTProfile.cpp index c21435dd8..557564db1 100644 --- a/tests/unittests/tiglWingCSTProfile.cpp +++ b/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{}; + +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 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 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,7 +260,7 @@ 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; @@ -222,8 +268,9 @@ TEST(CSTApprox, hypersonicAirfoil) 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 + ) +);