diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 4caab71b01b24c1d03c5fc3cc6a3c26e05d1576d..5e93359a035198ad31b7b1e2d4d5512425fc8fd1 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -119,12 +119,6 @@ std::string GEdge::getAdditionalInfoString()
return sstream.str();
}
-GPoint GEdge::closestPoint(const SPoint3 & queryPoint) const
-{
- Msg::Error("Closet point not implemented for this type of edge");
- return GPoint(0, 0, 0);
-}
-
bool GEdge::containsParam(double pt) const
{
Range<double> rg = parBounds(0);
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 834283d28b9dc98bc4d057110e5221fac7a84005..75fe34602a944bf76fa5836915711907e7e16bba 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -56,16 +56,9 @@ class GEdge : public GEntity {
// faces that this entity bounds
virtual std::list<GFace*> faces() const { return l_faces; }
- // get the parameter location for a point in space on the edge
- // (returns std::numeric_limits<double>::max() if failed)
- virtual double parFromPoint(const SPoint3 &) const = 0;
-
// get the point for the given parameter location
virtual GPoint point(double p) const = 0;
- // get the closest point on the edge to the given point
- virtual GPoint closestPoint(const SPoint3 & queryPoint) const;
-
// true if the edge contains the given parameter
virtual bool containsParam(double pt) const;
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index fe4353a52df44a1d2f352bec0db3f55f0fe74f67..3c02029ca90a6794912af4574ab7ba85d1199bba 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2878,18 +2878,11 @@ void ReplaceAllDuplicates()
}
-// Projection of a point on a curve or a surface
+// Projection of a point on a surface
+static Vertex *VERTEX;
static Curve *CURVE;
static Surface *SURFACE;
-static Vertex *VERTEX;
-
-static double min1d(double (*funct) (double), double *xmin)
-{
- // tolerance==0. allows for maximum as code in gsl_brent
- return brent(CURVE->ubeg, 0.5 * (CURVE->ubeg + CURVE->uend), CURVE->uend,
- funct, 0., xmin);
-}
static void projectPS(int N, double x[], double res[])
{
@@ -2908,35 +2901,6 @@ static void projectPS(int N, double x[], double res[])
(c.Pos.Z - VERTEX->Pos.Z) * dv.Pos.Z;
}
-static double projectPC(double u)
-{
- Vertex c = InterpolateCurve(CURVE, u, 0);
- return sqrt(SQU(c.Pos.X - VERTEX->Pos.X) +
- SQU(c.Pos.Y - VERTEX->Pos.Y) +
- SQU(c.Pos.Z - VERTEX->Pos.Z));
-}
-
-bool ProjectPointOnCurve(Curve *c, Vertex *v, Vertex *RES, Vertex *DER)
-{
- double xmin;
- CURVE = c;
- VERTEX = v;
- min1d(projectPC, &xmin);
- *RES = InterpolateCurve(CURVE, xmin, 0);
- *DER = InterpolateCurve(CURVE, xmin, 1);
- if(xmin > c->uend) {
- xmin = c->uend;
- *RES = InterpolateCurve(CURVE, c->uend, 0);
- *DER = InterpolateCurve(CURVE, c->uend, 1);
- }
- else if(xmin < c->ubeg) {
- xmin = c->ubeg;
- *RES = InterpolateCurve(CURVE, c->ubeg, 0);
- *DER = InterpolateCurve(CURVE, c->ubeg, 1);
- }
- return true;
-}
-
bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2])
{
double x[3] = { 0.5, u[0], u[1] };
diff --git a/Geo/MZoneBoundary.cpp b/Geo/MZoneBoundary.cpp
index f8c70e929835dfa398843288a20c061191fd63e2..b553a5f34a6ea9467490ae264f8d3c7ac52cfafe 100644
--- a/Geo/MZoneBoundary.cpp
+++ b/Geo/MZoneBoundary.cpp
@@ -205,7 +205,7 @@ void updateBoVec
* onlyFace - (I) If >= 0, only use this face to determine the
* interior vertex and normal to the mesh plane.
* returns - (O) 0 - success
- * 1 - parFromPoint() failed
+ * 1 - no parameter found on edge for vertex
*
*============================================================================*/
@@ -215,8 +215,8 @@ int edge_normal
&faces, SVector3 &boNormal, const int onlyFace = -1)
{
- const double par = gEdge->parFromPoint(vertex->point());
- if(par == std::numeric_limits<double>::max()) return 1;
+ double par;
+ if(!reparamMeshVertexOnEdge(vertex, gEdge, par)) return 1;
const SVector3 tangent(gEdge->firstDer(par));
// Tangent to the boundary face
@@ -639,8 +639,9 @@ void updateBoVec<3, MFace>
for(std::list<const GFace*>::const_iterator gFIt = useGFace.begin();
gFIt != useGFace.end(); ++gFIt) {
- const SPoint2 par = (*gFIt)->parFromPoint(vertex->point());
- if(par.x() == std::numeric_limits<double>::max()) //**?
+
+ SPoint2 par;
+ if(!reparamMeshVertexOnFace(vertex, (*gFIt), par))
goto getNormalFromElements; // :P After all that!
SVector3 boNormal = (*gFIt)->normal(par);
@@ -673,8 +674,8 @@ void updateBoVec<3, MFace>
{
const GFace *const gFace = static_cast<const GFace*>(ent);
- const SPoint2 par = gFace->parFromPoint(vertex->point());
- if(par.x() == std::numeric_limits<double>::max())
+ SPoint2 par;
+ if(!reparamMeshVertexOnFace(vertex, gFace, par))
goto getNormalFromElements;
SVector3 boNormal = static_cast<const GFace*>(ent)->normal(par);
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 70d23684bc350783cacc90e64374fec6cfbf1580..14b83453ac163d64dd0c22bec1e5ade459ae6fc1 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -146,12 +146,6 @@ SVector3 OCCEdge::firstDer(double par) const
return SVector3(d1.X(), d1.Y(), d1.Z());
}
-double OCCEdge::parFromPoint(const SPoint3 &pt) const
-{
- Msg::Error("parFromPoint not implemented for OCCEdge");
- return std::numeric_limits<double>::max();
-}
-
GEntity::GeomType OCCEdge::geomType() const
{
if(curve.IsNull()){
diff --git a/Geo/OCCEdge.h b/Geo/OCCEdge.h
index 2e17564849526c581554f2d2e823d771b5fa29e5..0b71d20aab2b82cadc1a445405bedd4ea6309329 100644
--- a/Geo/OCCEdge.h
+++ b/Geo/OCCEdge.h
@@ -35,7 +35,6 @@ class OCCEdge : public GEdge {
virtual SPoint2 reparamOnFace(GFace * face, double epar, int dir) const ;
ModelType getNativeType() const { return OpenCascadeModel; }
void * getNativePtr() const { return (void*)&c; }
- virtual double parFromPoint(const SPoint3 &pt) const;
virtual int minimumMeshSegments () const;
virtual int minimumDrawSegments () const;
bool is3D() const { return !curve.IsNull(); }
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 28b5fc15f80c72e3aa9b6b16ddf05fdfb3deecba..9d6dade36dbf298a62ecd251c25bd1663fd8c275 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -31,8 +31,3 @@ SVector3 discreteEdge::firstDer(double par) const
return SVector3();
}
-double discreteEdge::parFromPoint(const SPoint3 &pt) const
-{
- Msg::Error("Cannot compute parametric coordinate on discrete edge");
- return 0.;
-}
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index b690f8163acd65d3a4ce932677e73e88e5150525..77495657fc0add41bec0492f50d09d5e0c5be948 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -16,7 +16,6 @@ class discreteEdge : public GEdge {
virtual GeomType geomType() const { return DiscreteCurve; }
virtual GPoint point(double p) const;
virtual SVector3 firstDer(double par) const;
- virtual double parFromPoint(const SPoint3 &pt) const;
};
#endif
diff --git a/Geo/fourierEdge.cpp b/Geo/fourierEdge.cpp
index 95daf90540152d606d387ae6cf5afd8783621a5f..39d13e302fdb30f29eec8f833a87c1cd69691aa7 100644
--- a/Geo/fourierEdge.cpp
+++ b/Geo/fourierEdge.cpp
@@ -27,13 +27,6 @@ GPoint fourierEdge::point(double p) const
return GPoint(x,y,z);
}
-double fourierEdge::parFromPoint(const SPoint3 &pt) const
-{
- double p;
- edge->Inverse(pt.x(),pt.y(),pt.z(),p);
- return p;
-}
-
SVector3 fourierEdge::firstDer(double par) const
{
double x,y,z;
diff --git a/Geo/fourierEdge.h b/Geo/fourierEdge.h
index a50e82d5de0d4cd00d5fd4307e7a4fee294dc646..6d61b1c037dc6903e5dfdc5178fc9f2357288b2c 100644
--- a/Geo/fourierEdge.h
+++ b/Geo/fourierEdge.h
@@ -26,7 +26,6 @@ class fourierEdge : public GEdge {
virtual GeomType geomType() const { return ParametricCurve; }
virtual GPoint point(double p) const;
virtual SVector3 firstDer(double par) const;
- virtual double parFromPoint(const SPoint3 &pt) const;
virtual int minimumMeshSegments () const;
virtual int minimumDrawSegments () const;
ModelType getNativeType() const { return FourierModel; }
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 73a8f0952dadc8d1c78e668c079782fb519091f6..92a985975755bb0b0d7ec16efff953c28e2a5d92 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -40,36 +40,12 @@ GPoint gmshEdge::point(double par) const
return GPoint(a.Pos.X, a.Pos.Y, a.Pos.Z, this, par);
}
-GPoint gmshEdge::closestPoint(const SPoint3 &qp) const
-{
- Vertex v;
- Vertex a;
- Vertex der;
- v.Pos.X = qp.x();
- v.Pos.Y = qp.y();
- v.Pos.Z = qp.z();
- ProjectPointOnCurve(c, &v, &a, &der);
- return GPoint(a.Pos.X, a.Pos.Y, a.Pos.Z, this, a.u);
-}
-
SVector3 gmshEdge::firstDer(double par) const
{
Vertex a = InterpolateCurve(c, par, 1);
return SVector3(a.Pos.X, a.Pos.Y, a.Pos.Z);
}
-double gmshEdge::parFromPoint(const SPoint3 &pt) const
-{
- Vertex v;
- Vertex a;
- Vertex der;
- v.Pos.X = pt.x();
- v.Pos.Y = pt.y();
- v.Pos.Z = pt.z();
- ProjectPointOnCurve(c, &v, &a, &der);
- return a.u;
-}
-
GEntity::GeomType gmshEdge::geomType() const
{
switch (c->Typ){
diff --git a/Geo/gmshEdge.h b/Geo/gmshEdge.h
index fc81deabed7c66d43ac9907455777cd5a7f12453..684d64eace9458297481253b831178455864dde8 100644
--- a/Geo/gmshEdge.h
+++ b/Geo/gmshEdge.h
@@ -19,11 +19,9 @@ class gmshEdge : public GEdge {
virtual Range<double> parBounds(int i) const;
virtual GeomType geomType() const;
virtual GPoint point(double p) const;
- virtual GPoint closestPoint(const SPoint3 & queryPoint) const;
virtual SVector3 firstDer(double par) const;
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return c; }
- virtual double parFromPoint(const SPoint3 &pt) const;
virtual int minimumMeshSegments() const;
virtual int minimumDrawSegments() const;
virtual void resetMeshAttributes();
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index f72229c98f7ef440b81e9efb85da151c58dec834..851635ceb1bea57dec2a79979d63c873900c5ba6 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -186,7 +186,7 @@ GPoint gmshFace::point(double par1, double par2) const
GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2]) const
{
- if (s->Typ == MSH_SURF_PLAN && !s->geometry){
+ if (s->Typ == MSH_SURF_PLAN && !s->geometry){
double XP = qp.x();
double YP = qp.y();
double ZP = qp.z();
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index b6f63a77f10823e3cd49280fe46809feca377c80..9b293baf9872c148c87ad7dec6c3077d41ddea4d 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -36,121 +36,6 @@ int Init_Numeric()
return 1;
}
-static double (*nrfuncbrent) (double);
-
-static double fn1(double x, void *params)
-{
- double val = nrfuncbrent(x);
- return val;
-}
-
-static int gsl_min_fminimizer_set_with_values_MOD(gsl_min_fminimizer *s,
- gsl_function *f,
- double x_minimum, double f_minimum,
- double x_lower, double f_lower,
- double x_upper, double f_upper)
-{
- s->function = f;
- s->x_minimum = x_minimum;
- s->x_lower = x_lower;
- s->x_upper = x_upper;
- if (x_lower > x_upper){
- GSL_ERROR ("invalid interval (lower > upper)", GSL_EINVAL);
- }
- s->f_lower = f_lower;
- s->f_upper = f_upper;
- s->f_minimum = f_minimum;
- return (s->type->set) (s->state, s->function,
- x_minimum, f_minimum,
- x_lower, f_lower,
- x_upper, f_upper);
-}
-
-// Returns the minimum betwen ax and cx to a given tolerance tol using
-// brent's method.
-double brent(double ax, double bx, double cx,
- double (*f) (double), double tol, double *xmin)
-{
-#define MAXITER 100
- // The solver can stall at the following internal tolerance in brent - so
- // this is about the lowest possible tolerance.
- if(tol == 0.) tol = 10. * GSL_SQRT_DBL_EPSILON;
- const double tolsq = tol*tol;
- int status;
- int iter = 0;
- double a, b, c; // a < b < c
- double fa, fb, fc;
- const gsl_min_fminimizer_type *T;
- gsl_min_fminimizer *s;
- gsl_function F;
-
- // gsl wants a<b
- b = bx;
- if(ax < cx) {
- a = ax;
- c = cx;
- }
- else {
- a = cx;
- c = ax;
- }
-
- nrfuncbrent = f;
- F.function = &fn1;
- F.params = 0;
-
- fa = f(a);
- fb = f(b);
- fc = f(c);
- if(fb > fa || fb > fc) {
- if(a < c) {
- b = a;
- fb = fa;
- }
- else {
- b = c;
- fb = fc;
- }
- }
-
- // if a-c < tol, return func(b)
- if(gsl_min_test_interval(a, c, tolsq, tol) == GSL_SUCCESS) {
- *xmin = b;
- return (fb);
- }
-
- T = gsl_min_fminimizer_brent;
- s = gsl_min_fminimizer_alloc(T);
- // gsl_min_fminimizer_set(s, &F, b, a, c); // Standard
- // This mod avoids some error checks so we can have a minimum at an extent.
- gsl_min_fminimizer_set_with_values_MOD(s, &F, b, fb, a, fa, c, fc);
-
- do {
- iter++;
- status = gsl_min_fminimizer_iterate(s);
- if(status)
- break; // solver problem
- a = gsl_min_fminimizer_x_lower(s);
- c = gsl_min_fminimizer_x_upper(s);
- // fb = gsl_min_fminimizer_f_minimum(s);
- // fa = gsl_min_fminimizer_f_lower(s);
- // fc = gsl_min_fminimizer_f_upper(s);
- status = gsl_min_test_interval(a, c, tolsq, tol);
- // ... in case curve becomes too flat - not required?
- // std::fabs(fa - fc) < (std::fabs(fb) + tol)*tol
- }
- while(status == GSL_CONTINUE && iter < MAXITER);
- b = gsl_min_fminimizer_x_minimum(s);
-
- if(status != GSL_SUCCESS)
- Msg::Error("Minimization did not converge: f(%g) = %g", b, fn1(b, NULL));
-
- *xmin = b;
- gsl_min_fminimizer_free(s);
- return fn1(b, NULL);
-}
-
-
#define MAX_DIM_NEWT 10
#define MAXITER 100
#define PREC 1.e-8
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 7d3b29e374305f4eb4795decc701d0f21242bb7d..aebc115d68282ea09100e8bc33ff1da85e942c89 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -9,8 +9,6 @@
#include "NumericEmbedded.h"
int Init_Numeric();
-double brent(double ax, double bx, double cx,
- double (*f) (double), double tol, double *xmin);
void newt(double x[], int n, int *check,
void (*func) (int, double[], double[]));
void minimize_N(int N, double (*f) (double*, void *data),