diff --git a/bindings/python_internal/configuration.i b/bindings/python_internal/configuration.i index c606dd196..80e66a57e 100644 --- a/bindings/python_internal/configuration.i +++ b/bindings/python_internal/configuration.i @@ -155,7 +155,7 @@ #include "generated/CPACSAtaChapterList.h" #include "CTiglWingProfileNACA.h" #include "generated/CPACSNacaProfile.h" -#include "CTiglNACA4Calculator.h" +#include "CTiglNACACalculator.h" #include "CPACSNacaProfile.h" #include "CCPACSDeck.h" #include "generated/CPACSDecks.h" @@ -598,7 +598,7 @@ class CCPACSWingRibsPositioning; %include "generated/CPACSCurvePointListXYZ.h" %include "CCPACSCurvePointListXYZ.h" %include "generated/CPACSProfileGeometry.h" -%include "CTiglNACA4Calculator.h" +%include "CTiglNACACalculator.h" %include "generated/CPACSNacaProfile.h" %include "CCPACSWingProfile.h" %include "CCPACSFuselageProfile.h" diff --git a/src/wing/CTiglNACA4Calculator.cpp b/src/wing/CTiglNACA4Calculator.cpp deleted file mode 100644 index 69f8ed15a..000000000 --- a/src/wing/CTiglNACA4Calculator.cpp +++ /dev/null @@ -1,238 +0,0 @@ -/* -* Copyright (C) 2007-2026 German Aerospace Center (DLR/SC) -* -* Created: 2026-01-14 Hannah Gedler -* -* Licensed under the Apache License, Version 2.0 (the "License"); -* you may not use this file except in compliance with the License. -* You may obtain a copy of the License at -* -* http://www.apache.org/licenses/LICENSE-2.0 -* -* Unless required by applicable law or agreed to in writing, software -* distributed under the License is distributed on an "AS IS" BASIS, -* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -* See the License for the specific language governing permissions and -* limitations under the License. -*/ -/** -* @file -* @brief Calculator for NACA wing profile coordinates -*/ - -#include "gp_Vec2d.hxx" -#include "tiglmathfunctions.h" -#include "CTiglNACA4Calculator.h" -#include "CFunctionToBspline.h" -#include "CTiglError.h" - -namespace -{ - // Leading-edge reparametrization: x(t) = (1+eps)*t*t/(t+eps). - // Removes the thickness distribution's sqrt(x) derivative singularity at the leading - // edge (x(0)=0, x'(0)=0, same as a plain t*t substitution close to t=0), but - unlike a - // plain t*t, which reshapes the parametrization over the *entire* chord - relaxes back - // towards the identity x(t)~=t once t significantly exceeds eps, so the curve fit away - // from the leading edge (and its knot placement) stays close to the original, direct-x - // parametrization. A plain t*t was found to shift curve representation enough at - // e.g. x=0.25 to break a sibling-component boolean fuse in specific geometries; eps=0.02 - // keeps the identity-like region starting well before that. - double leParam(double t) - { - const double eps = 0.02; - return (1. + eps) * t * t / (t + eps); - } -} - -namespace tigl{ - - CTiglNACA4Calculator::CTiglNACA4Calculator(double max_camber, double max_camber_position, double max_profile_thickness, double trailing_edge_thickness) - : max_camber(max_camber/100) - , max_camber_position(max_camber_position/10) - , max_profile_thickness(max_profile_thickness/100) - , trailing_edge_thickness_half(trailing_edge_thickness/2) - { - if(this->max_camber > 1 || this->max_camber < 0){ - throw CTiglError("error in CTiglNACA4Calculator The argument max_camber must be between 0 and 9."); - } - if(this->max_camber_position > 1 || this->max_camber_position < 0){ - throw CTiglError("error in CTiglNACA4Calculator The argument max_camber_position must be between 0 and 9."); - } - if(this->max_profile_thickness > 1 || this->max_profile_thickness < 0){ - throw CTiglError("error in CTiglNACA4Calculator max_profile_thicknessmust be between 0 and 99."); - } - } - - CTiglNACA4Calculator::CTiglNACA4Calculator(::std::string const& naca_code , const double te_thickness) - : CTiglNACA4Calculator(0.0, - 0.0, - 0.0, - te_thickness) - { - if (naca_code.size() != 4){ - throw CTiglError("error in CTiglNACA4Calculator: naca_code is not four digits long"); - } - try{ - double m = static_cast(naca_code[0] - '0'); - double p = static_cast(naca_code[1] - '0'); - double t = static_cast(std::stoi(naca_code.substr(2,2))); - - *this = CTiglNACA4Calculator(m, p, t, te_thickness); - } - catch(...){ - throw CTiglError("error in CTiglNACA4Calculator: the naca_code format is not correct, it must to contain four digits and nothing else"); - } - - } - - double CTiglNACA4Calculator::get_trailing_edge_thickness() const - { - return trailing_edge_thickness_half; - } - - double CTiglNACA4Calculator::camberline(double x) const{ - - double m = this->max_camber; - double p = this->max_camber_position; - if(p == 0){ - return 0; - } - if(x < 0 || x > 1){ - throw CTiglError("error in CTiglNACA4Calculator::camberline: x must be between 0 and 1."); - } - if (x <= p) { - return (2*p*x - x*x)*m/(p*p);; - } - else { - return (1 - 2*p + 2*p*x - x*x)*m/((1-p)*(1-p));; - } - } - - gp_Vec2d CTiglNACA4Calculator::upper_curve(double x) const{ - double yt = profile_thickness(x); - double yc = camberline(x); - auto point = gp_Vec2d{x, yc}; - return point + yt*normal(x); - } - - gp_Vec2d CTiglNACA4Calculator::lower_curve(double x) const{ - double yt = profile_thickness(x); - double yc = camberline(x); - auto point = gp_Vec2d{x, yc}; - return point - yt*normal(x); - } - - double CTiglNACA4Calculator::profile_thickness(double x) const{ - double t = this->max_profile_thickness; - double e = -(trailing_edge_thickness_half/(5*t)) + (0.2969 - 0.1260 - 0.3516 + 0.2843); - return 5*t*(0.2969*sqrt(x) - 0.1260*x - 0.3516*(x*x)+0.2843*pow(x,3) - e*pow(x,4)); - } - - double CTiglNACA4Calculator::camberline_derivative(double x) const{ //c'(x) - double m = this->max_camber; - double p = this->max_camber_position; - - if(p == 0){ - return 0; - } - if(0 <= x && x <= p){ - return (2*p - 2*x)*m/(p*p); - } - else if(x > p){ - return ( 2*p - 2*x)*m/((1-p)*(1-p));; - } - else{ - throw CTiglError("error in CTiglNACA4Calculator::camberline_derivative: x must be between 0 and 1."); - } - } - - gp_Vec2d CTiglNACA4Calculator::normal(double x) const{ - gp_Vec2d normal = {-camberline_derivative(x), 1}; - double nrm = normal.Magnitude(); - if (nrm < 1e-12) { - return gp_Vec2d{0,1}; - } - return normal/nrm; - - } - - Handle(Geom_BSplineCurve) CTiglNACA4Calculator::upper_bspline() const{ - - CTiglNACA4UpperCurve upperCurve(*this); - - const double umin = 0.; - const double umax = 1.; - int degree = 3; - double tolerance=1e-5; - // CTiglNACA4UpperCurve/LowerCurve reparametrize near the leading edge (see - // leParam above), removing the thickness distribution's sqrt(x) derivative - // singularity there, so the adaptive Chebyshev fit now converges quickly - // everywhere (verified empirically: converges well below this depth, leaving no - // measurable tangent mismatch at any internal knot, including at the leading edge). - int maxDepth = 10; - - tigl::CFunctionToBspline converter(upperCurve, umin, umax, degree, tolerance, maxDepth); - return converter.Curve(); - } - - Handle(Geom_BSplineCurve) CTiglNACA4Calculator::lower_bspline() const{ - CTiglNACA4LowerCurve lowerCurve(*this); - - const double umin = 0.; - const double umax = 1.; - int degree = 3; - double tolerance=1e-5; - // CTiglNACA4UpperCurve/LowerCurve reparametrize near the leading edge (see - // leParam above), removing the thickness distribution's sqrt(x) derivative - // singularity there, so the adaptive Chebyshev fit now converges quickly - // everywhere (verified empirically: converges well below this depth, leaving no - // measurable tangent mismatch at any internal knot, including at the leading edge). - int maxDepth = 10; - - tigl::CFunctionToBspline converter(lowerCurve, umin, umax, degree, tolerance, maxDepth); - return converter.Curve(); - } - - CTiglNACA4UpperCurve::CTiglNACA4UpperCurve( CTiglNACA4Calculator const& calculator) - : MathFunc3d(), - calculator(calculator) - {} - - double CTiglNACA4UpperCurve::valueX(double t) { - // see leParam above: removes the thickness distribution's unbounded - // d(upper_curve)/dx at the leading edge (x=0) without perturbing the - // parametrization away from it. - gp_Vec2d vec = calculator.upper_curve(leParam(t)); - return vec.X(); - } - double CTiglNACA4UpperCurve::valueY(double t) { - return 0.0; - } - double CTiglNACA4UpperCurve::valueZ(double t) { - gp_Vec2d vec = calculator.upper_curve(leParam(t)); - return vec.Y(); - } - - CTiglNACA4LowerCurve::CTiglNACA4LowerCurve( CTiglNACA4Calculator const& calculator) - : MathFunc3d(), - calculator(calculator) - {} - - double CTiglNACA4LowerCurve::valueX(double t) { - // see CTiglNACA4UpperCurve::valueX/leParam for why t is reparametrized here - gp_Vec2d vec = calculator.lower_curve(leParam(t)); - return vec.X(); - } - - double CTiglNACA4LowerCurve::valueY(double t) { - return 0.0; - } - - double CTiglNACA4LowerCurve::valueZ(double t) { - gp_Vec2d vec = calculator.lower_curve(leParam(t)); - return vec.Y(); - } - -} //namespace tigl - - diff --git a/src/wing/CTiglNACACalculator.cpp b/src/wing/CTiglNACACalculator.cpp new file mode 100644 index 000000000..d4dec0a39 --- /dev/null +++ b/src/wing/CTiglNACACalculator.cpp @@ -0,0 +1,509 @@ +/* +* Copyright (C) 2007-2026 German Aerospace Center (DLR/SC) +* +* Created: 2026-01-14 Hannah Gedler +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +/** +* @file +* @brief Calculator for NACA wing profile coordinates +*/ + +#include "gp_Vec2d.hxx" +#include "tiglmathfunctions.h" +#include "CTiglNACACalculator.h" +#include "CFunctionToBspline.h" +#include "CTiglBSplineApproxInterp.h" +#include "CTiglError.h" +#include + + +namespace +{ + // Reparametrization x(t) = (1+eps)*t²/(t+eps) removes the leading-edge thickness singularity without distorting the rest of the chord like a plain t² would. + double leParam(double t) + { + const double eps = 0.02; + return (1. + eps) * t * t / (t + eps); + } +} + +namespace tigl{ + + NACA4DigitCode::NACA4DigitCode(std::string const& code) { + if (code.size() != 4) { + throw CTiglError("NACA4DigitCode: requires exactly 4 digits"); + } + + max_camber = (code[0] - '0') / 100.0; + max_camber_position = (code[1] - '0') / 10.0; + max_profile_thickness = std::stoi(code.substr(2, 2)) / 100.0; + } + + NACA5DigitCode::NACA5DigitCode(std::string const& code) { + if (code.size() != 5) { + throw CTiglError("NACA5DigitCode: requires exactly 5 digits"); + } + + max_camber = (code[0] - '0') / 100.0; + max_camber_position = (code[1] - '0') / 20.0; + reflex = code[2] - '0'; + max_profile_thickness = std::stoi(code.substr(3, 2)) / 100.0; + } + + CTiglNACACalculator::CTiglNACACalculator(const NACA4DigitCode& code, double trailing_edge_thickness) + : series_(Series::NACA4) + , max_camber(code.max_camber) + , max_camber_position(code.max_camber_position) + , max_profile_thickness(code.max_profile_thickness) + , trailing_edge_thickness_half(trailing_edge_thickness / 2.0) + { + if(this->max_camber > 1 || this->max_camber < 0){ + throw CTiglError("error in CTiglNACACalculator The argument max_camber must be between 0 and 9."); + } + if(this->max_camber_position > 1 || this->max_camber_position < 0){ + throw CTiglError("error in CTiglNACACalculator The argument max_camber_position must be between 0 and 9."); + } + if(this->max_profile_thickness > 1 || this->max_profile_thickness < 0){ + throw CTiglError("error in CTiglNACACalculator max_profile_thicknessmust be between 0 and 99."); + } + } + + CTiglNACACalculator::CTiglNACACalculator(const NACA5DigitCode& code, double trailing_edge_thickness) + : series_(Series::NACA5) + , max_camber(code.max_camber) + , max_camber_position(code.max_camber_position) + , reflex(code.reflex) + , max_profile_thickness(code.max_profile_thickness) + , trailing_edge_thickness_half(trailing_edge_thickness / 2.0) + { + if(this->max_camber > 1 || this->max_camber < 0){ + throw CTiglError("error in CTiglNACACalculator The argument max_camber must be between 0 and 9."); + } + if(this->max_camber_position > 1 || this->max_camber_position < 0){ + throw CTiglError("error in CTiglNACACalculator The argument max_camber_position must be between 0 and 9."); + } + if(this->max_profile_thickness > 1 || this->max_profile_thickness < 0){ + throw CTiglError("error in CTiglNACACalculator max_profile_thicknessmust be between 0 and 99."); + } + if(this->reflex != 0 && this->reflex != 1){ + throw CTiglError("error in CTiglNACACalculator The argument reflex must be 0 or 1."); + } + } + + double CTiglNACACalculator::get_trailing_edge_thickness() const + { + return trailing_edge_thickness_half; + } + + double CTiglNACACalculator::k1_const(double max_camber, double max_camber_position, double reflex) const{ + int max_camber_cl_whole = static_cast(max_camber*100); + int max_camber_position_whole = static_cast(max_camber_position*20); + int reflex_int = static_cast(reflex); + double k1; + std::string meanline_designation_str = std::to_string(max_camber_cl_whole) + std::to_string(max_camber_position_whole) + std::to_string(reflex_int); + int meanline_designation = std::stoi(meanline_designation_str); + switch(meanline_designation){ + case 210: + k1 =361.400; + break; + case 220: + k1 =51.640; + break; + case 230: + k1 = 15.957; + break; + case 240: + k1 = 6.643; + break; + case 250: + k1 = 3.230; + break; + case 211: + throw ::std::logic_error("error in CTiglNACACalculator::k1: this profile does not provide a constant for k1."); + break; + case 221: + k1 = 51.99; + break; + case 231: + k1 = 15.793; + break; + case 241: + k1 = 6.520; + break; + case 251: + k1 = 3.191; + break; + default: + throw ::std::logic_error("error in CTiglNACACalculator::k1: this profile does not provide a constant for k1."); + break; + } + return k1; + } + + double CTiglNACACalculator::m_const(double max_camber, double max_camber_position, double reflex) const{ + int max_camber_cl_whole = static_cast(max_camber*100); + int max_camber_position_whole = static_cast(max_camber_position*20); + int reflex_int = static_cast(reflex); + double m; + std::string meanline_designation_str = std::to_string(max_camber_cl_whole) + std::to_string(max_camber_position_whole) + std::to_string(reflex_int); + int meanline_designation = std::stoi(meanline_designation_str); + switch(meanline_designation){ + case 210: + m = 0.0580; + break; + case 220: + m = 0.1260; + break; + case 230: + m = 0.2025; + break; + case 240: + m = 0.2900; + break; + case 250: + m = 0.3910; + break; + case 211: + throw ::std::logic_error("error in CTiglNACACalculator::m: this profile does not provide a constant for m."); + break; + case 221: + m = 0.1300; + break; + case 231: + m = 0.2170; + break; + case 241: + m = 0.3180; + break; + case 251: + m = 0.4410; + break; + default: + throw ::std::logic_error("error in CTiglNACACalculator::m: this profile does not provide a constant for m."); + break; + } + return m; + } + + double CTiglNACACalculator::k2k1_const(double max_camber, double max_camber_position, double reflex) const{ + int max_camber_cl_whole = static_cast(max_camber*100); + int max_camber_position_whole = static_cast(max_camber_position*20); + int reflex_int = static_cast(reflex); + double k2k1; + std::string meanline_designation_str = std::to_string(max_camber_cl_whole) + std::to_string(max_camber_position_whole) + std::to_string(reflex_int); + int meanline_designation = std::stoi(meanline_designation_str); + switch(meanline_designation){ + case 210: + throw ::std::logic_error("error in CTiglNACACalculator::k2k1: this profile does not provide a constant for k2/k1."); + break; + case 220: + throw ::std::logic_error("error in CTiglNACACalculator::k2k1: this profile does not provide a constant for k2/k1."); + break; + case 230: + throw ::std::logic_error("error in CTiglNACACalculator::k2k1: this profile does not provide a constant for k2/k1m."); + break; + case 240: + throw ::std::logic_error("error in CTiglNACACalculator::k2k1: this profile does not provide a constant for k2/k1."); + break; + case 250: + throw ::std::logic_error("error in CTiglNACACalculator::k2k1: this profile does not provide a constant for k2/k1."); + break; + case 211: + throw ::std::logic_error("error in CTiglNACACalculator::k2k1: this profile does not provide a constant for k2/k1."); + break; + case 221: + k2k1 = 0.000764; + break; + case 231: + k2k1 = 0.00677; + break; + case 241: + k2k1 = 0.0303; + break; + case 251: + k2k1 = 0.1355; + break; + default: + throw ::std::logic_error("error in CTiglNACACalculator::k2k1: this profile does not provide a constant for m."); + break; + } + return k2k1; + } + + double CTiglNACACalculator::camberline(double x) const{ + + if(series_ == Series::NACA4){ + double m = this->max_camber; + double p = this->max_camber_position; + if(p == 0){ + return 0; + } + if(x < 0 || x > 1){ + throw CTiglError("error in CTiglNACACalculator::camberline::NACA4: x must be between 0 and 1."); + } + else if(x <= p){ + return (2*p*x - x*x)*m/(p*p); + } + else if(x > p){ + return (1 - 2*p + 2*p*x - x*x)*m/((1-p)*(1-p)); + } + } + else if(series_ == Series::NACA5){ + double s = this->max_camber; + + double p = this->max_camber_position; + + double q = this->reflex; + + double k1 = k1_const(s, p, q); + + double m = m_const(s, p, q); + + if(q == 0){ + + if(p == 0){ + return 0; + } + if(0 <= x && x <= m){ + double result = (k1/6)*(x*x*x-3*m*x*x+m*m*(3-m)*x); + return result; + } + else if(x > m){ + double result3 = ((k1*m*m*m)/6)*(1-x); + return result3; + } + } + else if(q == 1){ + double frack2k1 = k2k1_const(s,p,q); + + if(p == 0){ + return 0.0; + } + if(0 <= x && x <= m){ + double result1 = (k1/6)*(pow(x-m, 3) - frack2k1*pow(1-m, 3)*x - pow(m,3)*x + pow(m,3)); + + return result1; + } + else if(x > m){ + + double result2 = (k1/6)*(frack2k1*((x-m)*(x-m)*(x-m))-frack2k1*((1-m)*(1-m)*(1-m))*x-(m*m*m)*x+(m*m*m)); + return result2; + } + } + else{ + throw ::std::logic_error("error in CTiglNACACalculator::camberline::NACA5: x must be between 0 and 1."); + } + return 0; + } + return 0; + } + + + gp_Vec2d CTiglNACACalculator::upper_curve(double x) const{ + double yt = profile_thickness(x); + double yc = camberline(x); + auto point = gp_Vec2d{x, yc}; + gp_Vec2d point_calculated = point + yt*normal(x); + return point + yt*normal(x); + } + + gp_Vec2d CTiglNACACalculator::lower_curve(double x) const{ + double yt = profile_thickness(x); + double yc = camberline(x); + auto point = gp_Vec2d{x, yc}; + return point - yt*normal(x); + } + + double CTiglNACACalculator::profile_thickness(double x) const{ + double t = this->max_profile_thickness; + double e = -(trailing_edge_thickness_half/(5*t)) + (0.2969 - 0.1260 - 0.3516 + 0.2843); + return 5*t*(0.2969*sqrt(x) - 0.1260*x - 0.3516*(x*x)+0.2843*pow(x,3) - e*pow(x,4)); + } + + double CTiglNACACalculator::camberline_derivative(double x) const{ //c'(x) + + if(series_ == Series::NACA4){ + double m = this->max_camber; + double p = this->max_camber_position; + + if(p == 0){ + return 0; + } + if(0 <= x && x <= p){ + return (2*p - 2*x)*m/(p*p); + } + else if(x > p){ + return ( 2*p - 2*x)*m/((1-p)*(1-p));; + } + else{ + throw CTiglError("error in CTiglNACACalculator::camberline_derivative: x must be between 0 and 1."); + } + } + else if(series_ == Series::NACA5){ + double s = this->max_camber; + double p = this->max_camber_position; + double q = this->reflex; + double k1 = k1_const(s, p, q); + double m = m_const(s, p, q); + + if(p == 0){ + return 0; + } + if(q == 0){ + if(0 <= x && x <= m){ + double result6 = (k1/6)*(3*x*x-6*m*x+m*m*(3-m)); + return result6; + } + else if(x > m){ + return -((k1*m*m*m)/6); + } + } + else if(q==1){ + double frack2k1 = (3*((m-p)*(m-p))-m*m*m)/((1-m)*(1-m)*(1-m)); + if(0 <= x && x <= m){ + + double result5 = ((k1/6)*(3*((x-m)*(x-m))-frack2k1*((1-m)*(1-m)*(1-m))-m*m*m)); + return result5; + } + else if(x > m){ + + return (k1/6)*(3*frack2k1*((x-m)*(x-m))-frack2k1*((1-m)*(1-m)*(1-m))-m*m*m); + } + } + else{ + throw CTiglError("error in CTiglNACACalculator::camberline_derivative: x must be between 0 and 1."); + } + return 0; + } + return 0; + } + + gp_Vec2d CTiglNACACalculator::normal(double x) const{ + gp_Vec2d normal = {-camberline_derivative(x), 1}; + double nrm = normal.Magnitude(); + if (nrm < 1e-12) { + return gp_Vec2d{0,1}; + } + return normal/nrm; + } + + Handle(Geom_BSplineCurve) CTiglNACACalculator::upper_bspline() const{ + /* + int npnts = 200; + auto pnts = TColgp_Array1OfPnt(1, npnts); + for (int i=1; i<=npnts; ++i) { + double u = ((double)i-1)/(npnts-1); + //double x = 0.5 * (1.0 - cos(M_PI * u)); + double x = pow(0.5 * (1.0 - cos(M_PI * u)), 0.5); + auto p = upper_curve(x); + pnts.SetValue(i, gp_Pnt(p.X(), 0., p.Y())); + std::cout << p.X() << ", " << p.Y() << "\n"; + } + + int nControlPoints = 34; + int deg = 3; + bool continuous_if_closed = false; + auto builder = CTiglBSplineApproxInterp(pnts, nControlPoints, deg, continuous_if_closed); + builder.InterpolatePoint(0); + builder.InterpolatePoint(npnts-1); + auto result = builder.FitCurveOptimal(); + return result.curve; + */ + + CTiglNACAUpperCurve upperCurve(*this); + + const double umin = 0.; + const double umax = 1.; + int degree = 3; + double tolerance=1e-5; + int maxDepth = 10; + + tigl::CFunctionToBspline conv(upperCurve, umin, umax, degree, tolerance, maxDepth); + return conv.Curve(); + } + + Handle(Geom_BSplineCurve) CTiglNACACalculator::lower_bspline() const{ + /* + int npnts = 200; + auto pnts = TColgp_Array1OfPnt(1, npnts); + for (int i=1; i<=npnts; ++i) { + double u = ((double)i-1)/(npnts-1); + double x = pow(0.5 * (1.0 - cos(M_PI * u)), 0.5); + //double x = 0.5 * (1.0 - cos(M_PI * u)); + auto p = lower_curve(x); + pnts.SetValue(i, gp_Pnt(p.X(), 0., p.Y())); + } + + int nControlPoints = 34; + int deg = 3; + bool continuous_if_closed = false; + auto builder = CTiglBSplineApproxInterp(pnts, nControlPoints, deg, continuous_if_closed); + builder.InterpolatePoint(0); + builder.InterpolatePoint(npnts-1); + auto result = builder.FitCurveOptimal(); + return result.curve; + */ + + CTiglNACALowerCurve lowerCurve(*this); + + const double umin = 0.; + const double umax = 1.; + int degree = 3; + double tolerance=1e-5; + int maxDepth = 10; + + tigl::CFunctionToBspline conv(lowerCurve, umin, umax, degree, tolerance, maxDepth); + return conv.Curve(); + } + + CTiglNACAUpperCurve::CTiglNACAUpperCurve( CTiglNACACalculator const& calculator) + : MathFunc3d(), + calculator(calculator) + {} + + double CTiglNACAUpperCurve::valueX(double t) { + gp_Vec2d vec = calculator.upper_curve(leParam(t)); + return vec.X(); + } + double CTiglNACAUpperCurve::valueY(double t) { + return 0.0; + } + double CTiglNACAUpperCurve::valueZ(double t) { + gp_Vec2d vec = calculator.upper_curve(leParam(t)); + return vec.Y(); + } + + CTiglNACALowerCurve::CTiglNACALowerCurve( CTiglNACACalculator const& calculator) + : MathFunc3d(), + calculator(calculator) + {} + + double CTiglNACALowerCurve::valueX(double t) { + gp_Vec2d vec = calculator.lower_curve(leParam(t)); + return vec.X(); + } + + double CTiglNACALowerCurve::valueY(double t) { + return 0.0; + } + + double CTiglNACALowerCurve::valueZ(double t) { + gp_Vec2d vec = calculator.lower_curve(leParam(t)); + return vec.Y(); + } + +} //namespace tigl + + diff --git a/src/wing/CTiglNACA4Calculator.h b/src/wing/CTiglNACACalculator.h similarity index 73% rename from src/wing/CTiglNACA4Calculator.h rename to src/wing/CTiglNACACalculator.h index 1f5ab8310..d6c1cc2d9 100644 --- a/src/wing/CTiglNACA4Calculator.h +++ b/src/wing/CTiglNACACalculator.h @@ -26,23 +26,39 @@ #include #include "CFunctionToBspline.h" - namespace tigl{ -class CTiglNACA4Calculator{ + struct NACA4DigitCode { + double max_camber; + double max_camber_position; + double max_profile_thickness; + + TIGL_EXPORT explicit NACA4DigitCode(std::string const& code); +}; + + struct NACA5DigitCode { + double max_camber; + double max_camber_position; + double reflex; + double max_profile_thickness; + + TIGL_EXPORT explicit NACA5DigitCode(std::string const& code); +}; + +class CTiglNACACalculator{ public: /** - * @brief Construct a new CTiglNACA4Calculator object + * @brief Construct a new CTiglNACACalculator object * * @param max_camber * @param max_camber_position * @param max_profile_thickness */ - TIGL_EXPORT CTiglNACA4Calculator(double max_camber, double max_camber_position, double max_profile_thickness, double trailing_edge_thickness = 0); + TIGL_EXPORT CTiglNACACalculator(const NACA4DigitCode& code, double trailing_edge_thickness = 0.); + TIGL_EXPORT CTiglNACACalculator(const NACA5DigitCode& code, double trailing_edge_thickness = 0.); - explicit TIGL_EXPORT CTiglNACA4Calculator(const ::std::string& naca_code, const double te_thickness); /** * @brief getter for the trailing edge thickness @@ -82,6 +98,9 @@ class CTiglNACA4Calculator{ * @return double */ TIGL_EXPORT double profile_thickness(double x) const; + TIGL_EXPORT double m_const(double max_camber, double max_camber_position, double reflex) const; + TIGL_EXPORT double k1_const(double max_camber, double max_camber_position, double reflex) const; + TIGL_EXPORT double k2k1_const(double max_camber, double max_camber_position, double reflex) const; TIGL_EXPORT Handle(Geom_BSplineCurve) upper_bspline() const; @@ -104,14 +123,18 @@ class CTiglNACA4Calculator{ */ gp_Vec2d normal(double x) const; + + private: - double max_camber, max_camber_position, max_profile_thickness, trailing_edge_thickness_half; + double max_camber, max_camber_position, max_profile_thickness, trailing_edge_thickness_half, reflex; + enum class Series {Unknown, NACA4, NACA5}; + Series series_ = Series::Unknown; }; -class CTiglNACA4UpperCurve : public MathFunc3d { - public: - TIGL_EXPORT explicit CTiglNACA4UpperCurve( CTiglNACA4Calculator const& calculator); +class CTiglNACAUpperCurve : public MathFunc3d { + public: + TIGL_EXPORT explicit CTiglNACAUpperCurve( CTiglNACACalculator const& calculator); /** * @brief Get the X coordinate of the upper curve @@ -142,17 +165,17 @@ class CTiglNACA4UpperCurve : public MathFunc3d { double valueZ(double t) override; private: - CTiglNACA4Calculator const& calculator; + CTiglNACACalculator const& calculator; }; -class CTiglNACA4LowerCurve : public MathFunc3d { - public: - TIGL_EXPORT explicit CTiglNACA4LowerCurve( CTiglNACA4Calculator const& calculator); +class CTiglNACALowerCurve : public MathFunc3d { + public: + TIGL_EXPORT explicit CTiglNACALowerCurve( CTiglNACACalculator const& calculator); /** * @brief Get the X coordinate of the lower curve * - * Note: t is not the chord fraction x directly - see CTiglNACA4UpperCurve::valueX. + * Note: t is not the chord fraction x directly - see CTiglNACAUpperCurve::valueX. * * @param t * @return double @@ -176,6 +199,6 @@ class CTiglNACA4LowerCurve : public MathFunc3d { double valueZ(double t) override; private: - CTiglNACA4Calculator const& calculator; + CTiglNACACalculator const& calculator; }; } \ No newline at end of file diff --git a/src/wing/CTiglWingProfileNACA.cpp b/src/wing/CTiglWingProfileNACA.cpp index e5f5e92fe..21766d33b 100644 --- a/src/wing/CTiglWingProfileNACA.cpp +++ b/src/wing/CTiglWingProfileNACA.cpp @@ -17,7 +17,7 @@ */ /** * @file -* @brief builds the wing profile for NACA profiles based on the CTiglNACA4Calculator and returns the upper and lower wire, the trailing edge and the leading edge point +* @brief builds the wing profile for NACA profiles based on the CTiglNACACalculator and returns the upper and lower wire, the trailing edge and the leading edge point */ @@ -26,7 +26,7 @@ #include "ITiglWingProfileAlgo.h" #include "Cache.h" #include "geometry/CFunctionToBspline.h" -#include "CTiglNACA4Calculator.h" +#include "CTiglNACACalculator.h" #include "common/tiglcommonfunctions.h" @@ -36,7 +36,8 @@ #include #include #include "CWireToCurve.h" - +#include +#include namespace tigl { @@ -44,10 +45,19 @@ namespace tigl CTiglWingProfileNACA::CTiglWingProfileNACA(const CCPACSWingProfile& profile, const generated::CPACSNacaProfile& nacadef) : profileUID(profile.GetUID()) - , calculator(nacadef.GetNaca4DigitCode_choice1() ? CTiglNACA4Calculator(*nacadef.GetNaca4DigitCode_choice1(), nacadef.GetTrailingEdgeThickness() ? *nacadef.GetTrailingEdgeThickness() : 0.0) : throw CTiglError("ERROR in CTiglWingProfileNACA: Currently only 4 digit NACA codes implemented.") ) + , nacacode(std::in_place_type_t{}, "0000") + , te_thickness(nacadef.GetTrailingEdgeThickness() ? *nacadef.GetTrailingEdgeThickness() : 0.0) , wireCache(*this, &CTiglWingProfileNACA::BuildWires) { - + if (nacadef.GetNaca4DigitCode_choice1()) { + nacacode = NACA4DigitCode(*nacadef.GetNaca4DigitCode_choice1()); + } + else if (nacadef.GetNaca5DigitCode_choice2()) { + nacacode = NACA5DigitCode(*nacadef.GetNaca5DigitCode_choice2()); + } + else { + throw CTiglError("ERROR in CTiglWingProfileNACA: No valid NACA code provided. Expected NACA 4-digit or 5-digit code."); + } } void CTiglWingProfileNACA::Invalidate() const @@ -57,6 +67,11 @@ void CTiglWingProfileNACA::Invalidate() const void CTiglWingProfileNACA::BuildWires(WireCache& cache) const { + + CTiglNACACalculator calculator = std::visit( + [&](auto const& code) { return CTiglNACACalculator(code, te_thickness); }, + nacacode + ); auto upper_bspline = calculator.upper_bspline(); auto lower_bspline = calculator.lower_bspline(); @@ -71,7 +86,6 @@ void CTiglWingProfileNACA::BuildWires(WireCache& cache) const cache.lePoint = gp_Pnt(le_pnt.X(), 0.0, le_pnt.Y()); cache.tePoint = gp_Pnt(te_pnt.X(), 0.0, te_pnt.Y()); - // build trailing edge if (HasBluntTE()) { double upper_x_coord = upper_coord.X(); @@ -121,7 +135,6 @@ const TopoDS_Edge& CTiglWingProfileNACA::GetLowerWire(TiglShapeModifier mod) con // gets the upper and lower wing profile into on edge const TopoDS_Edge& CTiglWingProfileNACA::GetUpperLowerWire(TiglShapeModifier mod) const { - double te_thickness = calculator.get_trailing_edge_thickness(); if (mod == SHARP_TRAILINGEDGE && HasBluntTE()) { LOG(WARNING) << "Profile " << profileUID << ": SHARP_TRAILINGEDGE modifier requested but profile has blunt trailing edge (thickness: " << te_thickness << ")"; @@ -154,7 +167,6 @@ const gp_Pnt & CTiglWingProfileNACA::GetTEPoint() const bool CTiglWingProfileNACA::HasBluntTE() const { - return calculator.get_trailing_edge_thickness() > 0.; - + return te_thickness > 0.; } }//namespace tigl \ No newline at end of file diff --git a/src/wing/CTiglWingProfileNACA.h b/src/wing/CTiglWingProfileNACA.h index 776d2e1ae..7e145645f 100644 --- a/src/wing/CTiglWingProfileNACA.h +++ b/src/wing/CTiglWingProfileNACA.h @@ -25,7 +25,7 @@ #include "tigl_internal.h" #include "ITiglWingProfileAlgo.h" #include "Cache.h" -#include "CTiglNACA4Calculator.h" +#include "CTiglNACACalculator.h" #include #include @@ -111,12 +111,12 @@ namespace tigl{ private: struct WireCache { - TopoDS_Edge upperWire; /**< wire of the upper wing profile */ - TopoDS_Edge lowerWire; /**< wire of the lower wing profile */ + TopoDS_Edge upperWire; + TopoDS_Edge lowerWire; TopoDS_Edge trailingEdge; TopoDS_Edge upperLowerWire; - gp_Pnt lePoint; /**< Leading edge point */ - gp_Pnt tePoint; /**< Trailing edge point */ + gp_Pnt lePoint; + gp_Pnt tePoint; }; void InvalidateParent() const; @@ -126,7 +126,7 @@ namespace tigl{ private: std::string profileUID; - CTiglNACA4Calculator calculator; + std::variant nacacode; double te_thickness; Cache wireCache; }; diff --git a/src/wing/NACA5Calculator.cpp b/src/wing/NACA5Calculator.cpp new file mode 100644 index 000000000..51db11990 --- /dev/null +++ b/src/wing/NACA5Calculator.cpp @@ -0,0 +1,325 @@ + +/* +* Copyright (C) 2007-2026 German Aerospace Center (DLR/SC) +* +* Created: 2026-01-14 Hannah Gedler +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +/** +* @file +* @brief Calculator for NACA wing profile coordinates +*/ + +/** +#include "gp_Vec2d.hxx" +#include "tiglmathfunctions.h" +#include "NACA5Calculator.h" +#include "CFunctionToBspline.h" + +namespace tigl{ + + NACA5Calculator::NACA5Calculator(double max_camber_cl, double max_camber_position, double reflex, double max_profile_thickness, double trailing_edge_thickness) + : max_camber_cl(max_camber_cl/100) + , max_camber_position(max_camber_position/20) + , reflex(reflex) + , max_profile_thickness(max_profile_thickness/100) + , trailing_edge_thickness(trailing_edge_thickness/2) + { + + } + + NACA5Calculator::NACA5Calculator(::std::string const& naca_code , const double te_thickness) + : NACA5Calculator(static_cast(naca_code[0] - '0'), + static_cast(naca_code[1] - '0'), + static_cast(naca_code[2] - '0'), + static_cast(std::stoi(naca_code.substr(3,2))), + te_thickness) + {} + + double NACA5Calculator::get_trailing_edge_thickness() const + { + return trailing_edge_thickness; + } + + double NACA5Calculator::k1(double max_camber_cl, double max_camber_position, double reflex) const{ + int max_camber_cl_whole = static_cast(max_camber_cl*100); + int max_camber_position_whole = static_cast(max_camber_position*20); + int reflex = static_cast(reflex); + double k1; + //alles von chapter 8 table 8-6 + meanline_designation = std::to_string(max_camber_cl_whole) + std::to_string(max_camber_position_whole) + std::to_string(reflex); + switch(meanline_designation){ + case "210": + k1 =361.400; + case "220": + k1 =51.640; + case "230": + k1 = 15.957; + case "240": + k1 = 6.643; + case "250": + k1 = 3.230; + case "211": + throw ::std::logic_error("error in NACA5Calculator::k1: this profile does not provide a constant for k1."); //gibts hier nh ander elösung, mit dem ich rechnen kann? + case "221": + k1 = 51.99; + case "231": + k1 = 15.793; + case "241": + k1 = 6.520; + case "251": + k1 = 3.191; + } + return k1; + } + + double NACA5Calculator::m(double max_camber_cl, double max_camber_position, double reflex) const{ + int max_camber_cl_whole = static_cast(max_camber_cl*100); + int max_camber_position_whole = static_cast(max_camber_position*20); + int reflex = static_cast(reflex); + double m; + //alles von chapter 8 table 8-6 + meanline_designation = std::to_string(max_camber_cl_whole) + std::to_string(max_camber_position_whole) + std::to_string(reflex); + switch(meanline_designation){ + case "210": + m = 0.0580; + case "220": + m = 0.1260; + case "230": + m = 0.2025; + case "240": + m = 0.2900; + case "250": + m = 0.3910; + case "211": + throw ::std::logic_error("error in NACA5Calculator::m: this profile does not provide a constant for m."); //gibts hier nh ander elösung, mit dem ich rechnen kann? + case "221": + m = 0.1300; + case "231": + m = 0.2170; + case "241": + m = 0.3180; + case "251": + m = 0.4410; + } + return m; + } +/* + double NACA5Calculator::k2(double max_camber_cl, double max_camber_position, double reflex) const{ + int max_camber_cl_whole = static_cast(max_camber_cl*100); + int max_camber_position_whole = static_cast(max_camber_position*20); + int reflex = static_cast(reflex); + double k2; + //alles von chapter 8 table 8-6 + meanline_designation = std::to_string(max_camber_cl_whole) + std::to_string(max_camber_position_whole) + std::to_string(reflex); + switch(meanline_designation){ + case "210": + throw ::std::logic_error("error in NACA5Calculator::k2: this profile does not provide a constant for k2."); + case "220": + throw ::std::logic_error("error in NACA5Calculator::k2: this profile does not provide a constant for k2."); + case "230": + throw ::std::logic_error("error in NACA5Calculator::k2: this profile does not provide a constant for k2."); + case "240": + throw ::std::logic_error("error in NACA5Calculator::k2: this profile does not provide a constant for k2."); + case "250": + throw ::std::logic_error("error in NACA5Calculator::k2: this profile does not provide a constant for k2."); + case "211": + throw ::std::logic_error("error in NACA5Calculator::k2: this profile does not provide a constant for k2."); + case "221": + k2 = 7.64; + case "231": + k2 = 67.70; + case "241": + k2 = 303.0; + case "251": + k2 = 1355; + } + return k2; + } + */ +/** + double NACA5Calculator::camberline(double x) const{ //geändert + + double s = this->max_camber_cl; + double p = this->max_camber_position; + double q = this->reflex; + double k1 = k1(s, p, q) + double m = m(s, p, q) + double frack1k2 = (3*(m-p)*(m-p)-m*m*m)/((1-m)*(1-m)*(1-m)) + if(q = 0){ + if(p == 0){ + return 0; + } + if(0 <= x && x <= p){ + return (k1/6)(x*x*x-3*m*x*x+m*m*(3-m)*x); + } + else if(x > p){ + return ((k1*m*m*m)/6)*(1-x); + } + } + else if(q == 1){ + if(p == 0){ + return 0; + } + if(0 <= x && x <= p){ + return (k1/6)(((x-m)*(x-m)*(x-m))-(frack1k2*((1-m)*(1-m)*(1-m))*x)- (m*m*m*x+m*m*m)); + } + else if(x > p){ + return (k1/6)((frack1k2)*(x-m)*(x-m)*(x-m)-(frack1k2)*(1-m)*(1-m)*(1-m)*x-m*m*m*x+m*m*m); + } + } + else{ + throw ::std::logic_error("error in NACA5Calculator::camberline: x must be between 0 and 1."); + } + } + + gp_Vec2d NACA5Calculator::upper_curve(double x) const{ //gleich + double yt = profile_thickness(x); + double yc = camberline(x); + auto point = gp_Vec2d{x, yc}; + if(yt < 0){ + throw ::std::logic_error("error in NACA4Calculator::upper_curve: the profile_thickness must be positive or 0."); + } + return point + yt*normal(x); + } + + gp_Vec2d NACA5Calculator::lower_curve(double x) const{ //gleich + double yt = profile_thickness(x); + double yc = camberline(x); + auto point = gp_Vec2d{x, yc}; + if(yt < 0){ + throw ::std::logic_error("error in NACA4Calculator::upper_curve: the profile_thickness must be positive or 0."); + } + return point - yt*normal(x); + } + + double NACA5Calculator::profile_thickness(double x) const{ //bleibt gleich + double t = this->max_profile_thickness; + double e = trailing_edge_thickness_function(trailing_edge_thickness); + return 5*t*(0.2969*sqrt(x) - 0.1260*x - 0.3516*(x*x)+0.2843*pow(x,3) - e*pow(x,4)); + } + + double NACA5Calculator::camberline_derivative(double x) const{ //c'(x) geändert + double s = this->max_camber_cl; + double p = this->max_camber_position; + double q = this->reflex; + double k1 = k1(s, p, q) + double m = m(s, p, q) + double frack1k2 = (3*(m-p)*(m-p)-m*m*m)/((1-m)*(1-m)*(1-m)) + if(p == 0){ + return 0; + } + if(q == 0){ + if(0 <= x && x <= p){ + return (k1/6)(3*x*x-6*m*x+m*m*(3-m)); + } + else if(x > p){ + return -((k1*m*m*m)/6); + } + } + if(q==1){ + if(0 <= x && x <= p){ + return (k1/6)*(3*(x-m)*(x-m)-frack1k2*(1-m)*(1-m)*(1-m)-m*m*m); + } + else if(x > p){ + return (k1/6)(3*frack1k2*(x-m)*(x-m)-frack1k2*(1-m)*(1-m)*(1-m)-m*m*m); + } + } + else{ + throw ::std::logic_error("error in NACA4Calculator::camberline_derivative: x must be between 0 and 1."); + } + } + + gp_Vec2d NACA5Calculator::normal(double x) const{ //bleibt gleich + gp_Vec2d normal = {-camberline_derivative(x), 1}; + double nrm = normal.Magnitude(); + if (nrm < 1e-12) { + return gp_Vec2d{0,1}; + } + return normal/nrm; + } + + double NACA5Calculator::trailing_edge_thickness_function(double y) const{ //bleibt gleich + double t = max_profile_thickness; + if(y < 0){ + throw ::std::logic_error("error in NACA4Calculator::trailing_edge_thickness_function: trailing_edge_thickness must be bigger than 0."); + } + return -(y/(5*t)) + (0.2969 - 0.1260 - 0.3516 + 0.2843); + } + + Handle(Geom_BSplineCurve) NACA5Calculator::upper_bspline() const{ + + NACA4UpperCurve upperCurve(*this); + + const double umin = 0.; + const double umax = 1.; + int degree = 3; + double tolerance=1e-5; + int maxDepth = 10; + + tigl::CFunctionToBspline converter(upperCurve, umin, umax, degree, tolerance, maxDepth); + return converter.Curve(); + } + + Handle(Geom_BSplineCurve) NACA5Calculator::lower_bspline() const{ + NACA4LowerCurve lowerCurve(*this); + + const double umin = 0.; + const double umax = 1.; + int degree = 3; + double tolerance=1e-5; + int maxDepth = 10; + + tigl::CFunctionToBspline converter(lowerCurve, umin, umax, degree, tolerance, maxDepth); + return converter.Curve(); + } + + NACA4UpperCurve::NACA5UpperCurve( NACA4Calculator const& calculator) + : MathFunc3d(), + calculator(calculator) + {} + + double NACA5UpperCurve::valueX(double t) { + gp_Vec2d vec = calculator.upper_curve(t); + return vec.X(); + } + double NACA5UpperCurve::valueY(double t) { + return 0.0; + } + double NACA5UpperCurve::valueZ(double t) { + gp_Vec2d vec = calculator.upper_curve(t); + return vec.Y(); + } + + NACA5LowerCurve::NACA5LowerCurve( NACA5Calculator const& calculator) + : MathFunc3d(), + calculator(calculator) + {} + + double NACA5LowerCurve::valueX(double t) { + gp_Vec2d vec = calculator.lower_curve(t); + return vec.X(); + } + + double NACA5LowerCurve::valueY(double t) { + return 0.0; + } + + double NACA5LowerCurve::valueZ(double t) { + gp_Vec2d vec = calculator.lower_curve(t); + return vec.Y(); + } + +} //namespace tigl +**/ diff --git a/src/wing/NACA5Calculator.h b/src/wing/NACA5Calculator.h new file mode 100644 index 000000000..63b54b2ea --- /dev/null +++ b/src/wing/NACA5Calculator.h @@ -0,0 +1,110 @@ +/* +* Copyright (C) 2007-2026 German Aerospace Center (DLR/SC) +* +* Created: 2026-01-14 Hannah Gedler +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +/** +* @file +* @brief Header for Calculator for NACA wing profile coordinates +*/ + +/* + +#pragma once + +#include +#include +#include "CFunctionToBspline.h" + + +namespace tigl{ + + +class NACA5Calculator{ + + public: + + TIGL_EXPORT NACA5Calculator(double max_camber_cl = 0, double max_camber_position = 0, double reflex = 0,double max_profile_thickness = 12, double trailing_edge_thickness = 0); + + explicit TIGL_EXPORT NACA5Calculator(const ::std::string& naca_code, const double te_thickness); + + + TIGL_EXPORT double get_trailing_edge_thickness() const; + + + TIGL_EXPORT double camberline(double x) const; + + + TIGL_EXPORT gp_Vec2d upper_curve(double x) const; + + + TIGL_EXPORT gp_Vec2d lower_curve(double x) const; + + + TIGL_EXPORT double profile_thickness(double x) const; + + + TIGL_EXPORT double trailing_edge_thickness_function(double y) const; + + TIGL_EXPORT Handle(Geom_BSplineCurve) upper_bspline() const; + + TIGL_EXPORT Handle(Geom_BSplineCurve) lower_bspline() const; + + private: + + double camberline_derivative(double x) const; + + + gp_Vec2d normal(double x) const; + + + private: + double max_camber, max_camber_position, max_profile_thickness, trailing_edge_thickness; +}; + +class NACA4UpperCurve : public MathFunc3d { + public: + TIGL_EXPORT explicit NACA4UpperCurve( NACA4Calculator const& calculator); + + + double valueX(double t) override; + + + + double valueZ(double t) override; + + private: + NACA4Calculator const& calculator; +}; + +class NACA4LowerCurve : public MathFunc3d { + public: + TIGL_EXPORT explicit NACA4LowerCurve( NACA4Calculator const& calculator); + + + double valueX(double t) override; + + + double valueY(double t) override; + + + double valueZ(double t) override; + + private: + NACA4Calculator const& calculator; +}; +} + +*/ \ No newline at end of file diff --git a/tests/TestData/naca_test.cpacs.xml b/tests/TestData/naca_test.cpacs.xml index 3498e4264..3dd8b162e 100644 --- a/tests/TestData/naca_test.cpacs.xml +++ b/tests/TestData/naca_test.cpacs.xml @@ -508,6 +508,7 @@ NACA 4 Series Profile 0012 + 0.1 diff --git a/tests/python/test_naca4profile.py b/tests/python/test_naca4profile.py index 9a529d113..5ef3dd833 100644 --- a/tests/python/test_naca4profile.py +++ b/tests/python/test_naca4profile.py @@ -9,7 +9,8 @@ from tigl3.configuration import ( CCPACSConfigurationManager_get_instance, CPACSNacaProfile, - CTiglNACA4Calculator, + CTiglNACACalculator, + NACA4DigitCode, ) @@ -86,19 +87,14 @@ def test_trailing_edge_point(self) -> None: ) def test_upper_curve_type(self) -> None: - max_camber = 0 - camber_position = 0 - thickness = 9 trailing_edge_thickness = 0.000945 - calculator = CTiglNACA4Calculator( - max_camber, - camber_position, - thickness, + calculator = CTiglNACACalculator( + NACA4DigitCode("0009"), trailing_edge_thickness, ) - self.assertIsInstance(calculator, CTiglNACA4Calculator) + self.assertIsInstance(calculator, CTiglNACACalculator) upper_curve_point = calculator.upper_curve(0.5) diff --git a/tests/unittests/testNACA4Calculator.cpp b/tests/unittests/testNACACalculator.cpp similarity index 59% rename from tests/unittests/testNACA4Calculator.cpp rename to tests/unittests/testNACACalculator.cpp index cb2d8d1c9..c3f874043 100644 --- a/tests/unittests/testNACA4Calculator.cpp +++ b/tests/unittests/testNACACalculator.cpp @@ -20,14 +20,14 @@ * @brief Tests for NACA wing profiles */ -#include "test.h" // Brings in the GTest framework +#include "test.h" #include "testUtils.h" #include "tigl.h" #include "math/tiglmathfunctions.h" #include "CCPACSConfigurationManager.h" #include "CCPACSWingProfile.h" #include "gp_Pnt.hxx" -#include "CTiglNACA4Calculator.h" +#include "CTiglNACACalculator.h" #include "generated/CPACSNacaProfile.h" #include "CTiglUIDManager.h" #include "Debugging.h" @@ -44,13 +44,15 @@ #include "Geom_Curve.hxx" #include "Geom_BSplineCurve.hxx" #include "CTiglError.h" +#include #include -TEST(CTiglNACA4Calculator, naca2212_le_and_te_points){ - tigl::CTiglNACA4Calculator NACA4(2,2,12, 0.00252); + +TEST(CTiglNACACalculator, naca2212_le_and_te_points){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("2212"), 0.00252); ASSERT_TRUE(NACA4.profile_thickness(1) >= 0); gp_Vec2d result1 = NACA4.upper_curve(1); - EXPECT_NEAR(result1.X(), (1.00006), 1e-5); + EXPECT_NEAR(result1.X(), (1.00006), 1e-5); EXPECT_NEAR(result1.Y(), (0.00125843), 1e-8); gp_Vec2d result2 = NACA4.lower_curve(1); EXPECT_NEAR(result2.X(), (0.999937), 1e-6); @@ -65,9 +67,9 @@ TEST(CTiglNACA4Calculator, naca2212_le_and_te_points){ EXPECT_NEAR(result4.Y(), (0.0), 1e-7); } -TEST(CTiglNACA4Calculator, naca2212_le_and_te_points_with_class_lowerCurve){ - tigl::CTiglNACA4Calculator NACA4(2,2,12, 0.00252); - tigl::CTiglNACA4LowerCurve lowerCurve(NACA4); +TEST(CTiglNACACalculator, naca2212_le_and_te_points_with_class_lowerCurve){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("2212"), 0.00252); + tigl::CTiglNACALowerCurve lowerCurve(NACA4); EXPECT_NEAR(lowerCurve.valueX(1), (0.999937), 1e-6); EXPECT_NEAR(lowerCurve.valueY(1), 0.0, 1e-8); EXPECT_NEAR(lowerCurve.valueZ(1), (-0.00125843), 1e-8); @@ -79,8 +81,8 @@ TEST(CTiglNACA4Calculator, naca2212_le_and_te_points_with_class_lowerCurve){ EXPECT_NEAR(result4.Y(), (0.0), 1e-7); } -TEST(CTiglNACA4Calculator, naca0012_random_point){ - tigl::CTiglNACA4Calculator NACA4(0,0,12, 0.00252); +TEST(CTiglNACACalculator, naca0012_at05_and0){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("0012"), 0.00252); ASSERT_TRUE(NACA4.profile_thickness(0.5) >= 0); gp_Vec2d result1 = NACA4.upper_curve(0.5); EXPECT_NEAR(result1.X(), (0.5), 1e-5); @@ -88,7 +90,6 @@ TEST(CTiglNACA4Calculator, naca0012_random_point){ gp_Vec2d result2 = NACA4.lower_curve(0.5); EXPECT_NEAR(result2.X(), (0.5), 1e-6); EXPECT_NEAR(result2.Y(), (-0.0529403), 1e-7); - ASSERT_TRUE(NACA4.profile_thickness(0.) >= 0); gp_Vec2d result3 = NACA4.upper_curve(0.); EXPECT_NEAR(result3.X(), (0.0), 1e-5); @@ -98,8 +99,8 @@ TEST(CTiglNACA4Calculator, naca0012_random_point){ EXPECT_NEAR(result4.Y(), (0.0), 1e-7); } -TEST(CTiglNACA4Calculator, naca0009_random_point_and_le_point){ - tigl::CTiglNACA4Calculator NACA4(0,0,9, 0.00189); +TEST(CTiglNACACalculator, naca0009_random_point_and_le_point){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("0009"), 0.00189); gp_Vec2d result1 = NACA4.upper_curve(0.2); EXPECT_NEAR(result1.X(), (0.2), 1e-5); EXPECT_NEAR(result1.Y(), (0.0430316), 1e-7); @@ -114,8 +115,8 @@ TEST(CTiglNACA4Calculator, naca0009_random_point_and_le_point){ EXPECT_NEAR(result4.Y(), (0.0), 1e-7); } -TEST(CTiglNACA4Calculator, Nnaca6509_le_and_te_points){ - tigl::CTiglNACA4Calculator NACA4(6,5,9, 0.00189); +TEST(CTiglNACACalculator, Nnaca6509_le_and_te_points){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("6509"), 0.00189); gp_Vec2d result1 = NACA4.upper_curve(1); EXPECT_NEAR(result1.X(), (1.00022), 1e-5); EXPECT_NEAR(result1.Y(), (0.000918906), 1e-9); @@ -130,8 +131,8 @@ TEST(CTiglNACA4Calculator, Nnaca6509_le_and_te_points){ EXPECT_NEAR(result4.Y(), (0.0), 1e-7); } -TEST(CTiglNACA4Calculator, naca0012_max_profile_thickness){ - tigl::CTiglNACA4Calculator NACA4(0,0,12, 0.00126); +TEST(CTiglNACACalculator, naca0012_max_profile_thickness){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("0012"), 0.00126); double result1 = NACA4.profile_thickness(0.3); double left_result = NACA4.profile_thickness(0.299); double right_result = NACA4.profile_thickness(0.311); @@ -140,8 +141,8 @@ TEST(CTiglNACA4Calculator, naca0012_max_profile_thickness){ EXPECT_NEAR(result1, 0.06001216339, 1e-11); } -TEST(CTiglNACA4Calculator, naca0018_max_profile_thickness){ - tigl::CTiglNACA4Calculator NACA4(0,0,18, 0.00189); +TEST(CTiglNACACalculator, naca0018_max_profile_thickness){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("0018"), 0.00189); double result1 = NACA4.profile_thickness(0.3); double left_result = NACA4.profile_thickness(0.299); double right_result = NACA4.profile_thickness(0.311); @@ -150,48 +151,45 @@ TEST(CTiglNACA4Calculator, naca0018_max_profile_thickness){ EXPECT_NEAR(result1, 0.09001824509, 1e-12); } -TEST(CTiglNACA4Calculator, naca2212_camberline_at_te_and_le){ - tigl::CTiglNACA4Calculator NACA4(2,2,12, 0.00126); +TEST(CTiglNACACalculator, naca2212_camberline_at_te_and_le){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("2212"), 0.00126); double result1 = NACA4.camberline(0.0); EXPECT_NEAR(result1, 0, 1e-14); double result2 = NACA4.camberline(1.0); EXPECT_NEAR(result1, 0, 1e-14); } -TEST(CTiglNACA4Calculator, naca4509_camberline_at_te_and_le){ - tigl::CTiglNACA4Calculator NACA4(4,5,9, 0.000945); +TEST(CTiglNACACalculator, naca4509_camberline_at_te_and_le){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("4509"), 0.000945); double result1 = NACA4.camberline(0.0); EXPECT_NEAR(result1, 0, 1e-14); double result2 = NACA4.camberline(1.0); EXPECT_NEAR(result1, 0, 1e-14); } -TEST(CTiglNACA4Calculator, naca0015_camberline_at_te_and_le){ - tigl::CTiglNACA4Calculator NACA4(0,0,15, 0.001575); +TEST(CTiglNACACalculator, naca0015_camberline_at_te_and_le){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("0015"), 0.001575); double result1 = NACA4.camberline(0.0); EXPECT_NEAR(result1, 0, 1e-14); double result2 = NACA4.camberline(1.0); EXPECT_NEAR(result1, 0, 1e-14); } -TEST(CTiglNACA4Calculator, naca15030105_assertion_throw_normalization){ - - EXPECT_THROW(tigl::CTiglNACA4Calculator NACA4(150,30,105, 0.001575), tigl::CTiglError); +TEST(CTiglNACACalculator, naca15030105_assertion_throw_normalization){ + EXPECT_THROW(tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("15030105"), 0.001575), tigl::CTiglError); } -TEST(CTiglNACA4Calculator, naca03010_assertion_throw_normalization){ - - EXPECT_THROW(tigl::CTiglNACA4Calculator NACA4(0,30,10, 0.001575), tigl::CTiglError); +TEST(CTiglNACACalculator, naca030810_assertion_throw_normalization){ + EXPECT_THROW(tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("030810"), 0.001575), tigl::CTiglError); } -TEST(CTiglNACA4Calculator, naca03250_assertion_throw_normalization){ - - EXPECT_THROW(tigl::CTiglNACA4Calculator NACA4(0,3,250, 0.001575), tigl::CTiglError); +TEST(CTiglNACACalculator, naca0f3250_assertion_throw_normalization){ + EXPECT_THROW(tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("0f3250"), 0.001575), tigl::CTiglError); } -TEST(CTiglNACA4Calculator, naca0015_trailingedge_length){ - tigl::CTiglNACA4Calculator NACA4(0,0,15, 0.12); +TEST(CTiglNACACalculator, naca0015_trailingedge_length){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("0015"), 0.12); gp_Vec2d result1 = NACA4.upper_curve(1); double half_thickness1_y = result1.Y(); double half_thickness1_x = result1.X(); @@ -199,26 +197,23 @@ TEST(CTiglNACA4Calculator, naca0015_trailingedge_length){ double half_thickness2_y = result2.Y(); double half_thickness2_x = result2.X(); double thickness = sqrt((half_thickness1_x - half_thickness2_x)*(half_thickness1_x - half_thickness2_x)+(half_thickness1_y - half_thickness2_y)*(half_thickness1_y - half_thickness2_y)); - - EXPECT_NEAR(thickness, 0.12, 1e-14); } -TEST(CTiglNACA4Calculator, naca001515_trailingedge_length){ - tigl::CTiglNACA4Calculator NACA4(0,0,15, 0.15); +TEST(CTiglNACACalculator, naca001515_trailingedge_length){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("0015"), 0.15); gp_Vec2d result1 = NACA4.upper_curve(1); double half_thickness1_y = result1.Y(); double half_thickness1_x = result1.X(); gp_Vec2d result2 = NACA4.lower_curve(1); - double half_thickness2_y = result2.Y(); + double half_thickness2_y = result2.Y(); double half_thickness2_x = result2.X(); double thickness = sqrt((half_thickness1_x - half_thickness2_x)*(half_thickness1_x - half_thickness2_x)+(half_thickness1_y - half_thickness2_y)*(half_thickness1_y - half_thickness2_y)); - EXPECT_NEAR(thickness, 0.15, 1e-14); } -TEST(CTiglNACA4Calculator, naca2215_trailingedge_length){ - tigl::CTiglNACA4Calculator NACA4(2,2,15, 0.20); +TEST(CTiglNACACalculator, naca2215_trailingedge_length){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("2215"), 0.20); gp_Vec2d result1 = NACA4.upper_curve(1); double half_thickness1_y = result1.Y(); double half_thickness1_x = result1.X(); @@ -226,13 +221,11 @@ TEST(CTiglNACA4Calculator, naca2215_trailingedge_length){ double half_thickness2_y = result2.Y(); double half_thickness2_x = result2.X(); double thickness = sqrt((half_thickness1_x - half_thickness2_x)*(half_thickness1_x - half_thickness2_x)+(half_thickness1_y - half_thickness2_y)*(half_thickness1_y - half_thickness2_y)); - - EXPECT_NEAR(thickness, 0.20, 1e-14); } -TEST(CTiglNACA4Calculator, naca6415_trailingedge_length){ - tigl::CTiglNACA4Calculator NACA4(6,4,15, 0.13); +TEST(CTiglNACACalculator, naca6415_trailingedge_length){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("6415"), 0.13); gp_Vec2d result1 = NACA4.upper_curve(1); double half_thickness1_y = result1.Y(); double half_thickness1_x = result1.X(); @@ -240,14 +233,11 @@ TEST(CTiglNACA4Calculator, naca6415_trailingedge_length){ double half_thickness2_y = result2.Y(); double half_thickness2_x = result2.X(); double thickness = sqrt((half_thickness1_x - half_thickness2_x)*(half_thickness1_x - half_thickness2_x)+(half_thickness1_y - half_thickness2_y)*(half_thickness1_y - half_thickness2_y)); - - EXPECT_NEAR(thickness, 0.13, 1e-14); } - namespace { // Mirrors leParam() in CTiglNACA4Calculator.cpp: the leading-edge reparametrization - // x(t) = (1+eps)*t*t/(t+eps) used by CTiglNACA4UpperCurve/LowerCurve::valueX/valueZ, so + // x(t) = (1+eps)*t*t/(t+eps) used by CTiglNACAUpperCurve/LowerCurve::valueX/valueZ, so // these tests can independently compute the expected reparametrized chord fraction. double testLeParam(double t) { @@ -256,50 +246,31 @@ namespace { } } -TEST(CTiglNACA4Calculator, naca2212_upperCurve_ycoord_and_upper_curve_x_and_zcoord){ - tigl::CTiglNACA4Calculator NACA4(2,2,12, 15); - tigl::CTiglNACA4UpperCurve upperCurve(NACA4); +TEST(CTiglNACACalculator, naca2212_upperCurve_ycoord_and_upper_curve_x_and_zcoord){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("2212"), 15); + tigl::CTiglNACAUpperCurve upperCurve(NACA4); ASSERT_EQ(upperCurve.valueY(0.), 0.); ASSERT_EQ(upperCurve.valueY(0.5), 0.); ASSERT_EQ(upperCurve.valueY(1.), 0.); - // upperCurve.valueX/valueZ(t) evaluate the analytic curve at a reparametrized chord - // fraction x(t) = testLeParam(t), not x=t directly (see CTiglNACA4UpperCurve::valueX) gp_Vec2d pnt = NACA4.upper_curve(testLeParam(0.5)); ASSERT_EQ(upperCurve.valueX(0.5), pnt.X()); ASSERT_EQ(upperCurve.valueZ(0.5), pnt.Y()); } -TEST(CTiglNACA4Calculator, naca2212_lowerCurve_ycoord_and_lower_curve_x_and_zcoord){ - tigl::CTiglNACA4Calculator NACA4(2,2,12, 15); - tigl::CTiglNACA4LowerCurve lowerCurve(NACA4); - +TEST(CTiglNACACalculator, naca2212_lowerCurve_ycoord_and_lower_curve_x_and_zcoord){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("2212"), 15); + tigl::CTiglNACALowerCurve lowerCurve(NACA4); ASSERT_EQ(lowerCurve.valueY(0.), 0.); ASSERT_EQ(lowerCurve.valueY(0.5), 0.); ASSERT_EQ(lowerCurve.valueY(1.), 0.); - - // lowerCurve.valueX/valueZ(t) evaluate the analytic curve at a reparametrized chord - // fraction x(t) = testLeParam(t), not x=t directly (see CTiglNACA4UpperCurve::valueX) gp_Vec2d pnt = NACA4.lower_curve(testLeParam(0.5)); ASSERT_EQ(lowerCurve.valueX(0.5), pnt.X()); ASSERT_EQ(lowerCurve.valueZ(0.5), pnt.Y()); } -TEST(CTiglNACA4Calculator, naca2212_bspline_vs_lower_curve_coord) -{ - tigl::CTiglNACA4Calculator NACA4(2,2,12, 15); - Handle(Geom_BSplineCurve) lower_spline = NACA4.lower_bspline(); - - // the bspline's own parameter u corresponds to x=testLeParam(u), not x=u directly (see - // CTiglNACA4UpperCurve::valueX) - gp_Vec2d pnt = NACA4.lower_curve(testLeParam(0.5)); - gp_Pnt pnt2; - lower_spline->D0(0.5, pnt2); - ASSERT_NEAR(pnt2.X(), pnt.X(), 1e-4); - ASSERT_NEAR(pnt2.Z(), pnt.Y(), 1e-4); -} -TEST(CTiglNACA4Calculator, naca2212_export_bsplines){ - tigl::CTiglNACA4Calculator NACA4(2,2,12, .15); +TEST(CTiglNACACalculator, naca2212_assert_and_export_bsplines){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("2212"), .15); Handle(Geom_BSplineCurve) upperCurve = NACA4.upper_bspline(); Handle(Geom_BSplineCurve) lowerCurve = NACA4.lower_bspline(); ASSERT_FALSE(upperCurve.IsNull()); @@ -307,14 +278,13 @@ TEST(CTiglNACA4Calculator, naca2212_export_bsplines){ auto lowerEdge = BRepBuilderAPI_MakeEdge(lowerCurve).Edge(); ASSERT_FALSE(lowerEdge.IsNull()); BRepTools::Write(lowerEdge, "TestData/export/lowerEdgeTest.brep"); - auto upperEdge = BRepBuilderAPI_MakeEdge(upperCurve).Edge(); ASSERT_FALSE(upperEdge.IsNull()); BRepTools::Write(upperEdge, "TestData/export/upperEdgeTest.brep"); } -TEST(CTiglNACA4Calculator, naca0012_export_bsplines){ - tigl::CTiglNACA4Calculator NACA4(0,0,12, .015); +TEST(CTiglNACACalculator, naca0012_export_bsplines){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("0012"), .015); Handle(Geom_BSplineCurve) upperCurvesym = NACA4.upper_bspline(); Handle(Geom_BSplineCurve) lowerCurvesym = NACA4.lower_bspline(); ASSERT_FALSE(upperCurvesym.IsNull()); @@ -322,32 +292,24 @@ TEST(CTiglNACA4Calculator, naca0012_export_bsplines){ auto lowerEdge = BRepBuilderAPI_MakeEdge(lowerCurvesym).Edge(); ASSERT_FALSE(lowerEdge.IsNull()); BRepTools::Write(lowerEdge, "TestData/export/lowerEdgeTest_symetric.brep"); - auto upperEdge = BRepBuilderAPI_MakeEdge(upperCurvesym).Edge(); ASSERT_FALSE(upperEdge.IsNull()); BRepTools::Write(upperEdge, "TestData/export/upperEdgeTest_symetric.brep"); } -TEST(CTiglNACA4Calculator, naca2412_LePoint_TePoint){ - - // Create a UID manager and use it for the profile. Parent container is not needed for this unit test. +TEST(CTiglNACACalculator, naca2412_LePoint_TePoint){ tigl::CTiglUIDManager uidMgr; - // disambiguate nullptr for overloaded constructors by casting to the intended parent type tigl::CCPACSWingProfile cpacsProfile(static_cast(nullptr), &uidMgr); tigl::generated::CPACSNacaProfile nacadef(&cpacsProfile); - nacadef.SetNaca4DigitCode_choice1(boost::optional(std::string("2412"))); nacadef.SetTrailingEdgeThickness(boost::optional(0.15)); - tigl::CTiglWingProfileNACA profile(cpacsProfile, nacadef); - TopoDS_Edge upper = profile.GetUpperWire(); TopoDS_Edge lower = profile.GetLowerWire(); TopoDS_Edge te = profile.GetTrailingEdge(); gp_Pnt lePoint = profile.GetLEPoint(); gp_Pnt tePoint = profile.GetTEPoint(); - - tigl::CTiglNACA4Calculator NACA4(2,4,12, 0.15); + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("2412"), 0.15); gp_Pnt result1 = profile.GetLEPoint(); double result2 = NACA4.camberline(0.0); gp_Vec2d result3 = NACA4.upper_curve(0.0); @@ -355,40 +317,30 @@ TEST(CTiglNACA4Calculator, naca2412_LePoint_TePoint){ EXPECT_EQ(result2, result1.Z()); EXPECT_EQ(result1.Z(), result3.Y()); EXPECT_EQ(result1.Z(), result4.Y()); - EXPECT_FALSE(upper.IsNull()); EXPECT_FALSE(lower.IsNull()); EXPECT_FALSE(te.IsNull()); - TopoDS_Edge ul = profile.GetUpperLowerWire(); EXPECT_FALSE(ul.IsNull()); TopoDS_Edge ulSharp = profile.GetUpperLowerWire(SHARP_TRAILINGEDGE); EXPECT_FALSE(ulSharp.IsNull()); TopoDS_Edge ulBlunt = profile.GetUpperLowerWire(BLUNT_TRAILINGEDGE); EXPECT_FALSE(ulBlunt.IsNull()); - } -TEST(CTiglNACA4Calculator, naca0012_LePoint_TePoint){ - - // Create a UID manager and use it for the profile. Parent container is not needed for this unit test. +TEST(CTiglNACACalculator, naca0012_LePoint_TePoint){ tigl::CTiglUIDManager uidMgr; - // disambiguate nullptr for overloaded constructors by casting to the intended parent type tigl::CCPACSWingProfile cpacsProfile(static_cast(nullptr), &uidMgr); tigl::generated::CPACSNacaProfile nacadef(&cpacsProfile); - nacadef.SetNaca4DigitCode_choice1(boost::optional(std::string("0012"))); nacadef.SetTrailingEdgeThickness(boost::optional(0.015)); - tigl::CTiglWingProfileNACA profile(cpacsProfile, nacadef); - TopoDS_Edge upper = profile.GetUpperWire(); TopoDS_Edge lower = profile.GetLowerWire(); TopoDS_Edge te = profile.GetTrailingEdge(); gp_Pnt lePoint = profile.GetLEPoint(); gp_Pnt tePoint = profile.GetTEPoint(); - - tigl::CTiglNACA4Calculator NACA4(0,0,12, 0.015); + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("0012"), 0.015); gp_Pnt result1 = profile.GetLEPoint(); double result2 = NACA4.camberline(0.0); gp_Vec2d result3 = NACA4.upper_curve(0.0); @@ -396,36 +348,30 @@ TEST(CTiglNACA4Calculator, naca0012_LePoint_TePoint){ EXPECT_EQ(result2, result1.Z()); EXPECT_EQ(result1.Z(), result3.Y()); EXPECT_EQ(result1.Z(), result4.Y()); - EXPECT_FALSE(upper.IsNull()); EXPECT_FALSE(lower.IsNull()); EXPECT_FALSE(te.IsNull()); - TopoDS_Edge ul = profile.GetUpperLowerWire(); EXPECT_FALSE(ul.IsNull()); TopoDS_Edge ulSharp = profile.GetUpperLowerWire(SHARP_TRAILINGEDGE); EXPECT_FALSE(ulSharp.IsNull()); TopoDS_Edge ulBlunt = profile.GetUpperLowerWire(BLUNT_TRAILINGEDGE); EXPECT_FALSE(ulBlunt.IsNull()); - } -TEST(CTiglNACA4Calculator, naca0012_trailingEdge_absent_when_zero_thickness){ +TEST(CTiglNACACalculator, naca0012_trailingEdge_absent_when_zero_thickness){ tigl::CTiglUIDManager uidMgr; tigl::CCPACSWingProfile cpacsProfile(static_cast(nullptr), &uidMgr); tigl::generated::CPACSNacaProfile nacadef(&cpacsProfile); const std::string code = "0012"; nacadef.SetNaca4DigitCode_choice1(boost::optional(code)); nacadef.SetTrailingEdgeThickness(boost::optional(0.0)); - tigl::CTiglWingProfileNACA profile(cpacsProfile, nacadef); - TopoDS_Edge te = profile.GetTrailingEdge(); EXPECT_TRUE(te.IsNull()); } -TEST(CTiglNACA4Calculator, naca2412_edge_counter){ - +TEST(CTiglNACACalculator, naca2412_edge_counter){ tigl::CTiglUIDManager uidMgr; tigl::CCPACSWingProfile cpacsProfile(static_cast(nullptr), &uidMgr); tigl::generated::CPACSNacaProfile nacadef(&cpacsProfile); @@ -433,26 +379,17 @@ TEST(CTiglNACA4Calculator, naca2412_edge_counter){ const double te_thickness = 0.15; nacadef.SetNaca4DigitCode_choice1(boost::optional(code)); nacadef.SetTrailingEdgeThickness(boost::optional(te_thickness)); - tigl::CTiglWingProfileNACA profile(cpacsProfile, nacadef); - - TopoDS_Edge upper = profile.GetUpperWire(); TopoDS_Edge lower = profile.GetLowerWire(); TopoDS_Edge te = profile.GetTrailingEdge(); - - ASSERT_TRUE(BRepCheck_Analyzer(upper).IsValid()); ASSERT_TRUE(BRepCheck_Analyzer(lower).IsValid()); - - - tigl::CTiglNACA4Calculator calc(code, te_thickness); + tigl::CTiglNACACalculator calc(tigl::NACA4DigitCode(code), te_thickness); EXPECT_EQ(profile.HasBluntTE(), calc.get_trailing_edge_thickness() > 0.0); - if (!te.IsNull()) { ASSERT_TRUE(BRepCheck_Analyzer(te).IsValid()); } - BRepBuilderAPI_MakeWire closedWireBuilder; closedWireBuilder.Add(lower); closedWireBuilder.Add(upper); @@ -463,7 +400,6 @@ TEST(CTiglNACA4Calculator, naca2412_edge_counter){ ASSERT_TRUE(closedWireBuilder.IsDone()); TopoDS_Wire wire = closedWireBuilder.Wire(); ASSERT_TRUE(BRepCheck_Analyzer(wire).IsValid()); - int edgeCount = 0; for (TopExp_Explorer ex(wire, TopAbs_EDGE); ex.More(); ex.Next()) ++edgeCount; if (te.IsNull()) { @@ -473,77 +409,181 @@ TEST(CTiglNACA4Calculator, naca2412_edge_counter){ } } -TEST(CTiglNACA4Calculator, upper_curve_does_not_throw_for_valid_x_range){ - tigl::CTiglNACA4Calculator naca(2, 4, 12, 0.00252); - for (double x = 0.0; x <= 1.0; x += 0.05) { - EXPECT_NO_THROW((void)naca.upper_curve(x)); - } +TEST(CTiglNACACalculator, naca22012_le_and_te_points){ + tigl::CTiglNACACalculator NACA5(tigl::NACA5DigitCode("22112"), 0.00252); + gp_Vec2d result1 = NACA5.upper_curve(1); + EXPECT_NEAR(result1.Y(), (0.00125923), 1e-6); } +TEST(CTiglNACACalculator, naca23012_export_bsplines){ + tigl::CTiglNACACalculator NACA5(tigl::NACA5DigitCode("22018"), 0.0); + Handle(Geom_BSplineCurve) upperCurve = NACA5.upper_bspline(); + Handle(Geom_BSplineCurve) lowerCurve = NACA5.lower_bspline(); + ASSERT_FALSE(upperCurve.IsNull()); + ASSERT_FALSE(lowerCurve.IsNull()); + auto lowerEdge = BRepBuilderAPI_MakeEdge(lowerCurve).Edge(); + ASSERT_FALSE(lowerEdge.IsNull()); + BRepTools::Write(lowerEdge, "TestData/export/lowerEdgeTest5_22012.brep"); + + auto upperEdge = BRepBuilderAPI_MakeEdge(upperCurve).Edge(); + ASSERT_FALSE(upperEdge.IsNull()); + BRepTools::Write(upperEdge, "TestData/export/upperEdgeTest5_22012.brep"); +} +TEST(CTiglNACACalculator, naca23012_python){ + tigl::CTiglNACACalculator NACA5(tigl::NACA5DigitCode("22018"), 0.0); + Handle(Geom_BSplineCurve) upperCurve = NACA5.upper_bspline(); + Handle(Geom_BSplineCurve) lowerCurve = NACA5.lower_bspline(); + ASSERT_FALSE(upperCurve.IsNull()); + ASSERT_FALSE(lowerCurve.IsNull()); + auto lowerEdge = BRepBuilderAPI_MakeEdge(lowerCurve).Edge(); + ASSERT_FALSE(lowerEdge.IsNull()); + BRepTools::Write(lowerEdge, "TestData/export/lowerEdgeTest5_22012.brep"); + auto upperEdge = BRepBuilderAPI_MakeEdge(upperCurve).Edge(); + ASSERT_FALSE(upperEdge.IsNull()); + BRepTools::Write(upperEdge, "TestData/export/upperEdgeTest5_22012.brep"); +} -TEST(CTiglNACA4Calculator, naca2412_getUpperLowerWire) { +TEST(CTiglNACACalculator, naca22112_le_and_te_points_with_class_lowerCurve){ + tigl::CTiglNACACalculator NACA5(tigl::NACA5DigitCode("22112"), 0.00252); + tigl::CTiglNACALowerCurve lowerCurve(NACA5); + EXPECT_NEAR(lowerCurve.valueX(1), (0.99999099), 1e-5); + EXPECT_NEAR(lowerCurve.valueY(1), 0.0, 1e-8); + EXPECT_NEAR(lowerCurve.valueZ(1), (-0.00125997), 1e-7);/ +} + +TEST(CTiglNACACalculator, naca24112_le_and_te_points_with_class_lowerCurve){ + tigl::CTiglNACACalculator NACA5(tigl::NACA5DigitCode("24112"), 0.00252); + tigl::CTiglNACALowerCurve lowerCurve(NACA5); + EXPECT_NEAR(lowerCurve.valueX(1), (0.99999999), 1e-6); + EXPECT_NEAR(lowerCurve.valueY(1), 0.0, 1e-8); + EXPECT_NEAR(lowerCurve.valueZ(1), (-0.00126000), 1e-6); +} + +TEST(CTiglNACACalculator, naca22112_export_bsplines){ + tigl::CTiglNACACalculator NACA5(tigl::NACA5DigitCode("22112"), 0.00252); + Handle(Geom_BSplineCurve) upperCurve = NACA5.upper_bspline(); + Handle(Geom_BSplineCurve) lowerCurve = NACA5.lower_bspline(); + ASSERT_FALSE(upperCurve.IsNull()); + ASSERT_FALSE(lowerCurve.IsNull()); + auto lowerEdge = BRepBuilderAPI_MakeEdge(lowerCurve).Edge(); + ASSERT_FALSE(lowerEdge.IsNull()); + BRepTools::Write(lowerEdge, "TestData/export/lowerEdgeTest5.brep"); + auto upperEdge = BRepBuilderAPI_MakeEdge(upperCurve).Edge(); + ASSERT_FALSE(upperEdge.IsNull()); + BRepTools::Write(upperEdge, "TestData/export/upperEdgeTest5.brep"); +} + +TEST(CTiglNACACalculator, naca24112_export_bsplines){ + tigl::CTiglNACACalculator NACA5(tigl::NACA5DigitCode("24112"), 0.00252); + Handle(Geom_BSplineCurve) upperCurve = NACA5.upper_bspline(); + Handle(Geom_BSplineCurve) lowerCurve = NACA5.lower_bspline(); + ASSERT_FALSE(upperCurve.IsNull()); + ASSERT_FALSE(lowerCurve.IsNull()); + auto lowerEdge = BRepBuilderAPI_MakeEdge(lowerCurve).Edge(); + ASSERT_FALSE(lowerEdge.IsNull()); + BRepTools::Write(lowerEdge, "TestData/export/lowerEdgeTest5_24112.brep"); + auto upperEdge = BRepBuilderAPI_MakeEdge(upperCurve).Edge(); + ASSERT_FALSE(upperEdge.IsNull()); + BRepTools::Write(upperEdge, "TestData/export/upperEdgeTest5_24112.brep"); +} + +TEST(CTiglNACACalculator, naca24112_upper_curve_does_not_throw_for_valid_x_range){ + tigl::CTiglNACACalculator NACA5(tigl::NACA5DigitCode("24112"), 0.00252); + for(double x =0; x<1; x+=0.05){ + auto z = NACA5.upper_curve(x); + } +} + +TEST(CTiglNACACalculator, upper_curve_does_not_throw_for_valid_x_range){ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("2412"), 0.00252); + for (double x = 0.0; x <= 1.0; x += 0.05) { + EXPECT_NO_THROW((void)NACA4.upper_curve(x)); + } +} + +TEST(CTiglNACACalculator, naca2412_getUpperLowerWire) { tigl::CTiglUIDManager uidMgr; tigl::CCPACSWingProfile cpacsProfile(static_cast(nullptr), &uidMgr); tigl::generated::CPACSNacaProfile nacadef(&cpacsProfile); - nacadef.SetNaca4DigitCode_choice1(boost::optional(std::string("2412"))); nacadef.SetTrailingEdgeThickness(boost::optional(0.15)); - tigl::CTiglWingProfileNACA profile(cpacsProfile, nacadef); - TopoDS_Edge ul = profile.GetUpperLowerWire(); EXPECT_FALSE(ul.IsNull()); - TopoDS_Edge ulSharp = profile.GetUpperLowerWire(SHARP_TRAILINGEDGE); EXPECT_FALSE(ulSharp.IsNull()); - TopoDS_Edge ulBlunt = profile.GetUpperLowerWire(BLUNT_TRAILINGEDGE); EXPECT_FALSE(ulBlunt.IsNull()); - EXPECT_TRUE(profile.HasBluntTE()); + Standard_Real u1, u2; + Handle(Geom_Curve) ulCurve = BRep_Tool::Curve(ul, u1, u2); + ASSERT_FALSE(ulCurve.IsNull()); +} +TEST(CTiglNACACalculator, naca23012_getUpperLowerWire) { + tigl::CTiglUIDManager uidMgr; + tigl::CCPACSWingProfile cpacsProfile(static_cast(nullptr), &uidMgr); + tigl::generated::CPACSNacaProfile nacadef(&cpacsProfile); + nacadef.SetNaca5DigitCode_choice2(boost::optional(std::string("23012"))); + nacadef.SetTrailingEdgeThickness(boost::optional(0.15)); + tigl::CTiglWingProfileNACA profile(cpacsProfile, nacadef); + TopoDS_Edge ul = profile.GetUpperLowerWire(); + EXPECT_FALSE(ul.IsNull()); + TopoDS_Edge ulSharp = profile.GetUpperLowerWire(SHARP_TRAILINGEDGE); + EXPECT_FALSE(ulSharp.IsNull()); + TopoDS_Edge ulBlunt = profile.GetUpperLowerWire(BLUNT_TRAILINGEDGE); + EXPECT_FALSE(ulBlunt.IsNull()); + EXPECT_TRUE(profile.HasBluntTE()); Standard_Real u1, u2; Handle(Geom_Curve) ulCurve = BRep_Tool::Curve(ul, u1, u2); ASSERT_FALSE(ulCurve.IsNull()); } -// Regression test for CFunctionToBspline::concatC1 (used internally by -// CTiglNACA4Calculator::upper_bspline/lower_bspline). A strongly cambered profile forces the -// adaptive Chebyshev fit to produce many segments, exercising the multi-segment -// concatenation path. Away from the leading/trailing edge (where the thickness formula's -// sqrt(x) term makes the true tangent direction change very fast / become vertical - an -// inherent feature of the NACA4 shape, not a bug), every internal knot join should now be -// honestly continuous, since concatC1 checks continuity instead of blindly declaring it. -TEST(CTiglNACA4Calculator, naca6415_upper_lower_bspline_c1_continuous_everywhere) -{ - // Strongly cambered profile: exercises the adaptive multi-segment path in - // CFunctionToBspline. CTiglNACA4UpperCurve/LowerCurve reparametrize near the leading - // edge (see leParam in CTiglNACA4Calculator.cpp), which removes the thickness - // distribution's sqrt(x) derivative singularity there - so unlike before that - // reparametrization, every internal knot (including right at the leading/trailing edge) - // should now be honestly C1 continuous, with no margin needed to exclude a "genuinely - // near-vertical" region. - tigl::CTiglNACA4Calculator NACA4(6, 4, 15, 0.13); +TEST(CTiglNACACalculator, naca5digit_22018_bsplinePoles_fromXML) { + const char* xmlfile = "/localdata2/gedl_ha/code/tigl/tests/TestData/naca_5_test.xml"; + TixiDocumentHandle tixiHandle = -1; + TiglCPACSConfigurationHandle tiglHandle = -1; + EXPECT_EQ(tixiOpenDocument(xmlfile, &tixiHandle), SUCCESS); + ASSERT_TRUE(tixiHandle >= 0); + EXPECT_EQ(tiglOpenCPACSConfiguration(tixiHandle, "", &tiglHandle), TIGL_SUCCESS); + char* nacaCodeStr = NULL; + std::string xpath = "//wingAirfoil[@uID='NACA0009']/nacaProfile/naca5DigitCode"; + EXPECT_EQ(tixiGetTextElement(tixiHandle, xpath.c_str(), &nacaCodeStr), SUCCESS); + ASSERT_STREQ(nacaCodeStr, "22018"); + tigl::CTiglNACACalculator nacaCalc(tigl::NACA5DigitCode(nacaCodeStr), 0.0); + Handle(Geom_BSplineCurve) upperCurve = nacaCalc.upper_bspline(); + Handle(Geom_BSplineCurve) lowerCurve = nacaCalc.lower_bspline(); + ASSERT_FALSE(upperCurve.IsNull()); + ASSERT_FALSE(lowerCurve.IsNull()); + auto upperEdge = BRepBuilderAPI_MakeEdge(upperCurve).Edge(); + auto lowerEdge = BRepBuilderAPI_MakeEdge(lowerCurve).Edge(); + ASSERT_FALSE(upperEdge.IsNull()); + ASSERT_FALSE(lowerEdge.IsNull()); + BRepTools::Write(upperEdge, "TestData/export/upperEdgeTest5_22018_fromXML.brep"); + BRepTools::Write(lowerEdge, "TestData/export/lowerEdgeTest5_22018_fromXML.brep"); + EXPECT_EQ(tiglCloseCPACSConfiguration(tiglHandle), TIGL_SUCCESS); + EXPECT_EQ(tixiCloseDocument(tixiHandle), SUCCESS); +} + + +TEST(CTiglNACACalculator, naca6415_upper_lower_bspline_c1_continuous_everywhere) +{ + tigl::CTiglNACACalculator NACA4(tigl::NACA4DigitCode("6415"), 0.13); auto checkContinuity = [](const Handle(Geom_BSplineCurve)& curve) { - // make sure this actually exercises the multi-segment concatenation path ASSERT_GT(curve->NbKnots(), 2); - const double eps = 1e-7; for (int i = 2; i < curve->NbKnots(); ++i) { double u = curve->Knot(i); - gp_Pnt pLeft, pRight; gp_Vec dLeft, dRight; curve->D1(u - eps, pLeft, dLeft); curve->D1(u + eps, pRight, dRight); - EXPECT_NEAR(pLeft.Distance(pRight), 0.0, 1e-6) << "position jump at knot " << u; EXPECT_NEAR(dLeft.Angle(dRight), 0.0, 1e-3) << "tangent kink at knot " << u; } }; - checkContinuity(NACA4.upper_bspline()); checkContinuity(NACA4.lower_bspline()); }