From d289aae9e11645e5cb83cd5610c091450a33a77b Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Tue, 23 Apr 2013 00:18:25 +0000
Subject: [PATCH] Added analytical gradients and Hessian for NACA00 level-set
---
Geo/gmshLevelset.cpp | 103 ++++++++++++++++++++++++++++++++++++-------
Geo/gmshLevelset.h | 12 ++++-
2 files changed, 99 insertions(+), 16 deletions(-)
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 6d82c1cd2c..97be674aee 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -1229,41 +1229,114 @@ gLevelsetConrod::gLevelsetConrod(const gLevelsetConrod &lv)
// Level-set for NACA0012 airfoil, last coeff. modified for zero-thickness trailing edge
// cf. http://en.wikipedia.org/wiki/NACA_airfoil
-double gLevelsetNACA00::operator() (const double x, const double y, const double z) const
+gLevelsetNACA00::gLevelsetNACA00(double x0, double y0, double c, double t)
+ : _x0(x0), _y0(y0), _c(c), _t(t)
+{
+ _hasDerivatives = true;
+}
+
+
+void gLevelsetNACA00::getClosestBndPoint(const double x, const double y, const double z,
+ double &xb, double &yb, double &curvRad, bool &in) 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
+ 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.21125*_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
+ if (xt-_c > 1.21125*_t*yt) { // Behind line normal to airfoil at trailing edge, closest
+ xb = _x0+_c; // boundary point is trailing edge...
+ yb = _y0;
+ curvRad = 0.;
+ }
+ else { // ...otherwise Newton-Raphson to find closest boundary point
const double fact = 5.*_t*_c;
- double xb = std::max(xt,tolr);
+ double xtb = std::max(xt,tolr), ytb;
+ double dyb, ddyb;
for (int it=0; it<maxIter; it++) {
- const double xbr = xb/_c, sxbr = sqrt(xbr), xbr32 = xbr*sxbr,
+ const double xbr = xtb/_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.1036*xbr4);
- const double dyb = fact*(0.14845/sxbr-0.4144*xbr3+0.8529*xbr2-0.7032*xbr-0.126)/_c;
- const double ddyb = fact*(-0.074225/xbr32-1.2432*xbr2+1.7058*xbr-0.7032)/(_c*_c);
- const double xx = xt-xb, yy = yt-yb;
+ ytb = fact*(0.2969*sxbr-0.1260*xbr-0.3516*xbr2+0.2843*xbr3-0.1036*xbr4);
+ dyb = fact*(0.14845/sxbr-0.4144*xbr3+0.8529*xbr2-0.7032*xbr-0.126)/_c;
+ ddyb = fact*(-0.074225/xbr32-1.2432*xbr2+1.7058*xbr-0.7032)/(_c*_c);
+ const double xx = xt-xtb, yy = yt-ytb;
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;
+ else xtb -= mIncr;
+ if (xtb < tolr) xtb = tolr;
}
+ xb = _x0+xtb;
+ yb = (y >= _y0) ? _y0+ytb : -(_y0+ytb);
+ const double norm = sqrt(1.+dyb*dyb);
+ curvRad = norm*norm*norm/fabs(ddyb);
}
+}
+
+
+double gLevelsetNACA00::operator() (const double x, const double y, const double z) const
+{
+
+ double xb, yb, curvRadb;
+ bool in;
+
+ getClosestBndPoint(x, y, z, xb, yb, curvRadb, in);
+ const double xx = x-xb, yy = y-yb, distSq = xx*xx+yy*yy;
+
return in ? -sqrt(distSq) : sqrt(distSq);
}
+
+
+void gLevelsetNACA00::gradient (double x, double y, double z,
+ double & dfdx, double & dfdy, double & dfdz) const
+{
+
+ double xb, yb, curvRadb;
+ bool in;
+
+ getClosestBndPoint(x, y, z, xb, yb, curvRadb, in);
+ const double xx = x-xb, yy = y-yb, distSq = xx*xx+yy*yy;
+ const double dist = in ? -sqrt(distSq) : sqrt(distSq);
+
+ dfdx = xx/dist;
+ dfdy = yy/dist;
+ dfdz = 0.;
+
+}
+
+
+void gLevelsetNACA00::hessian (double x, double y, double z,
+ double & dfdxx, double & dfdxy, double & dfdxz,
+ double & dfdyx, double & dfdyy, double & dfdyz,
+ double & dfdzx, double & dfdzy, double & dfdzz) const
+{
+
+
+ double xb, yb, curvRadb;
+ bool in;
+
+ getClosestBndPoint(x, y, z, xb, yb, curvRadb, in);
+
+ const double xx = x-xb, yy = y-yb, distSq = xx*xx+yy*yy;
+ const double dist = in ? -sqrt(distSq) : sqrt(distSq);
+ const double oneOverCurvRad = 1./(curvRadb+dist);
+
+ dfdxx = yy*yy*oneOverCurvRad;
+ dfdxy = -xx*yy*oneOverCurvRad;
+ dfdxz = 0.;
+ dfdyx = dfdxy;
+ dfdyy = xx*xx*oneOverCurvRad;
+ dfdyz = 0.;
+ dfdzx = 0.;
+ dfdzy = 0.;
+ dfdzz = 0.;
+
+}
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 609af08581..d40cb9c239 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -621,10 +621,20 @@ 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 x0, double y0, double c, double t);
~gLevelsetNACA00() {}
double operator () (const double x, const double y, const double z) const;
+// std::vector<double> getValueAndGradients(const double x, const double y, const double z) const;
+ void gradient (double x, double y, double z,
+ double & dfdx, double & dfdy, double & dfdz) const;
+ void hessian (double x, double y, double z,
+ double & dfdxx, double & dfdxy, double & dfdxz,
+ double & dfdyx, double & dfdyy, double & dfdyz,
+ double & dfdzx, double & dfdzy, double & dfdzz ) const;
int type() const { return UNKNOWN; }
+private:
+ void getClosestBndPoint(const double x, const double y, const double z,
+ double &xb, double &yb, double &curvRad, bool &in) const;
};
#endif
--
GitLab