From e4a4128462ad62206f49aa91d1ef2f3ad5afcc89 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Thu, 25 Apr 2013 13:01:05 +0000
Subject: [PATCH] Added analytical gradients and Hessian for sphere level-set
---
Geo/gmshLevelset.cpp | 40 ++++++++++++++++++++++++++++++++++++++++
Geo/gmshLevelset.h | 8 +++++++-
2 files changed, 47 insertions(+), 1 deletion(-)
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 97be674aee..d62dd21e44 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -314,6 +314,46 @@ gLevelset::gLevelset(const gLevelset &lv)
tag_ = lv.tag_;
}
+gLevelsetSphere::gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag=1)
+ : gLevelsetPrimitive(tag), xc(x), yc(y), zc(z), r(R)
+{
+_hasDerivatives = true;
+}
+
+void gLevelsetSphere::gradient (double x, double y, double z,
+ double & dfdx, double & dfdy, double & dfdz) const
+{
+
+ const double xx = x-xc, yy = y-yc, zz = z-zc, dist = sqrt(xx*xx+yy*yy+zz*zz);
+
+ dfdx = xx/dist;
+ dfdy = yy/dist;
+ dfdz = zz/dist;
+
+}
+
+void gLevelsetSphere::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
+{
+
+ const double xx = x-xc, yy = y-yc, zz = z-zc;
+ const double distSq = xx*xx+yy*yy, fact = 1./(distSq*sqrt(distSq));
+
+ dfdxx = (zz*zz+yy*yy)*fact;
+ dfdxy = -xx*yy*fact;
+ dfdxz = -xx*zz*fact;
+ dfdyx = dfdxy;
+ dfdyy = (xx*xx+zz*zz)*fact;
+ dfdyz = -yy*zz*fact;
+ dfdzx = dfdxz;
+ dfdzy = dfdyz;
+ dfdzz = (xx*xx+yy*yy)*fact;
+
+}
+
+
gLevelsetPlane::gLevelsetPlane(const std::vector<double> &pt,
const std::vector<double> &norm, int tag)
: gLevelsetPrimitive(tag)
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index d40cb9c239..d6716c50b1 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -128,12 +128,18 @@ class gLevelsetSphere : public gLevelsetPrimitive
protected:
double xc, yc, zc, r;
public:
- gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag=1) : gLevelsetPrimitive(tag), xc(x), yc(y), zc(z), r(R) {}
+ gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag);
virtual double operator () (const double x, const double y, const double z) const{
if(r >= 0.)
return sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) + (zc - z) * (zc - z)) - r;
return (- r - sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) + (zc - z) * (zc - z)));
}
+ 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 SPHERE;}
};
--
GitLab