diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index cba31a1e63ef217e0523499b83bc8d57e2dcf26d..4b60a3297ae2250b39281ba94be4553a48543ce6 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -237,6 +237,8 @@ fullMatrix<double> GRbf::generateRbfMat(int p,
int m = nodes2.size1();
int n = nodes1.size1();
fullMatrix<double> rbfMat(m,n);
+
+
//fullVector<double > epsilon;
//computeEpsilon(nodes1, epsilon, inUN);
double eps = epsilonXYZ;
@@ -250,6 +252,31 @@ fullMatrix<double> GRbf::generateRbfMat(int p,
}
}
+
+
+ // //////////////////////
+ // double min_dist = 10;
+ // for (int i = 0; i < m; i++) {
+ // for (int j = i+1; j < n; j++) {
+ // double dx = nodes1(i,0)-nodes1(j,0);
+ // double dy = nodes1(i,1)-nodes1(j,1);
+ // double dz = nodes1(i,2)-nodes1(j,2);
+ // double dist_node = sqrt(dx*dx+dy*dy+dz*dz);
+ // if (dist_node<min_dist) min_dist = dist_node;
+ // }
+ // }
+ // double eps = 0.5/min_dist;
+ // for (int i = 0; i < m; i++) {
+ // for (int j = 0; j < n; j++) {
+ // double dx = nodes2(i,0)-nodes1(j,0);
+ // double dy = nodes2(i,1)-nodes1(j,1);
+ // double dz = nodes2(i,2)-nodes1(j,2);
+ // rbfMat(i,j) = evalRadialFnDer(p,dx,dy,dz,eps);
+ // }
+ // }
+ // //////////////////////
+
+
return rbfMat;
}
@@ -502,7 +529,7 @@ void GRbf::RbfLapSurface_local_CPM(bool isLow,
//NEW METHOD #1 CPM GLOBAL HIGH
void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
- fullMatrix<double> &Oper){
+ fullMatrix<double> &Oper){
int numNodes = cntrs.size1();
int nnTot = 3*numNodes;
Oper.resize(nnTot,nnTot);
@@ -681,12 +708,14 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
Oper(iFix,iFix) = 1.0;
F(iFix,0) = cos(theta);
F(iFix,1) = sin(theta);
+
}
Oper.invertInPlace();
Oper.mult(F,UV);
//ANN UVtree + dist_min
+
double dist_min = 1.e6;
#if defined (HAVE_ANN)
UVnodes = annAllocPts(nbNodes, 3);
@@ -775,9 +804,9 @@ bool GRbf::UVStoXYZ_global(const double u_eval, const double v_eval,
for (int j = 0; j< 3;j++)
nodes_eval(0,j) = nodes_eval(0,j)
- - Jac(j,0)*(u_vec_eval(0,0) - u_temp(0,0))
- - Jac(j,1)*(u_vec_eval(0,1) - u_temp(0,1))
- - Jac(j,2)*(0.0 - u_temp(0,2));
+ - Jac(j,0)*(-u_vec_eval(0,0) + u_temp(0,0))
+ - Jac(j,1)*(-u_vec_eval(0,1) + u_temp(0,1))
+ - Jac(j,2)*(-0.0 + u_temp(0,2));
// Find the corresponding u, v, s
evalRbfDer(0,extendedX,nodes_eval,u_vec,u_temp);
@@ -833,11 +862,11 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
u_vec(i+num_neighbours,0) = UV(index[i]+nbNodes,0);
u_vec(i+num_neighbours,1) = UV(index[i]+nbNodes,1);
- u_vec(i+num_neighbours,2) = surfInterp(index[i]+nbNodes,0)*deltaUV;
+ u_vec(i+num_neighbours,2) = surfInterp(index[i]+nbNodes,0);//*deltaUV;
u_vec(i+2*num_neighbours,0) = UV(index[i]+2*nbNodes,0);
u_vec(i+2*num_neighbours,1) = UV(index[i]+2*nbNodes,1);
- u_vec(i+2*num_neighbours,2) = surfInterp(index[i]+2*nbNodes,0)*deltaUV;
+ u_vec(i+2*num_neighbours,2) = surfInterp(index[i]+2*nbNodes,0);//*deltaUV;
//THE LOCAL XYZ
@@ -861,13 +890,13 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
u_vec_eval(0,2) = 0.0;
//we will use a local interpolation to find the corresponding XYZ point to (u_eval,v_eval).
- evalRbfDer(0, u_vec, u_vec_eval,xyz_local, nodes_eval, 1);
+ //evalRbfDer(0, u_vec, u_vec_eval,xyz_local, nodes_eval, 1);
//printf("nodes eval =%g %g %g \n", nodes_eval(0,0), nodes_eval(0,1), nodes_eval(0,2));
//exit(1);
- // nodes_eval(0,0) = extendedX(index[0],0);
- // nodes_eval(0,1) = extendedX(index[0],1);
- // nodes_eval(0,2) = extendedX(index[0],2);
+ nodes_eval(0,0) = extendedX(index[0],0);
+ nodes_eval(0,1) = extendedX(index[0],1);
+ nodes_eval(0,2) = extendedX(index[0],2);
u_temp(0,0) = u_eval;
u_temp(0,1) = v_eval;
@@ -902,19 +931,24 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
double normSq = sqrt(u_temp(0,2)*u_temp(0,2));
norm_s = fabs(log10(normSq));
-
incr++;
}
if (norm_s < 5 ){
printf("Newton not converged for point (uv)=(%g,%g -> norm_s =%g ) XYZ =%g %g %g \n", u_eval, v_eval, norm_s, nodes_eval(0,0), nodes_eval(0,1), nodes_eval(0,2));
converged = false;
+
+
+ // if Newton diverges, call the routine again with
+ // num_neighbors = num_neighbors + 5;
+ // UVStoXYZ(const double u_eval, const double v_eval,
+ // double &XX, double &YY, double &ZZ,
+ // SVector3 &dXdu, SVector3& dXdv)
}
XX = nodes_eval(0,0);
YY = nodes_eval(0,1);
ZZ = nodes_eval(0,2);
dXdu = SVector3(Jac(0,0), Jac(1,0), Jac(2,0));
dXdv = SVector3(Jac(0,1), Jac(1,1), Jac(2,1));
-
return converged;
}