From 9f283000b78d0aa389287d7e92c495a333c4d32e Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Mon, 11 Feb 2013 18:35:49 +0000
Subject: [PATCH] Added NACA00 level set
---
Geo/gmshLevelset.cpp | 40 ++++++++++++++++++++++++++++++++++++++++
Geo/gmshLevelset.h | 10 ++++++++++
2 files changed, 50 insertions(+)
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index a7dbd96eaf..5fdadd4934 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -1157,3 +1157,43 @@ gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1,
gLevelsetConrod::gLevelsetConrod(const gLevelsetConrod &lv)
: gLevelsetImproved(lv){}
+
+
+double gLevelsetNACA00::operator() (const double x, const double y, const double z) const
+{
+
+ static const int maxIter = 100;
+ static const double tol = 1.e-8;
+
+ const double tolr = tol/_c; // Tolerance (scaled bu chord)
+ double distSq; // Square of dist. between point and boundary
+ bool in = false; // Whether the point is inside the airfoil
+
+ const double xt = x-_x0, yt = fabs(y-_y0); // Point translated according to airfoil origin and symmetry
+
+ if (xt-_c > 1.16925*_t*yt) // Behind line normal to airfoil at trailing edge, closest
+ distSq = (xt-_c)*(xt-_c)+yt*yt; // boundary point is trailing edge
+ else { // Otherwise Newton-Raphson to find closest boundary point
+ const double fact = 5.*_t*_c;
+ double xb = std::max(xt,tolr);
+ for (int it=0; it<maxIter; it++) {
+ const double xbr = xb/_c, sxbr = sqrt(xbr), xbr32 = xbr*sxbr,
+ xbr2 = xbr*xbr, xbr3 = xbr2*xbr, xbr4 = xbr2*xbr2;
+ double yb = fact*(0.2969*sxbr-0.1260*xbr-0.3516*xbr2+0.2843*xbr3-0.1015*xbr4);
+ const double dyb = fact*(0.14845/sxbr-0.406*xbr3+0.8529*xbr2-0.7032*xbr-0.126)/_c;
+ const double ddyb = fact*(-0.074225/xbr32-1.218*xbr2+1.7058*xbr-0.7032)/(_c*_c);
+ const double xx = xt-xb, yy = yt-yb;
+ in = (xt > 0) && (yy < 0);
+ distSq = xx*xx+yy*yy;
+ const double dDistSq = -2.*(xx+dyb*yy);
+ const double ddDistSq = 2.*(1.-ddyb*yy+dyb*dyb);
+ const double mIncr = dDistSq/ddDistSq;
+ if (fabs(mIncr) < tolr) break;
+ else xb -= mIncr;
+ if (xb < tolr) xb = tolr;
+ }
+ }
+
+ return in ? -sqrt(distSq) : sqrt(distSq);
+
+}
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 476992e55f..d4cb0f2d33 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -600,4 +600,14 @@ public:
/* int type2() const {return DISK;} */
/* }; */
+class gLevelsetNACA00: public gLevelsetPrimitive
+{
+ double _x0, _y0, _c, _t;
+public:
+ gLevelsetNACA00(double x0, double y0, double c, double t) : _x0(x0), _y0(y0), _c(c), _t(t) {}
+ ~gLevelsetNACA00() {}
+ double operator () (const double x, const double y, const double z) const;
+ int type() const { return UNKNOWN; }
+};
+
#endif
--
GitLab