From df5d9b8c123c52516beea6c6a759266dd0a50c2c Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 14 Feb 2014 09:08:44 +0000
Subject: [PATCH] modify newton_fd
---
Geo/Geo.cpp | 10 ++++---
Geo/intersectCurveSurface.cpp | 50 +++++++++++++++++++----------------
Geo/intersectCurveSurface.h | 29 +++++++++++---------
Numeric/Numeric.cpp | 8 +++---
Numeric/Numeric.h | 2 +-
5 files changed, 54 insertions(+), 45 deletions(-)
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 41e285a3f9..1a766e79bd 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -3517,7 +3517,7 @@ struct PointSurface{
Surface *s;
};
-static void projectPS(fullVector<double> &x, fullVector<double> &res, void *data)
+static bool projectPS(fullVector<double> &x, fullVector<double> &res, void *data)
{
PointSurface *ps = (PointSurface*)data;
Vertex c = InterpolateSurface(ps->s, x(0), x(1), 0, 0);
@@ -3531,6 +3531,7 @@ static void projectPS(fullVector<double> &x, fullVector<double> &res, void *data
(c.Pos.X - ps->p->Pos.X) * dv.Pos.X +
(c.Pos.Y - ps->p->Pos.Y) * dv.Pos.Y +
(c.Pos.Z - ps->p->Pos.Z) * dv.Pos.Z;
+ return true;
}
bool ProjectPointOnSurface(Surface *s, Vertex &p, double uv[2])
@@ -3725,7 +3726,7 @@ bool SplitCurve(int line_id, List_T *vertices_id, List_T *shapes)
}
}
List_Delete(Surfs);
-
+
// replace original curve by the new curves in physical groups
for(int i = 0; i < List_Nbr(GModel::current()->getGEOInternals()->PhysicalGroups); i++){
PhysicalGroup *p = *(PhysicalGroup**)List_Pointer
@@ -3742,7 +3743,7 @@ bool SplitCurve(int line_id, List_T *vertices_id, List_T *shapes)
}
}
}
-
+
DeleteShape(c->Typ, c->Num);
List_Delete(new_list);
List_Delete(rshapes);
@@ -3757,7 +3758,7 @@ struct CurveSurface{
Surface *s;
};
-static void intersectCS(fullVector<double> &uvt, fullVector<double> &res, void *data)
+static bool intersectCS(fullVector<double> &uvt, fullVector<double> &res, void *data)
{
CurveSurface *cs = (CurveSurface*)data;
Vertex vs = InterpolateSurface(cs->s, uvt(0), uvt(1), 0, 0);
@@ -3765,6 +3766,7 @@ static void intersectCS(fullVector<double> &uvt, fullVector<double> &res, void *
res(0) = vs.Pos.X - vc.Pos.X;
res(1) = vs.Pos.Y - vc.Pos.Y;
res(2) = vs.Pos.Z - vc.Pos.Z;
+ return true;
}
diff --git a/Geo/intersectCurveSurface.cpp b/Geo/intersectCurveSurface.cpp
index c125bef4ab..0d1ba82405 100644
--- a/Geo/intersectCurveSurface.cpp
+++ b/Geo/intersectCurveSurface.cpp
@@ -1,20 +1,24 @@
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
#include "intersectCurveSurface.h"
#include "Numeric.h"
-static void _kaboom(fullVector<double> &uvt,
- fullVector<double> &res, void *_data);
+static bool _kaboom(fullVector<double> &uvt,
+ fullVector<double> &res, void *_data);
-struct intersectCurveSurfaceData
+struct intersectCurveSurfaceData
{
const curveFunctor &c;
const surfaceFunctor &s;
const double epsilon;
- intersectCurveSurfaceData (const curveFunctor & _c,
- const surfaceFunctor & _s,
- const double &eps) : c(_c),s(_s),epsilon(eps)
- { }
-
- bool apply (double newPoint[3]){
+ intersectCurveSurfaceData(const curveFunctor & _c,
+ const surfaceFunctor & _s,
+ const double &eps) : c(_c),s(_s),epsilon(eps){}
+ bool apply(double newPoint[3])
+ {
try {
fullVector<double> uvt(3);
uvt(0) = newPoint[0];
@@ -24,8 +28,8 @@ struct intersectCurveSurfaceData
_kaboom(uvt,res,this);
// printf("start with %12.5E\n",res.norm());
if (res.norm() < epsilon)return true;
-
-
+
+
if(newton_fd(_kaboom, uvt, this)){
// printf("--- CONVERGED -----------\n");
newPoint[0] = uvt(0);
@@ -33,31 +37,31 @@ struct intersectCurveSurfaceData
newPoint[2] = uvt(2);
// printf("newton done\n");
return true;
- }
+ }
}
catch (...){
- // printf("intersect curve surface failed !\n");
+ // printf("intersect curve surface failed !\n");
}
- // printf("newton failsed\n");
+ // printf("newton failed\n");
return false;
}
};
-void _kaboom(fullVector<double> &uvt,
- fullVector<double> &res, void *_data){
+static bool _kaboom(fullVector<double> &uvt,
+ fullVector<double> &res, void *_data)
+{
intersectCurveSurfaceData *data = (intersectCurveSurfaceData*)_data;
- // printf("coucuo1 %g %g\n",uvt(0),uvt(1));
- SPoint3 s = data->s(uvt(0),uvt(1));
- // printf("coucuo1\n");
+ SPoint3 s = data->s(uvt(0), uvt(1));
SPoint3 c = data->c(uvt(2));
- // printf("coucuo1\n");
res(0) = s.x() - c.x();
res(1) = s.y() - c.y();
res(2) = s.z() - c.z();
- // printf("%g %g %g vs %g %g %g \n",s.x(),s.y(),s.z(),c.x(),c.y(),c.z());
+ return true;
}
-int intersectCurveSurface (curveFunctor &c, surfaceFunctor & s, double uvt[3], double epsilon){
- intersectCurveSurfaceData data(c,s,epsilon);
+int intersectCurveSurface(curveFunctor &c, surfaceFunctor & s, double uvt[3],
+ double epsilon)
+{
+ intersectCurveSurfaceData data(c, s, epsilon);
return data.apply(uvt);
}
diff --git a/Geo/intersectCurveSurface.h b/Geo/intersectCurveSurface.h
index c06c10cdf3..2c70b5fb74 100644
--- a/Geo/intersectCurveSurface.h
+++ b/Geo/intersectCurveSurface.h
@@ -8,20 +8,20 @@
// Intersection of a curve and a surface using newton's method
// FD are used to compute derivatives of the parametrizations
+#include <math.h>
#include "SPoint3.h"
#include "SVector3.h"
#include "GFace.h"
#include "GEdge.h"
-#include "math.h"
-class surfaceFunctor
+class surfaceFunctor
{
public :
virtual ~surfaceFunctor(){}
virtual SPoint3 operator () (double u, double v) const = 0;
};
-class curveFunctor
+class curveFunctor
{
public :
virtual ~curveFunctor(){}
@@ -57,10 +57,11 @@ class curveFunctorCircle : public curveFunctor
SVector3 middle;
double d;
public :
- curveFunctorCircle(const SVector3 & _n1, const SVector3 & _n2, const SVector3 & _middle, const double & _d) :
- n1(_n1),n2(_n2),middle(_middle),d(_d) {
- }
- virtual SPoint3 operator () (double t) const {
+ curveFunctorCircle(const SVector3 & _n1, const SVector3 & _n2,
+ const SVector3 & _middle, const double & _d) :
+ n1(_n1), n2(_n2), middle(_middle), d(_d){}
+ virtual SPoint3 operator () (double t) const
+ {
SVector3 dir = (n1*cos(t)+n2*sin(t))*d;
return SPoint3(middle.x() + dir.x(),
middle.y() + dir.y(),
@@ -73,17 +74,21 @@ class surfaceFunctorPlane : public surfaceFunctor
const SPoint3 p;
const SVector3 v1,v2;
public :
- surfaceFunctorPlane(const SPoint3 &_p, const SVector3 &_v1, const SVector3 &_v2) : p(_p),v1(_v1),v2(_v2) {}
- virtual SPoint3 operator () (double u, double v) const {
+ surfaceFunctorPlane(const SPoint3 &_p, const SVector3 &_v1, const SVector3 &_v2)
+ : p(_p),v1(_v1),v2(_v2) {}
+ virtual SPoint3 operator () (double u, double v) const
+ {
return SPoint3 (p.x() + u * v1.x() + v * v2.x(),
- p.y() + u * v1.y() + v * v2.y(),
- p.z() + u * v1.z() + v * v2.z()) ;
+ p.y() + u * v1.y() + v * v2.y(),
+ p.z() + u * v1.z() + v * v2.z()) ;
}
};
+
// intersects the curve and the surface using Newton.
// the initial guess should be a good guess
// returns 1 --> OK
// returns 0 --> NOT CONVERGED
-int intersectCurveSurface (curveFunctor &c, surfaceFunctor & s, double uvt[3], double epsilon);
+int intersectCurveSurface(curveFunctor &c, surfaceFunctor & s, double uvt[3],
+ double epsilon);
#endif
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 5827e151d2..9b26a8a5fb 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -656,7 +656,7 @@ void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
}
}
-bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
+bool newton_fd(bool (*func)(fullVector<double> &, fullVector<double> &, void *),
fullVector<double> &x, void *data, double relax, double tolx)
{
const int MAXIT = 10;
@@ -667,7 +667,7 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
fullVector<double> f(N), feps(N), dx(N);
for (int iter = 0; iter < MAXIT; iter++){
- func(x, f, data);
+ if(!func(x, f, data)) return false;
bool isZero = false;
for (int k=0; k<N; k++) {
@@ -681,10 +681,9 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
double h = EPS * fabs(x(j));
if(h == 0.) h = EPS;
x(j) += h;
- func(x, feps, data);
+ if(!func(x, feps, data)) return false;
for (int i = 0; i < N; i++){
J(i, j) = (feps(i) - f(i)) / h;
- // printf("J(%d,%d) = %12.5E \n",i,j,J(i,j));
}
x(j) -= h;
}
@@ -700,7 +699,6 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
for (int i = 0; i < N; i++)
x(i) -= relax * dx(i);
- // printf("dx = %12.5E\n",dx.norm());
if(dx.norm() < tolx) return true;
}
return false;
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 3e100539c6..a19affad2e 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -106,7 +106,7 @@ void gradSimplex(double *x, double *y, double *z, double *v, double *grad);
double ComputeVonMises(double *val);
double ComputeScalarRep(int numComp, double *val);
void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
-bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
+bool newton_fd(bool (*func)(fullVector<double> &, fullVector<double> &, void *),
fullVector<double> &x, void *data, double relax=1., double tolx=1.e-6);
double minimize_grad_fd(double (*func)(fullVector<double> &, void *),
fullVector<double> &x, void *data);
--
GitLab