diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index a724c9beb94cd85c90619a5008b74550aaa92646..b96ace6065ddfbe7fdbb0c2d235af1dbb5f53030 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -175,29 +175,14 @@ SPoint2 GEdge::reparamOnFace(const GFace *face, double epar,int dir) const
double GEdge::curvature(double par) const
{
- double eps1 = 1.e-5;
- double eps2 = 1.e-5;
-
- Range<double> r = parBounds(0);
- if(r.low() == par) eps2 = 0;
- if(r.high() == par) eps1 = 0;
-
- SVector3 n1 = firstDer(par - eps1);
- SVector3 n2 = firstDer(par + eps2);
-
- GPoint P1 = point(par - eps1);
- GPoint P2 = point(par + eps2);
-
- double D = sqrt ((P1.x() - P2.x()) * (P1.x() - P2.x()) +
- (P1.y() - P2.y()) * (P1.y() - P2.y()) +
- (P1.z() - P2.z()) * (P1.z() - P2.z()));
+ SVector3 d1 = firstDer(par);
+ SVector3 d2 = secondDer(par);
+
+ double one_over_norm = 1. / norm(d1);
- n1.normalize();
- n2.normalize();
- const double one_over_D = 1. / D;
- SVector3 d = one_over_D * (n2 - n1);
- return norm(d);
+ SVector3 cross_prod = crossprod(d1,d2);
+ return ( norm(cross_prod) * one_over_norm * one_over_norm * one_over_norm );
}
double GEdge::length(const double &u0, const double &u1, const int nbQuadPoints)
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 2d9772dd049f5ab376ad579751768c35eef41560..d715c6f93ebfdde34d5032b0d524f93c16248abf 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -14,6 +14,7 @@
#include "VertexArray.h"
#include "GmshMatrix.h"
#include "Numeric.h"
+#include "EigSolve.h"
#if defined(HAVE_GMSH_EMBEDDED)
#include "GmshEmbedded.h"
@@ -441,7 +442,7 @@ void GFace::getMeanPlaneData(double plan[3][3]) const
plan[i][j] = meanPlane.plan[i][j];
}
-double GFace::curvature(const SPoint2 ¶m) const
+double GFace::curvatureDiv(const SPoint2 ¶m) const
{
if (geomType() == Plane) return 0;
@@ -449,7 +450,7 @@ double GFace::curvature(const SPoint2 ¶m) const
// curv = div n = dnx/dx + dny/dy + dnz/dz
// dnx/dx = dnx/du du/dx + dnx/dv dv/dx
- const double eps = 1.e-3;
+ const double eps = 1.e-5;
Pair<SVector3,SVector3> der = firstDer(param);
@@ -462,20 +463,82 @@ double GFace::curvature(const SPoint2 ¶m) const
du.normalize();
dv.normalize();
- SVector3 n1 = normal(SPoint2(param.x() - eps, param.y()));
- SVector3 n2 = normal(SPoint2(param.x() + eps, param.y()));
- SVector3 n3 = normal(SPoint2(param.x(), param.y() - eps));
- SVector3 n4 = normal(SPoint2(param.x(), param.y() + eps));
+ SVector3 n1,n2,n3,n4;
+ if( param.x() - eps < 0.0 ) {
+ n1 = normal(SPoint2(param.x(), param.y()));
+ n2 = normal(SPoint2(param.x() + eps, param.y()));
+ } else {
+ n1 = normal(SPoint2(param.x() - eps, param.y()));
+ n2 = normal(SPoint2(param.x(), param.y()));
+ }
+ if( param.y() - eps < 0.0 ) {
+ n3 = normal(SPoint2(param.x(), param.y() ));
+ n4 = normal(SPoint2(param.x(), param.y() + eps));
+ } else {
+ n3 = normal(SPoint2(param.x(), param.y() - eps));
+ n4 = normal(SPoint2(param.x(), param.y() ));
+ }
+
+ SVector3 dndu = 100000 * (n2 - n1);
+ SVector3 dndv = 100000 * (n4 - n3);
- SVector3 dndu = 500 * (n2 - n1);
- SVector3 dndv = 500 * (n4 - n3);
+ SVector3 dudu = SVector3();
+ SVector3 dvdv = SVector3();
+ SVector3 dudv = SVector3();
+ secondDer(param,&dudu,&dvdv,&dudv);
double ddu = dot(dndu,du);
double ddv = dot(dndv,dv);
- double c = std::max(fabs(ddu),fabs(ddv))/detJ;
- // Msg::Info("c = %g detJ %g", c, detJ);
- return c;
+ return ( fabs(ddu) + fabs(ddv) ) / detJ;
+}
+
+double GFace::curvatureMax(const SPoint2 ¶m) const
+{
+ if (geomType() == Plane) return 0.;
+
+ double eigVal[2], eigVec[8];
+ getMetricEigenVectors(param, eigVal, eigVec);
+
+ return fabs(eigVal[1]);
+}
+
+double GFace::curvatures(const SPoint2 ¶m,
+ SVector3 *dirMax,
+ SVector3 *dirMin,
+ double *curvMax,
+ double *curvMin) const
+{
+ Pair<SVector3,SVector3> D1 = firstDer(param);
+
+ if (geomType() == Plane)
+ {
+ *dirMax = D1.first();
+ *dirMin = D1.second();
+ *curvMax = 0.;
+ *curvMin = 0.;
+ return 0.;
+ }
+
+ if (geomType() == Sphere)
+ {
+ *dirMax = D1.first();
+ *dirMin = D1.second();
+ *curvMax = curvatureDiv(param);
+ *curvMin = *curvMax;
+ return *curvMax;
+ }
+
+ double eigVal[2], eigVec[8];
+ getMetricEigenVectors(param, eigVal, eigVec);
+
+ // curvatures and main directions
+ *curvMax = fabs(eigVal[1]);
+ *curvMin = fabs(eigVal[0]);
+ *dirMax = eigVec[1] * D1.first() + eigVec[3] * D1.second();
+ *dirMin = eigVec[0] * D1.first() + eigVec[2] * D1.second();
+
+ return *curvMax;
}
double GFace::getMetricEigenvalue(const SPoint2 &)
@@ -484,6 +547,73 @@ double GFace::getMetricEigenvalue(const SPoint2 &)
return 0.;
}
+// eigen values are absolute values and sorted from min to max of absolute values
+// eigen vectors are the corresponding COLUMNS of eigVec
+void GFace::getMetricEigenVectors(const SPoint2 ¶m,
+ double eigVal[2],
+ double eigVec[4]) const
+{
+ // first derivatives
+ Pair<SVector3,SVector3> D1 = firstDer(param);
+ SVector3 du = D1.first();
+ SVector3 dv = D1.second();
+ SVector3 nor = crossprod(du, dv); nor.normalize();
+
+ // second derivatives
+ SVector3 dudu = SVector3();
+ SVector3 dvdv = SVector3();
+ SVector3 dudv = SVector3();
+ secondDer(param,&dudu,&dvdv,&dudv);
+
+ // first form
+ double form1[2][2];
+ form1[0][0] = normSq(du);
+ form1[1][1] = normSq(dv);
+ form1[0][1] = form1[1][0] = dot(du,dv);
+
+ // second form
+ double form2[2][2];
+ form2[0][0] = dot(nor,dudu);
+ form2[1][1] = dot(nor,dvdv);
+ form2[0][1] = form2[1][0] = dot(nor,dudv);
+
+ // inverse of first form
+ double inv_form1[2][2];
+ double inv_det_form1 = 1. / ( form1[0][0] * form1[1][1] - form1[1][0] * form1[0][1] );
+ inv_form1[0][0] = inv_det_form1 * form1[1][1];
+ inv_form1[1][1] = inv_det_form1 * form1[0][0];
+ inv_form1[1][0] = inv_form1[0][1] = -1 * inv_det_form1 * form1[0][1];
+
+ // N = (inverse of form1) X (form2)
+ double N[4]; // { N00 N01 N10 N11 }
+ N[0] = inv_form1[0][0] * form2[0][0] + inv_form1[0][1] * form2[1][0];
+ N[1] = inv_form1[0][0] * form2[0][1] + inv_form1[0][1] * form2[1][1];
+ N[2] = inv_form1[1][0] * form2[0][0] + inv_form1[1][1] * form2[1][0];
+ N[3] = inv_form1[1][0] * form2[0][1] + inv_form1[1][1] * form2[1][1];
+
+ // eigen values and vectors of N
+ int work1[2];
+ double work2[2];
+ double eigValI[2];
+ if ( EigSolve(2,2,N,eigVal,eigValI,eigVec,work1,work2) != 1 ) {
+ Msg::Warning("Problem in eigen vectors computation");
+ printf(" N: %f %f %f %f\n",N[0],N[1],N[2],N[3]);
+ printf(" * Eigen values:\n %f + i * %f\n %f + i * %f\n",
+ eigVal[0],eigValI[0],eigVal[1],eigValI[1]);
+ printf(" * Eigen vectors (trust it only if eigen values are real):\n");
+ printf(" ( %f, %f ),\n ( %f, %f ).\n",
+ eigVec[0],eigVec[2],eigVec[1],eigVec[3]);
+ throw;
+ }
+ if ( fabs(eigValI[0]) > 1.e-12 || fabs(eigValI[1]) > 1.e-12 ) {
+ Msg::Error("Found imaginary eigenvalues");
+ }
+
+ eigVal[0] = fabs(eigVal[0]);
+ eigVal[1] = fabs(eigVal[1]);
+ EigSort(2, eigVal, eigValI, eigVec);
+}
+
void GFace::XYZtoUV(const double X, const double Y, const double Z,
double &U, double &V, const double relax,
const bool onSurface) const
diff --git a/Geo/GFace.h b/Geo/GFace.h
index b1033498a254a1ab97b50b1cfb62f94c2f008624..814bb3db4589678a6b8157577e7fe6ad7307e3a8 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -128,6 +128,11 @@ class GFace : public GEntity
// stereographic mappings of the sphere that is used in 2D mesh
// generation !
virtual double getMetricEigenvalue(const SPoint2 &);
+
+ // eigen values are absolute values and sorted from min to max of absolute values
+ // eigen vectors are the COLUMNS of eigVec
+ virtual void getMetricEigenVectors(const SPoint2 ¶m,
+ double eigVal[2], double eigVec[4]) const;
// return the parmater location on the face given a point in space
// that is on the face
@@ -145,8 +150,23 @@ class GFace : public GEntity
// return the first derivate of the face at the parameter location
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const = 0;
- // return the curvature i.e. the divergence of the normal
- virtual double curvature(const SPoint2 ¶m) const;
+ // compute the second derivates of the face at the parameter location
+ // (default implementation by central differences)
+ // the derivates have to be allocated before calling this function
+ virtual void secondDer(const SPoint2 ¶m,
+ SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const = 0;
+
+ // return the curvature computed as the divergence of the normal
+ virtual double curvatureDiv(const SPoint2 ¶m) const;
+
+ // return the maximum curvature at a point
+ virtual double curvatureMax(const SPoint2 ¶m) const;
+
+ // compute the min and max curvatures and the corresponding directions
+ // return the max curvature
+ // outputs have to be allocated before calling this function
+ virtual double curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
+ double *curvMax, double *curvMin) const;
// return a type-specific additional information string
virtual std::string getAdditionalInfoString();
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index a98014f468ff7d4cc37e98911e8a792f359e38f4..ac25984062206961790560a660b8a7d7384f4a62 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -559,6 +559,12 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
return Pair<SVector3, SVector3>(dXdu,dXdv);
}
+void GFaceCompound::secondDer(const SPoint2 ¶m,
+ SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
+{
+ Msg::Error("Computation of the second derivatives not implemented for compound faces");
+}
+
static void GFaceCompoundBB(void *a, double*mmin, double*mmax)
{
GFaceCompoundTriangle *t = (GFaceCompoundTriangle *)a;
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index d4ce6647aaa219e0a9fdd41e4eab80d68c1b6a50..e1587a3de0451acdb084d63e8ba58450af837b40 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -70,7 +70,8 @@ public:
virtual ~GFaceCompound();
Range<double> parBounds(int i) const { return Range<double>(0, 1); }
virtual GPoint point(double par1, double par2) const;
- virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
virtual GEntity::GeomType geomType() const { return CompoundSurface; }
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return 0; }
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 3b7f66ad7cc9c30dc67d8734367f2dc32b6c3da1..a43fdf65f65c044e60efc0d479d2430d7668be04 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -601,7 +601,7 @@ static Vertex InterpolateExtrudedSurface(Surface *s, double u, double v)
Vertex InterpolateSurface(Surface *s, double u, double v, int derivee, int u_v)
{
- if(derivee) {
+ if( derivee == 1 ) {
double eps = 1.e-6;
Vertex D[4];
if(u_v == 1) {
@@ -629,6 +629,44 @@ Vertex InterpolateSurface(Surface *s, double u, double v, int derivee, int u_v)
(D[1].Pos.Z - D[0].Pos.Z) / eps);
}
+ else if ( derivee == 2 ) {
+ double eps = 1.e-6;
+ Vertex D[2];
+ if(u_v == 1) { // dudu
+ if(u - eps < 0.0) {
+ D[0] = InterpolateSurface(s, u, v, 1, 1);
+ D[1] = InterpolateSurface(s, u + eps, v, 1, 1);
+ }
+ else {
+ D[0] = InterpolateSurface(s, u - eps, v, 1, 1);
+ D[1] = InterpolateSurface(s, u, v, 1, 1);
+ }
+ }
+ else if(u_v == 2) { // dvdv
+ if(v - eps < 0.0) {
+ D[0] = InterpolateSurface(s, u, v, 1, 2);
+ D[1] = InterpolateSurface(s, u, v + eps, 1, 2);
+ }
+ else {
+ D[0] = InterpolateSurface(s, u, v - eps, 1, 2);
+ D[1] = InterpolateSurface(s, u, v, 1, 2);
+ }
+ }
+ else { // dudv
+ if(v - eps < 0.0) {
+ D[0] = InterpolateSurface(s, u, v, 1, 1);
+ D[1] = InterpolateSurface(s, u, v + eps, 1, 1);
+ }
+ else {
+ D[0] = InterpolateSurface(s, u, v - eps, 1, 1);
+ D[1] = InterpolateSurface(s, u, v, 1, 1);
+ }
+ }
+ return Vertex((D[1].Pos.X - D[0].Pos.X) / eps,
+ (D[1].Pos.Y - D[0].Pos.Y) / eps,
+ (D[1].Pos.Z - D[0].Pos.Z) / eps);
+ }
+
if(s->geometry){
SPoint3 p = s->geometry->point(u, v);
return Vertex(p.x(), p.y(), p.z());
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 9fda93d3a9b5c735eaf795e213c23a0a2b84b5ec..bfb3493b3e6037332cb46e96bde1a6f4011f247e 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -122,6 +122,18 @@ Pair<SVector3,SVector3> OCCFace::firstDer(const SPoint2 ¶m) const
SVector3(dv.X(), dv.Y(), dv.Z()));
}
+void OCCFace::secondDer(const SPoint2 ¶m,
+ SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
+{
+ gp_Pnt pnt;
+ gp_Vec du, dv, duu, dvv, duv;
+ occface->D2(param.x(), param.y(), pnt, du, dv, duu, dvv, duv);
+
+ *dudu = SVector3(duu.X(), duu.Y(), duu.Z());
+ *dvdv = SVector3(dvv.X(), dvv.Y(), dvv.Z());
+ *dudv = SVector3(duv.X(), duv.Y(), duv.Z());
+}
+
GPoint OCCFace::point(double par1, double par2) const
{
double pp[2] = {par1, par2};
@@ -187,7 +199,7 @@ GEntity::GeomType OCCFace::geomType() const
return Unknown;
}
-double OCCFace::curvature(const SPoint2 ¶m) const
+double OCCFace::curvatureMax(const SPoint2 ¶m) const
{
const double eps = 1.e-12;
BRepAdaptor_Surface sf(s, Standard_True);
@@ -200,6 +212,37 @@ double OCCFace::curvature(const SPoint2 ¶m) const
return std::max(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature()));
}
+double OCCFace::curvatures(const SPoint2 ¶m,
+ SVector3 *dirMax,
+ SVector3 *dirMin,
+ double *curvMax,
+ double *curvMin) const
+{
+ const double eps = 1.e-12;
+ BRepAdaptor_Surface sf(s, Standard_True);
+ BRepLProp_SLProps prop(sf, 2, eps);
+ prop.SetParameters (param.x(),param.y());
+
+ if (!prop.IsCurvatureDefined()){
+ return -1.;
+ }
+
+ *curvMax = std::max(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature()));
+ *curvMin = std::min(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature()));
+
+ gp_Dir dMax = gp_Dir();
+ gp_Dir dMin = gp_Dir();
+ prop.CurvatureDirections(dMax,dMin);
+ (*dirMax)[0] = dMax.X();
+ (*dirMax)[1] = dMax.Y();
+ (*dirMax)[2] = dMax.Z();
+ (*dirMin)[0] = dMin.X();
+ (*dirMin)[1] = dMin.Y();
+ (*dirMin)[2] = dMin.Z();
+
+ return *curvMax;
+}
+
bool OCCFace::containsPoint(const SPoint3 &pt) const
{
if(geomType() == Plane){
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index e2c38ceea6ad07e07b673d72c8505cb392de16b7..8a7c41fa077c2626441ba639520e7a812f580c28 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -30,12 +30,15 @@ class OCCFace : public GFace {
virtual GPoint closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const;
virtual bool containsPoint(const SPoint3 &pt) const;
virtual SVector3 normal(const SPoint2 ¶m) const;
- virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ virtual void secondDer(const SPoint2 ¶m, SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
virtual GEntity::GeomType geomType() const;
ModelType getNativeType() const { return OpenCascadeModel; }
void * getNativePtr() const { return (void*)&s; }
virtual SPoint2 parFromPoint(const SPoint3 &) const;
- virtual double curvature(const SPoint2 ¶m) const;
+ virtual double curvatureMax(const SPoint2 ¶m) const;
+ virtual double curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
+ double *curvMax, double *curvMin) const;
surface_params getSurfaceParams() const;
};
diff --git a/Geo/OCCVertex.cpp b/Geo/OCCVertex.cpp
index 16f52706f3158962a17bb4386fb67e9b263667fd..2c75f8cbc52d62412171512aecb5c44aa44d14cf 100644
--- a/Geo/OCCVertex.cpp
+++ b/Geo/OCCVertex.cpp
@@ -45,7 +45,7 @@ double max_surf_curvature(const GVertex *gv, double x, double y, double z,
double curv = 1.e-22;
while(it != faces.end()){
SPoint2 par = gv->reparamOnFace((*it), 1);
- double cc = (*it)->curvature(par);
+ double cc = (*it)->curvatureDiv(par);
if(cc > 0) curv = std::max(curv, cc);
++it;
}
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index 42d5a64f202a9c975b12059e8dfb95737404c8fe..f27143814fb73e7f9eff4c95c43d551a6d731d41 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -89,12 +89,17 @@ class SMetric3 {
setMat(m3);
return *this;
}
- void eig (gmshMatrix<double> &V, gmshVector<double> &S) {
+ void eig (gmshMatrix<double> &V, gmshVector<double> &S) const {
gmshMatrix<double> me(3,3),right(3,3);
gmshVector<double> im(3);
getMat(me);
me.eig(V,S,im,right);
}
+ void print() const {
+ printf(" %f\t%f\t%f\n",_val[0],_val[1],_val[3]);
+ printf(" %f\t%f\t%f\n",_val[1],_val[2],_val[4]);
+ printf(" %f\t%f\t%f\n",_val[3],_val[4],_val[5]);
+ }
};
// scalar product with respect to the metric
diff --git a/Geo/SVector3.h b/Geo/SVector3.h
index 9102ba2844503565dc1306d7d16e8ac48202c831..149180993587ce740af28f3a16228a17dfac836c 100644
--- a/Geo/SVector3.h
+++ b/Geo/SVector3.h
@@ -7,6 +7,7 @@
#define _SVECTOR3_H_
#include "SPoint3.h"
+#include <string>
// concrete class for vector of size 3
class SVector3 {
@@ -21,10 +22,12 @@ class SVector3 {
SVector3(double x, double y, double z) : P(x, y, z) {}
SVector3(double v) : P(v, v, v) {}
SVector3(const double *array) : P(array) {}
+ SVector3(const SVector3& v) : P(v.P) {}
inline double x(void) const { return P.x(); }
inline double y(void) const { return P.y(); }
inline double z(void) const { return P.z(); }
inline double norm() { return sqrt(P[0] * P[0] + P[1] * P[1] + P[2] * P[2]); }
+ inline double normSq() { return (P[0] * P[0] + P[1] * P[1] + P[2] * P[2]); }
double normalize()
{
double n = norm(); if(n){ P[0] /= n; P[1] /= n; P[2] /= n; }
@@ -62,6 +65,8 @@ class SVector3 {
return *this;
}
operator double *() { return P; }
+ void print(std::string name="") const
+ { printf("Vector \'%s\': %f %f %f\n",name.c_str(),P[0],P[1],P[2]); }
};
inline double dot(const SVector3 &a, const SVector3 &b)
@@ -70,6 +75,9 @@ inline double dot(const SVector3 &a, const SVector3 &b)
inline double norm(const SVector3 &v)
{ return sqrt(dot(v, v)); }
+inline double normSq(const SVector3 &v)
+{ return dot(v, v); }
+
inline SVector3 crossprod(const SVector3 &a, const SVector3 &b)
{ return SVector3(a.y() * b.z() - b.y() * a.z(),
-(a.x() * b.z() - b.x() * a.z()),
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index d6927f745a554858c9fa8b81e6ad261fc6be1f0b..ead70a6610540706030a0e53aecb9ff79d69d9b6 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -100,3 +100,9 @@ Pair<SVector3, SVector3> discreteFace::firstDer(const SPoint2 ¶m) const
Msg::Error("Cannot evaluate derivative on discrete face");
return Pair<SVector3, SVector3>();
}
+
+void discreteFace::secondDer(const SPoint2 ¶m,
+ SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
+{
+ Msg::Error("Cannot evaluate derivative on discrete face");
+}
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index bca7f3cac8d39b18e8de1a089096cd945ad6f20d..2e1b4309eafc7da4e46ce3723c8fb86482f144a4 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -19,6 +19,8 @@ class discreteFace : public GFace {
virtual SVector3 normal(const SPoint2 ¶m) const;
virtual GEntity::GeomType geomType() const { return DiscreteSurface; }
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ virtual void secondDer(const SPoint2 ¶m,
+ SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
void setBoundEdges( std::vector<discreteEdge*> edges );
};
diff --git a/Geo/fourierFace.cpp b/Geo/fourierFace.cpp
index a564ab9a6aafca0395b78ab044ae13de2a5a031c..4203b1acfff3ae25441845a2a17a6cf1cff4657c 100644
--- a/Geo/fourierFace.cpp
+++ b/Geo/fourierFace.cpp
@@ -71,4 +71,10 @@ Pair<SVector3, SVector3> fourierFace::firstDer(const SPoint2 ¶m) const
return Pair<SVector3, SVector3>();
}
+void fourierFace::secondDer(const SPoint2 ¶m,
+ SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
+{
+ Msg::Error("Computation of the second derivatives not implemented for fourer face");
+}
+
#endif
diff --git a/Geo/fourierFace.h b/Geo/fourierFace.h
index 09412d30dcaf301faaaad7f78b249b6dbe10e1d9..2c5cfb017167f20fbbde5b72751b73c52f102498 100644
--- a/Geo/fourierFace.h
+++ b/Geo/fourierFace.h
@@ -29,6 +29,7 @@ class fourierFace : public GFace {
virtual SVector3 normal(const SPoint2 ¶m) const;
virtual GEntity::GeomType geomType() const;
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
ModelType getNativeType() const { return FourierModel; }
void * getNativePtr() const { return face; }
};
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 3d46a89b6bada39a85124941aa2477fbea8a1665..923e1e87334a50acf7d1e0b1ab4d1da21d3a0239 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -165,6 +165,24 @@ Pair<SVector3,SVector3> gmshFace::firstDer(const SPoint2 ¶m) const
}
}
+void gmshFace::secondDer(const SPoint2 ¶m,
+ SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
+{
+ if(s->Typ == MSH_SURF_PLAN && !s->geometry){
+ *dudu = SVector3(0.,0.,0.);
+ *dvdv = SVector3(0.,0.,0.);
+ *dudv = SVector3(0.,0.,0.);
+ }
+ else{
+ Vertex vuu = InterpolateSurface(s, param[0], param[1], 2, 1);
+ Vertex vvv = InterpolateSurface(s, param[0], param[1], 2, 2);
+ Vertex vuv = InterpolateSurface(s, param[0], param[1], 2, 3);
+ *dudu = SVector3(vuu.Pos.X,vuu.Pos.Y,vuu.Pos.Z);
+ *dvdv = SVector3(vvv.Pos.X,vvv.Pos.Y,vvv.Pos.Z);
+ *dudv = SVector3(vuv.Pos.X,vuv.Pos.Y,vuv.Pos.Z);
+ }
+}
+
GPoint gmshFace::point(double par1, double par2) const
{
double pp[2] = {par1, par2};
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index cd2fb927180c25160956219737993377ba2a682a..c5c1aba83fd7c987cee9478bd9ccade6a4bcf8d9 100644
--- a/Geo/gmshFace.h
+++ b/Geo/gmshFace.h
@@ -26,6 +26,7 @@ class gmshFace : public GFace {
virtual double getMetricEigenvalue(const SPoint2 &);
virtual SVector3 normal(const SPoint2 ¶m) const;
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ virtual void secondDer(const SPoint2 ¶m, SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
virtual GEntity::GeomType geomType() const;
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return s; }
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 8e8b718777b6aa92d9965a5e67f64d4992f2cb78..5183002a463138f0dd4b88073c5f70e058b5cbee 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -55,19 +55,19 @@ void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace
pts[1]->X, pts[1]->Y, pts[1]->Z,
pts[2]->X, pts[2]->Y, pts[2]->Z,
pts[3]->X, pts[3]->Y, pts[3]->Z,
- gf->curvature(SPoint2(pts[0]->u, pts[0]->v)),
- gf->curvature(SPoint2(pts[1]->u, pts[1]->v)),
- gf->curvature(SPoint2(pts[2]->u, pts[2]->v)),
- gf->curvature(SPoint2(pts[3]->u, pts[3]->v)));
+ gf->curvatureDiv(SPoint2(pts[0]->u, pts[0]->v)),
+ gf->curvatureDiv(SPoint2(pts[1]->u, pts[1]->v)),
+ gf->curvatureDiv(SPoint2(pts[2]->u, pts[2]->v)),
+ gf->curvatureDiv(SPoint2(pts[3]->u, pts[3]->v)));
}
else{
fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
pts[0]->X, pts[0]->Y, pts[0]->Z,
pts[1]->X, pts[1]->Y, pts[1]->Z,
pts[2]->X, pts[2]->Y, pts[2]->Z,
- gf->curvature(SPoint2(pts[0]->u, pts[0]->v)),
- gf->curvature(SPoint2(pts[1]->u, pts[1]->v)),
- gf->curvature(SPoint2(pts[2]->u, pts[2]->v)));
+ gf->curvatureDiv(SPoint2(pts[0]->u, pts[0]->v)),
+ gf->curvatureDiv(SPoint2(pts[1]->u, pts[1]->v)),
+ gf->curvatureDiv(SPoint2(pts[2]->u, pts[2]->v)));
}
}
}
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 5cb04d200006b2a8f22ec6f9d502aec77b7d41fe..6e68114106d52f3247775c1f7a18b9f29a1e2bd4 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -41,7 +41,7 @@ static double max_surf_curvature(const GEdge *ge, double u)
std::list<GFace *>::iterator it = faces.begin();
while(it != faces.end()){
SPoint2 par = ge->reparamOnFace((*it), u, 1);
- double cc = (*it)->curvature(par);
+ double cc = (*it)->curvatureDiv(par);
val = std::max(cc, val);
++it;
}
@@ -90,7 +90,7 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
case 2:
{
GFace *gf = (GFace *)ge;
- Crv = gf->curvature(SPoint2(U, V));
+ Crv = gf->curvatureDiv(SPoint2(U, V));
}
break;
}