diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 3c02029ca90a6794912af4574ab7ba85d1199bba..ddf640c0d31d4550e56a16b0fd89ba1aa938fabf 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2877,51 +2877,40 @@ void ReplaceAllDuplicates()
ReplaceDuplicateSurfaces();
}
-
// Projection of a point on a surface
-static Vertex *VERTEX;
-static Curve *CURVE;
-static Surface *SURFACE;
-
-static void projectPS(int N, double x[], double res[])
-{
- //x[1] = u x[2] = v
- Vertex du, dv, c;
- c = InterpolateSurface(SURFACE, x[1], x[2], 0, 0);
- du = InterpolateSurface(SURFACE, x[1], x[2], 1, 1);
- dv = InterpolateSurface(SURFACE, x[1], x[2], 1, 2);
- res[1] =
- (c.Pos.X - VERTEX->Pos.X) * du.Pos.X +
- (c.Pos.Y - VERTEX->Pos.Y) * du.Pos.Y +
- (c.Pos.Z - VERTEX->Pos.Z) * du.Pos.Z;
- res[2] =
- (c.Pos.X - VERTEX->Pos.X) * dv.Pos.X +
- (c.Pos.Y - VERTEX->Pos.Y) * dv.Pos.Y +
- (c.Pos.Z - VERTEX->Pos.Z) * dv.Pos.Z;
-}
-
-bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2])
-{
- double x[3] = { 0.5, u[0], u[1] };
- int check;
- SURFACE = s;
- VERTEX = &p;
- newt(x, 2, &check, projectPS);
-
- Vertex vv = InterpolateSurface(s, x[1], x[2], 0, 0);
- double res[3];
- projectPS(2, x, res);
- double resid = sqrt(res[1] * res[1] + res[2] * res[2]);
-
- p.Pos.X = vv.Pos.X;
- p.Pos.Y = vv.Pos.Y;
- p.Pos.Z = vv.Pos.Z;
- u[0] = x[1];
- u[1] = x[2];
- if(resid > 1.e-6)
- return false;
- return true;
+struct PointSurface{
+ Vertex *p;
+ Surface *s;
+};
+
+static void projectPS(Double_Vector &x, Double_Vector &res, void *data)
+{
+ PointSurface *ps = (PointSurface*)data;
+ Vertex c = InterpolateSurface(ps->s, x(0), x(1), 0, 0);
+ Vertex du = InterpolateSurface(ps->s, x(0), x(1), 1, 1);
+ Vertex dv = InterpolateSurface(ps->s, x(0), x(1), 1, 2);
+ res(0) =
+ (c.Pos.X - ps->p->Pos.X) * du.Pos.X +
+ (c.Pos.Y - ps->p->Pos.Y) * du.Pos.Y +
+ (c.Pos.Z - ps->p->Pos.Z) * du.Pos.Z;
+ res(1) =
+ (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;
+}
+
+bool ProjectPointOnSurface(Surface *s, Vertex &p, double uv[2])
+{
+ Double_Vector x(2);
+ x(0) = uv[0];
+ x(1) = uv[1];
+ PointSurface ps = {&p, s};
+ if(newton_fd(projectPS, x, &ps, 1.)){
+ p = InterpolateSurface(s, x(0), x(1), 0, 0);
+ return true;
+ }
+ return false;
}
// Split line
@@ -3039,25 +3028,19 @@ bool SplitCurve(int line_id, List_T *vertices_id, List_T *shapes)
// Intersect a curve with a surface
-static void intersectCS(int N, double x[], double res[])
-{
- // (x[1], x[2]) = surface params, x[3] = curve param
- Vertex s = InterpolateSurface(SURFACE, x[1], x[2], 0, 0);
- Vertex c = InterpolateCurve(CURVE, x[3], 0);
- res[1] = s.Pos.X - c.Pos.X;
- res[2] = s.Pos.Y - c.Pos.Y;
- res[3] = s.Pos.Z - c.Pos.Z;
-}
+struct CurveSurface{
+ Curve *c;
+ Surface *s;
+};
-static bool IntersectCurveSurface(Curve *c, Surface *s, double x[4])
+static void intersectCS(Double_Vector &uvt, Double_Vector &res, void *data)
{
- int check;
- SURFACE = s;
- CURVE = c;
- newt(x, 3, &check, intersectCS);
-
- if(check) return false;
- return true;
+ CurveSurface *cs = (CurveSurface*)data;
+ Vertex vs = InterpolateSurface(cs->s, uvt(0), uvt(1), 0, 0);
+ Vertex vc = InterpolateCurve(cs->c, uvt(2), 0);
+ res(0) = vs.Pos.X - vc.Pos.X;
+ res(1) = vs.Pos.Y - vc.Pos.Y;
+ res(2) = vs.Pos.Z - vc.Pos.Z;
}
bool IntersectCurvesWithSurface(List_T *curve_ids, int surface_id, List_T *shapes)
@@ -3072,15 +3055,19 @@ bool IntersectCurvesWithSurface(List_T *curve_ids, int surface_id, List_T *shape
List_Read(curve_ids, i, &curve_id);
Curve *c = FindCurve((int)curve_id);
if(c){
- double x[4] = {0., 0.5, 0.5, 0.5};
- if(IntersectCurveSurface(c, s, x)){
- Vertex p = InterpolateCurve(c, x[3], 0);
+ CurveSurface cs = {c, s};
+ Double_Vector uvt(3);
+ uvt(0) = 0.5;
+ uvt(1) = 0.5;
+ uvt(2) = 0.5;
+ if(newton_fd(intersectCS, uvt, &cs, 1.)){
+ Vertex p = InterpolateCurve(c, uvt(2), 0);
Vertex *v = Create_Vertex(NEWPOINT(), p.Pos.X, p.Pos.Y, p.Pos.Z, p.lc, p.u);
Tree_Insert(GModel::current()->getGEOInternals()->Points, &v);
- Shape s;
- s.Type = MSH_POINT;
- s.Num = v->Num;
- List_Add(shapes, &s);
+ Shape sh;
+ sh.Type = MSH_POINT;
+ sh.Num = v->Num;
+ List_Add(shapes, &sh);
}
}
else{
diff --git a/Geo/Geo.h b/Geo/Geo.h
index d210c9b46f6910b3ab69514394f709cab23378a3..c9456c96be64df98c3ae4318ecbff735013c6cb2 100644
--- a/Geo/Geo.h
+++ b/Geo/Geo.h
@@ -275,8 +275,7 @@ void ProtudeXYZ(double &x, double &y, double &z, ExtrudeParams *e);
void ReplaceAllDuplicates();
-bool ProjectPointOnCurve(Curve *c, Vertex *v, Vertex *RES, Vertex *DER);
-bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2]);
+bool ProjectPointOnSurface(Surface *s, Vertex &p, double uv[2]);
bool IntersectCurvesWithSurface(List_T *curve_ids, int surface_id, List_T *shapes);
bool SplitCurve(int line_id, List_T *vertices_id, List_T *shapes);
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 9b293baf9872c148c87ad7dec6c3077d41ddea4d..b130fabf67cadc2c0e72660e27a8f5bd9d8bdcf1 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -36,90 +36,6 @@ int Init_Numeric()
return 1;
}
-#define MAX_DIM_NEWT 10
-#define MAXITER 100
-#define PREC 1.e-8
-
-static int nrdim;
-static double nru[MAX_DIM_NEWT], nrv[MAX_DIM_NEWT];
-static void (*nrfunc) (int n, double x[], double y[]);
-struct gsl_dummy{ int i; };
-
-static void convert_vector_from_gsl(const gsl_vector * gv, double *v)
-{
- int i, m;
- m = gv->size;
- for(i = 0; i < m; i++) {
- v[i + 1] = gsl_vector_get(gv, i);
- }
-}
-
-static void convert_vector_to_gsl(double *v, int n, gsl_vector * gv)
-{
- int i;
- for(i = 0; i < n; i++) {
- gsl_vector_set(gv, i, v[i + 1]);
- }
-}
-
-static int gslfunc(const gsl_vector * xx, void *params, gsl_vector * f)
-{
- convert_vector_from_gsl(xx, nru);
- (*nrfunc) (nrdim, nru, nrv);
- // Msg::Info("f(%lf,%lf) = %lf %lf\n",nru[1],nru[2],nrv[1],nrv[2]);
- convert_vector_to_gsl(nrv, nrdim, f);
- return GSL_SUCCESS;
-}
-
-// Warning: for compatibility with the old newt from NR, x[] and the
-// arguments of func() are indexed from 1 to n...
-void newt(double x[], int n, int *check,
- void (*func) (int, double[], double[]))
-{
- const gsl_multiroot_fsolver_type *T;
- gsl_multiroot_fsolver *s;
- int status;
- size_t iter = 0;
- struct gsl_dummy p = { 1 };
- gsl_multiroot_function f = { &gslfunc, n, &p };
- gsl_vector *xx = gsl_vector_alloc(n);
-
- if(n > MAX_DIM_NEWT - 1)
- Msg::Fatal("Maximum Newton dimension exceeded\n");
- nrdim = n;
-
- nrfunc = func;
- convert_vector_to_gsl(x, n, xx);
-
- T = gsl_multiroot_fsolver_hybrid;
- s = gsl_multiroot_fsolver_alloc(T, n);
- gsl_multiroot_fsolver_set(s, &f, xx);
-
- do {
- iter++;
- status = gsl_multiroot_fsolver_iterate(s);
- // Msg::Info("status %d %d %d %lf %lf\n",
- // status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1));
- if(status)
- break; // solver problem
- status = gsl_multiroot_test_residual(s->f, n * PREC);
- }
- while(status == GSL_CONTINUE && iter < MAXITER);
-
- if(status == GSL_CONTINUE) {
- *check = 1; // problem !!!
- }
- else {
- // Msg::Info("status %d %d %d %lf %lf\n",
- // status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1));
- convert_vector_from_gsl(s->x, x);
- *check = 0; // converged
- }
-
- gsl_multiroot_fsolver_free(s);
- gsl_vector_free(xx);
-}
-
static double (*f_statN) (double*, void *data);
static void (*df_statN) (double*, double*, double &, void *data);
@@ -196,18 +112,6 @@ int Init_Numeric()
return 1;
}
-double brent(double ax, double bx, double cx,
- double (*f) (double), double tol, double *xmin)
-{
- Msg::Error("Gmsh must be compiled with GSL support for brent");
-}
-
-void newt(double x[], int n, int *check,
- void (*func) (int, double[], double[]))
-{
- Msg::Error("Gmsh must be compiled with GSL support for newt");
-}
-
void minimize_N(int N, double (*f) (double*, void *data),
void (*df) (double*, double*, double &, void *data) ,
void *data,int niter,
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index aebc115d68282ea09100e8bc33ff1da85e942c89..0348210f6e65af796ba7ee74852fa48970346bb4 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -9,8 +9,6 @@
#include "NumericEmbedded.h"
int Init_Numeric();
-void newt(double x[], int n, int *check,
- void (*func) (int, double[], double[]));
void minimize_N(int N, double (*f) (double*, void *data),
void (*df) (double*, double*, double &, void *data),
void *data,int niter,
diff --git a/Numeric/NumericEmbedded.cpp b/Numeric/NumericEmbedded.cpp
index c625836a55866847fe6065109536eeb1f6a9c197..219d1167a913e4818084713e32774e560023aaa9 100644
--- a/Numeric/NumericEmbedded.cpp
+++ b/Numeric/NumericEmbedded.cpp
@@ -8,7 +8,6 @@
// Msg)
#include "NumericEmbedded.h"
-#include "GmshMatrix.h"
#include "GmshMessage.h"
#define SQU(a) ((a)*(a))
@@ -435,7 +434,7 @@ void FindCubicRoots(const double coef[4], double real[3], double imag[3])
real[2] = -term1 + r13*cos((dum1 + 4.0*M_PI)/3.0);
}
-void eigsort(double d[3])
+void eigsort(double d[3])
{
int k, j, i;
double p;
@@ -487,3 +486,42 @@ void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
}
}
}
+
+bool newton_fd(void (*func)(Double_Vector &, Double_Vector &, void *),
+ Double_Vector &x, void *data, double relax)
+{
+ const double PRECISION = 1.e-6;
+ const int MAXIT = 10;
+ const double EPS = 1.e-4;
+ const int N = x.size();
+
+ Double_Matrix J(N, N);
+ Double_Vector r(N), rp(N), dx(N);
+
+ int iter = 1;
+ while (iter < MAXIT){
+ iter++;
+ func(x, r, data);
+
+ for (int j = 0; j < N; j++){
+ double h = EPS * fabs(x(j));
+ if(h == 0.) h = EPS;
+ x(j) += h;
+ func(x, rp, data);
+ for (int i = 0; i < N; i++)
+ J(i, j) = (rp(i) - r(i)) / h;
+ x(j) -= h;
+ }
+
+ if (N == 1)
+ dx(0) = r(0) / J(0, 0);
+ else
+ J.lu_solve(r, dx);
+
+ for (int i = 0; i < N; i++)
+ x(i) -= relax * dx(i);
+
+ if(dx.norm() < PRECISION) return true;
+ }
+ return false;
+}
diff --git a/Numeric/NumericEmbedded.h b/Numeric/NumericEmbedded.h
index ceba0e734ba6f32e3a494f9e511613aa0d64b2d3..fb457310ff5a5dbbc529125ae438de23adc90143 100644
--- a/Numeric/NumericEmbedded.h
+++ b/Numeric/NumericEmbedded.h
@@ -7,6 +7,7 @@
#define _NUMERIC_EMBEDDED_H_
#include <math.h>
+#include "GmshMatrix.h"
#define myhypot(a,b) (sqrt((a)*(a)+(b)*(b)))
#define sign(x) (((x)>=0)?1:-1)
@@ -69,5 +70,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)(Double_Vector &, Double_Vector &, void *),
+ Double_Vector &x, void *data, double relax);
+
#endif