diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 07d9eca7436048a0d640d7bc2552aa8a8a991ef9..3ba644df303b8d7fafa29a4d7dbdd6ff84070b10 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -40,28 +40,7 @@ static void myresid(int N, GEdge *ge, double *u, fullVector<double> &r)
for (int i = 0; i < N - 2; i++) r(i) = L[i + 1] - L[i];
}
-static bool computeEquidistantParameters1(GEdge *ge, double u0, double uN, int N,
- double *u, double underRelax)
-{
- GPoint p0 = ge->point(u0);
- GPoint p1 = ge->point(uN);
- double du = 1. / (N - 1);
- u[0] = u0;
- // printf("starting with %g %g %g\n",p0.x(),p0.y(),u0);
- // printf("ending with %g %g %g\n",p1.x(),p1.y(),uN);
- for (int i = 1; i < N; i++){
- SPoint3 pi (p0.x() + i * du * (p1.x()-p0.x()),
- p0.y() + i * du * (p1.y()-p0.y()),
- p0.z() + i * du * (p1.z()-p0.z()));
- double t;
- GPoint gp = ge->closestPoint(pi, t);
- u[i] = gp.u();
- // printf("going to %g %g u %g\n",pi.x(),pi.y(),gp.u());
- }
- return true;
-}
-
-static bool computeEquidistantParameters0(GEdge *ge, double u0, double uN, int N,
+static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
double *u, double underRelax)
{
const double PRECISION = 1.e-6;
@@ -121,59 +100,11 @@ static bool computeEquidistantParameters0(GEdge *ge, double u0, double uN, int N
return false;
}
-// 1 = geodesics
-static int method_for_computing_intermediary_points = 0;
-static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
- double *u, double underRelax){
- if (method_for_computing_intermediary_points == 0) // use linear abscissa
- return computeEquidistantParameters0(ge,u0,uN,N,u,underRelax);
- else if (method_for_computing_intermediary_points == 1) // use projection
- return computeEquidistantParameters1(ge,u0,uN,N,u,underRelax);
- return false;
-}
-
-static double mylength(GFace *gf, int i, double *u, double *v)
-{
- return gf->length(SPoint2(u[i], v[i]), SPoint2(u[i + 1], v[i + 1]), 10);
-}
-
-static void myresid(int N, GFace *gf, double *u, double *v, fullVector<double> &r)
-{
- double L[100];
- for (int i = 0; i < N - 1; i++) L[i] = mylength(gf, i, u, v);
- for (int i = 0; i < N - 2; i++) r(i) = L[i + 1] - L[i];
-}
-
-static bool computeEquidistantParameters1(GFace *gf, double u0, double uN,
- double v0, double vN, int N,
- double *u, double *v)
-{
- GPoint p0 = gf->point(u0,v0);
- GPoint p1 = gf->point(uN,vN);
- double du = 1. / (N - 1);
- u[0] = u0;
- u[0] = u0;
- v[0] = v0;
- for (int i = 1; i < N; i++){
- SPoint3 pi(p0.x() + i * du * (p1.x()-p0.x()),
- p0.y() + i * du * (p1.y()-p0.y()),
- p0.z() + i * du * (p1.z()-p0.z()));
- SPoint2 t;
- GPoint gp = gf->closestPoint(pi, t);
- u[i] = gp.u();
- v[i] = gp.v();
- }
- return true;
-}
-static bool computeEquidistantParameters0(GFace *gf, double u0, double uN,
+static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
double v0, double vN, int N,
double *u, double *v)
{
- const double PRECISION = 1.e-6;
- const int MAX_ITER = 50;
- const double eps = 1.e-4;
-
double t[100];
// initialize the points by equal subdivision of geodesics
u[0] = u0;
@@ -191,67 +122,6 @@ static bool computeEquidistantParameters0(GFace *gf, double u0, double uN,
return true;
- // create the tangent matrix
- const int M = N - 2;
- fullMatrix<double> J(M, M);
- fullVector<double> DU(M);
- fullVector<double> R(M);
- fullVector<double> Rp(M);
-
- int iter = 1;
-
- while (iter < MAX_ITER){
- iter++;
- myresid(N, gf, u, v, R);
-
- for (int i = 0; i < M; i++){
- t[i + 1] += eps;
- double tempu = u[i + 1];
- double tempv = v[i + 1];
- SPoint2 p = gf->geodesic(SPoint2(u0, v0), SPoint2(uN, vN), t[i + 1]);
- u[i + 1] = p.x();
- v[i + 1] = p.y();
- myresid(N, gf, u, v, Rp);
- for (int j = 0; j < M; j++){
- J(i, j) = (Rp(j) - R(j)) / eps;
- }
- t[i + 1] -= eps;
- u[i + 1] = tempu;
- v[i + 1] = tempv;
- }
-
- if (M == 1)
- DU(0) = R(0) / J(0, 0);
- else
- J.luSolve(R, DU);
-
- for (int i = 0; i < M; i++){
- t[i + 1] -= DU(i);
- SPoint2 p = gf->geodesic(SPoint2(u0, v0), SPoint2(uN, vN), t[i + 1]);
- u[i + 1] = p.x();
- v[i + 1] = p.y();
- }
- double newt_norm = DU.norm();
- if (newt_norm < PRECISION) return true;
- }
- // FAILED, use equidistant in param space
- for (int i = 1; i < N; i++){
- t[i] = (double)i / (N - 1);
- SPoint2 p = gf->geodesic(SPoint2(u0, v0), SPoint2(uN, vN), t[i]);
- u[i] = p.x();
- v[i] = p.y();
- }
- return false;
-}
-
-static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
- double v0, double vN, int N,
- double *u, double *v){
- if (method_for_computing_intermediary_points == 0) // use linear abscissa
- return computeEquidistantParameters0(gf,u0,uN,v0,vN,N,u,v);
- else if (method_for_computing_intermediary_points == 1) // use projection
- return computeEquidistantParameters1(gf,u0,uN,v0,vN,N,u,v);
- return false;
}
static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,