From 3b2154ad5c43d17b62ebb0c8fe1dfc9185e173e8 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Fri, 17 Jun 2011 07:28:49 +0000
Subject: [PATCH] CPM high is working fine now
---
Geo/GRbf.cpp | 73 +++++++++++++++++++++-------------------------------
Geo/GRbf.h | 2 +-
2 files changed, 30 insertions(+), 45 deletions(-)
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index b07b0576ab..56b108c179 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -548,7 +548,6 @@ 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,
@@ -576,8 +575,7 @@ void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
evalRbfDer(23,extX,cntrs,surf,syz);
evalRbfDer(33,extX,cntrs,surf,szz);
- szz.print("szz");
- // Finds differentiation matrices
+ // Finds differentiation matrices
A=generateRbfMat(0,extX,extX,0);
Ax=generateRbfMat(1,extX,cntrs,0);
Ay=generateRbfMat(2,extX,cntrs,0);
@@ -591,6 +589,7 @@ void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
Alap=generateRbfMat(222,extX,cntrs,0);
// Fills up the operator matrix
+ AOper.resize(nnTot, nnTot);
for (int i=0;i<numNodes ; ++i){
for (int j=0;j<nnTot ; ++j){
AOper(i,j) = Alap(i,j);
@@ -770,10 +769,13 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
}
-void GRbf::UVStoXYZ_global(const double u_eval, const double v_eval,
+bool GRbf::UVStoXYZ_global(const double u_eval, const double v_eval,
double &XX, double &YY, double &ZZ,
- SVector3 &dXdu, SVector3& dXdv)
+ SVector3 &dXdu, SVector3& dXdv)
{
+
+ bool converged = true;
+
int numNodes = 3*nbNodes; //WARNING: THIS IS KO FOR PROJECTION ??
double norm_s = 0.0;
fullMatrix<double> u_temp(1,3);
@@ -787,7 +789,7 @@ void GRbf::UVStoXYZ_global(const double u_eval, const double v_eval,
u_vec(i,1) = UV(i,1);
u_vec(i,2) = surfInterp(i,0);
}
-
+
// u_vec_eval = [u_eval v_eval s_eval]
fullMatrix<double> u_vec_eval(1, 3);
u_vec_eval(0,0) = u_eval;
@@ -797,15 +799,13 @@ void GRbf::UVStoXYZ_global(const double u_eval, const double v_eval,
//Start with the closest point (u,v)
#if defined (HAVE_ANN)
double uvw[3] = { u_eval, v_eval, 0.0 };
- ANNidxArray index1;
- ANNdistArray dist1;
- UVkdtree->annkSearch(uvw, 1, index1, dist1);
+ UVkdtree->annkSearch(uvw, num_neighbours, index, dist);
#endif
- nodes_eval(0,0) = extendedX(index1[0],0);
- nodes_eval(0,1) = extendedX(index1[0],1);
- nodes_eval(0,2) = extendedX(index1[0],2);
- u_temp(0,0) = UV(index1[0],0);
- u_temp(0,1) = UV(index1[0],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);
+ u_temp(0,0) = UV(index[0],0);
+ u_temp(0,1) = UV(index[0],1);
u_temp(0,2) = 0.0;
int incr =0;
@@ -839,8 +839,10 @@ void GRbf::UVStoXYZ_global(const double u_eval, const double v_eval,
incr++;
}
- if (norm_s < 5 )
+ if (norm_s < 5 ){
printf("Newton not converged for point (uv)=(%g,%g -> norm_s =%g )\n", u_eval, v_eval, norm_s);
+ converged = false;
+ }
XX = nodes_eval(0,0);
YY = nodes_eval(0,1);
@@ -849,25 +851,14 @@ void GRbf::UVStoXYZ_global(const double u_eval, const double v_eval,
dXdu = SVector3(Jac(0,0), Jac(1,0), Jac(2,0));
dXdv = SVector3(Jac(0,1), Jac(1,1), Jac(2,1));
- delete [] index1;
- delete [] dist1;
-
+ return converged;
+
}
bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
double &XX, double &YY, double &ZZ,
SVector3 &dXdu, SVector3& dXdv){
-
- fullMatrix<double> vec(3*nbNodes,3);
- for (int i=0; i<3*nbNodes;i++){
- vec(i,0) = normals(i,0);
- vec(i,1) = normals(i,1);
- vec(i,2) = normals(i,2);
- printf("%g %g %g\n", vec(i,0),vec(i,1),vec(i,2));
- }
- exit(1);
-
//Finds the U,V,S (in the 0-level set) that are the 'num_neighbours' closest to u_eval and v_eval.
//Thus in total, we're working with '3*num_neighbours' nodes
//Say that the vector 'index' gives us the indices of the closest points
@@ -876,9 +867,7 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
#if defined (HAVE_ANN)
double uvw[3] = { u_eval, v_eval, 0.0 };
- ANNidxArray index2 = new ANNidx[num_neighbours];
- ANNdistArray dist2 = new ANNdist[num_neighbours];
- UVkdtree->annkSearch(uvw, num_neighbours, index2, dist2);
+ UVkdtree->annkSearch(uvw, num_neighbours, index, dist);
#endif
fullMatrix<double> ux(1,3), uy(1,3), uz(1,3), u_temp(1,3);
@@ -886,7 +875,6 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
fullMatrix<double> u_vec(3*num_neighbours,3);
fullMatrix<double> xyz_local(3*num_neighbours,3);
-
// u_vec = [u v s] : These are the u,v,s that are on the surface *and* outside of it also!
for (int i = 0; i < num_neighbours; i++){
@@ -905,17 +893,17 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
//THE LOCAL XYZ
- xyz_local(i,0) = extendedX(index2[i],0);
- xyz_local(i,1) = extendedX(index2[i],1);
- xyz_local(i,2) = extendedX(index2[i],2);
+ xyz_local(i,0) = extendedX(index[i],0);
+ xyz_local(i,1) = extendedX(index[i],1);
+ xyz_local(i,2) = extendedX(index[i],2);
- xyz_local(i+num_neighbours,0) = extendedX(index2[i]+nbNodes,0);
- xyz_local(i+num_neighbours,1) = extendedX(index2[i]+nbNodes,1);
- xyz_local(i+num_neighbours,2) = extendedX(index2[i]+nbNodes,2);
+ xyz_local(i+num_neighbours,0) = extendedX(index[i]+nbNodes,0);
+ xyz_local(i+num_neighbours,1) = extendedX(index[i]+nbNodes,1);
+ xyz_local(i+num_neighbours,2) = extendedX(index[i]+nbNodes,2);
- xyz_local(i+2*num_neighbours,0) = extendedX(index2[i]+2*nbNodes,0);
- xyz_local(i+2*num_neighbours,1) = extendedX(index2[i]+2*nbNodes,1);
- xyz_local(i+2*num_neighbours,2) = extendedX(index2[i]+2*nbNodes,2);
+ xyz_local(i+2*num_neighbours,0) = extendedX(index[i]+2*nbNodes,0);
+ xyz_local(i+2*num_neighbours,1) = extendedX(index[i]+2*nbNodes,1);
+ xyz_local(i+2*num_neighbours,2) = extendedX(index[i]+2*nbNodes,2);
}
// u_vec_eval = [u_eval v_eval s_eval]
@@ -978,9 +966,6 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
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));
-
- delete [] index2;
- delete [] dist2;
return converged;
}
diff --git a/Geo/GRbf.h b/Geo/GRbf.h
index cad135e2b4..39a0d57205 100644
--- a/Geo/GRbf.h
+++ b/Geo/GRbf.h
@@ -132,7 +132,7 @@ class GRbf {
const fullMatrix<double> &node,
fullMatrix<double> &curvature);
- void UVStoXYZ_global(const double u_eval, const double v_eval,
+ bool UVStoXYZ_global(const double u_eval, const double v_eval,
double &XX, double &YY, double &ZZ,
SVector3 &dXdu, SVector3& dxdv);
--
GitLab