From 275f8bd5b5a024dc242a5a9fd7b51b024bd96663 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Fri, 17 Jun 2011 12:00:14 +0000
Subject: [PATCH] Added new CPM high global method
---
Geo/GFaceCompound.cpp | 2 +
Geo/GRbf.cpp | 215 ++++++++++++++++++++----------------------
Geo/GRbf.h | 7 ++
3 files changed, 109 insertions(+), 115 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 6df0242e27..eee4bbbb40 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -676,6 +676,8 @@ bool GFaceCompound::parametrize() const
fullMatrix<double> Oper(3*allNodes.size(),3*allNodes.size());
_rbf = new GRbf(epsilon, delta, radius, variableEps, radFunInd, _normals, allNodes);
+ //_rbf->test(_rbf->getXYZ());
+
//_rbf->RbfLapSurface_global_CPM_low(_rbf->getXYZ(), _rbf->getN(), Oper);
//_rbf->RbfLapSurface_local_CPM(true, _rbf->getXYZ(), _rbf->getN(), Oper);
_rbf->RbfLapSurface_global_CPM_high(_rbf->getXYZ(), _rbf->getN(), Oper);
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index 56b108c179..cba31a1e63 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -243,9 +243,9 @@ fullMatrix<double> GRbf::generateRbfMat(int p,
if (inUV) eps = epsilonUV;
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);
+ 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);
}
}
@@ -499,60 +499,10 @@ void GRbf::RbfLapSurface_local_CPM(bool isLow,
}
}
-// void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
-// const fullMatrix<double> &normals,
-// fullMatrix<double> &Oper){
-
-// int numNodes = cntrs.size1();
-// int nnTot = 3*numNodes;
-// Oper.resize(nnTot,nnTot);
-
-// fullMatrix<double> sx, sy, sz, sxx, sxy, sxz,syy, syz, szz;
-// fullMatrix<double> Dx, Dy, Dz, Dxx, Dxy, Dxz, Dyy, Dyz, Dzz, Lap, extX, surf;
-
-// setup_level_set(cntrs,normals,extX,surf);
-// if (!isLocal) extendedX = extX;
-// if (!isLocal) surfInterp = surf;
-
-// // Find derivatives of the surface interpolant
-// evalRbfDer(1,extX,cntrs,surf,sx);
-// evalRbfDer(2,extX,cntrs,surf,sy);
-// evalRbfDer(3,extX,cntrs,surf,sz);
-// evalRbfDer(11,extX,cntrs,surf,sxx);
-// evalRbfDer(12,extX,cntrs,surf,sxy);
-// evalRbfDer(13,extX,cntrs,surf,sxz);
-// evalRbfDer(22,extX,cntrs,surf,syy);
-// evalRbfDer(23,extX,cntrs,surf,syz);
-// evalRbfDer(33,extX,cntrs,surf,szz);
-
-// // Finds differentiation matrices
-// RbfOp(1,extX,cntrs,Dx);
-// RbfOp(2,extX,cntrs,Dy);
-// RbfOp(3,extX,cntrs,Dz);
-// RbfOp(11,extX,cntrs,Dxx);
-// RbfOp(12,extX,cntrs,Dxy);
-// RbfOp(13,extX,cntrs,Dxz);
-// RbfOp(22,extX,cntrs,Dyy);
-// RbfOp(23,extX,cntrs,Dyz);
-// RbfOp(33,extX,cntrs,Dzz);
-// RbfOp(222,extX,cntrs,Lap);
-
-// // Fills up the operator matrix
-// for (int i=0;i<numNodes ; ++i){
-// for (int j=0;j<nnTot ; ++j){
-// Oper(i,j) = Lap(i,j);
-// Oper(i+numNodes,j)=sx(i,0)*Dx(i,j)+sy(i,0)*Dy(i,j)+sz(i,0)*Dz(i,j);
-// Oper(i+2*numNodes,j)=sx(i,0)*sx(i,0)*Dxx(i,j)+sy(i,0)*sy(i,0)*Dyy(i,j)+sz(i,0)*sz(i,0)*Dzz(i,j)+2*sx(i,0)*sy(i,0)*Dxy(i,j)+2*sx(i,0)*sz(i,0)*Dxz(i,j)+2*sy(i,0)*sz(i,0)*Dyz(i,j)+(sx(i,0)*sxx(i,0)+sy(i,0)*sxy(i,0)+sz(i,0)*sxz(i,0))*Dx(i,j)+(sx(i,0)*sxy(i,0)+sy(i,0)*syy(i,0)+sz(i,0)*syz(i,0))*Dy(i,j)+(sx(i,0)*sxz(i,0)+sy(i,0)*syz(i,0)+sz(i,0)*szz(i,0))*Dz(i,j);
-// }
-// }
-
-// }
-
//NEW METHOD #1 CPM GLOBAL HIGH
void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
fullMatrix<double> &Oper){
-
int numNodes = cntrs.size1();
int nnTot = 3*numNodes;
Oper.resize(nnTot,nnTot);
@@ -602,65 +552,64 @@ void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
}
-// NEW METHOD #2 CPM GLOBAL HIGH
-// Produces a nxn differentiation matrix (like the projection method)
-// So the local method for this is the local projection method
-// void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
-// const fullMatrix<double> &normals,
-// fullMatrix<double> &Oper){
-
-// int numNodes = cntrs.size1();
-// int nnTot = 3*numNodes;
-// Oper.resize(nnTot,nnTot);
-
-// fullMatrix<double> sx, sy, sz, sxx, sxy, sxz,syy, syz, szz;
-// fullMatrix<double> A, Ax, Ay, Az, Axx, Axy, Axz, Ayy, Ayz, Azz, Alap, AOper, Temp, extX, surf;
-
-// setup_level_set(cntrs,normals,extX,surf);
-// if (!isLocal) extendedX = extX;
-// if (!isLocal) surfInterp = surf;
-
-// // Find derivatives of the surface interpolant
-// evalRbfDer(1,extX,cntrs,surf,sx);
-// evalRbfDer(2,extX,cntrs,surf,sy);
-// evalRbfDer(3,extX,cntrs,surf,sz);
-// evalRbfDer(11,extX,cntrs,surf,sxx);
-// evalRbfDer(12,extX,cntrs,surf,sxy);
-// evalRbfDer(13,extX,cntrs,surf,sxz);
-// evalRbfDer(22,extX,cntrs,surf,syy);
-// evalRbfDer(23,extX,cntrs,surf,syz);
-// evalRbfDer(33,extX,cntrs,surf,szz);
+//NEW METHOD #2 CPM GLOBAL HIGH
+//Produces a nxn differentiation matrix (like the projection method)
+//So the local method for this is the local projection method
+void GRbf::RbfLapSurface_global_CPM_high_2(const fullMatrix<double> &cntrs,
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &Oper){
+
+ int numNodes = cntrs.size1();
+ int nnTot = 3*numNodes;
+ Oper.resize(numNodes,numNodes);
+
+ fullMatrix<double> sx(numNodes,1), sy(numNodes,1), sz(numNodes,1), sxx(numNodes,1), sxy(numNodes,1), sxz(numNodes,1),syy(numNodes,1), syz(numNodes,1), szz(numNodes,1);
+ fullMatrix<double> A(nnTot,nnTot), Ax(numNodes,nnTot), Ay(numNodes,nnTot), Az(numNodes,nnTot), Axx(numNodes,nnTot), Axy(numNodes,nnTot), Axz(numNodes,nnTot), Ayy(numNodes,nnTot), Ayz(numNodes,nnTot), Azz(numNodes,nnTot), Alap(numNodes,nnTot), AOper(nnTot,nnTot), Temp(numNodes,nnTot), extX(nnTot,3), surf(nnTot,1);
+
+ setup_level_set(cntrs,normals,extX,surf);
+ if (!isLocal) extendedX = extX;
+ if (!isLocal) surfInterp = surf;
+
+ // Find derivatives of the surface interpolant
+ evalRbfDer(1,extX,cntrs,surf,sx);
+ evalRbfDer(2,extX,cntrs,surf,sy);
+ evalRbfDer(3,extX,cntrs,surf,sz);
+ evalRbfDer(11,extX,cntrs,surf,sxx);
+ evalRbfDer(12,extX,cntrs,surf,sxy);
+ evalRbfDer(13,extX,cntrs,surf,sxz);
+ evalRbfDer(22,extX,cntrs,surf,syy);
+ evalRbfDer(23,extX,cntrs,surf,syz);
+ evalRbfDer(33,extX,cntrs,surf,szz);
-// // Finds differentiation matrices
-// A=generateRbfMat(0,extX,extX,0);
-// Ax=generateRbfMat(1,extX,cntrs,0);
-// Ay=generateRbfMat(2,extX,cntrs,0);
-// Az=generateRbfMat(3,extX,cntrs,0);
-// Axx=generateRbfMat(11,extX,cntrs,0);
-// Axy=generateRbfMat(12,extX,cntrs,0);
-// Axz=generateRbfMat(13,extX,cntrs,0);
-// Ayy=generateRbfMat(22,extX,cntrs,0);
-// Ayz=generateRbfMat(23,extX,cntrs,0);
-// Azz=generateRbfMat(33,extX,cntrs,0);
-// Alap=generateRbfMat(222,extX,cntrs,0);
-
-// // Fills up the operator matrix
-// for (int i=0;i<numNodes ; ++i){
-// for (int j=0;j<nnTot ; ++j){
-// AOper(i,j) = A(i,j);
-// AOper(i+numNodes,j)=sx(i,0)*Ax(i,j)+sy(i,0)*Ay(i,j)+sz(i,0)*Az(i,j);
-// AOper(i+2*numNodes,j)=sx(i,0)*sx(i,0)*Axx(i,j)+sy(i,0)*sy(i,0)*Ayy(i,j)+sz(i,0)*sz(i,0)*Azz(i,j)+2*sx(i,0)*sy(i,0)*Axy(i,j)+2*sx(i,0)*sz(i,0)*Axz(i,j)+2*sy(i,0)*sz(i,0)*Ayz(i,j)+(sx(i,0)*sxx(i,0)+sy(i,0)*sxy(i,0)+sz(i,0)*sxz(i,0))*Ax(i,j)+(sx(i,0)*sxy(i,0)+sy(i,0)*syy(i,0)+sz(i,0)*syz(i,0))*Ay(i,j)+(sx(i,0)*sxz(i,0)+sy(i,0)*syz(i,0)+sz(i,0)*szz(i,0))*Az(i,j);
-// }
-// }
-// AOper.invertInPlace();
-// Temp.gemm(Alap, AOper, 1.0, 0.0);
-
-// for (int i=0;i<numNodes ; ++i){
-// for (int j=0;j<numNodes ; ++j){
-// Oper(i,j) = Temp(i,j);
-// }
-// }
-// }
+ // Finds differentiation matrices
+ A=generateRbfMat(0,extX,extX,0);
+ Ax=generateRbfMat(1,extX,cntrs,0);
+ Ay=generateRbfMat(2,extX,cntrs,0);
+ Az=generateRbfMat(3,extX,cntrs,0);
+ Axx=generateRbfMat(11,extX,cntrs,0);
+ Axy=generateRbfMat(12,extX,cntrs,0);
+ Axz=generateRbfMat(13,extX,cntrs,0);
+ Ayy=generateRbfMat(22,extX,cntrs,0);
+ Ayz=generateRbfMat(23,extX,cntrs,0);
+ Azz=generateRbfMat(33,extX,cntrs,0);
+ Alap=generateRbfMat(222,extX,cntrs,0);
+
+ // Fills up the operator matrix
+ for (int i=0;i<numNodes ; ++i){
+ for (int j=0;j<nnTot ; ++j){
+ AOper(i,j) = A(i,j);
+ AOper(i+numNodes,j)=sx(i,0)*Ax(i,j)+sy(i,0)*Ay(i,j)+sz(i,0)*Az(i,j);
+ AOper(i+2*numNodes,j)=sx(i,0)*sx(i,0)*Axx(i,j)+sy(i,0)*sy(i,0)*Ayy(i,j)+sz(i,0)*sz(i,0)*Azz(i,j)+2*sx(i,0)*sy(i,0)*Axy(i,j)+2*sx(i,0)*sz(i,0)*Axz(i,j)+2*sy(i,0)*sz(i,0)*Ayz(i,j)+(sx(i,0)*sxx(i,0)+sy(i,0)*sxy(i,0)+sz(i,0)*sxz(i,0))*Ax(i,j)+(sx(i,0)*sxy(i,0)+sy(i,0)*syy(i,0)+sz(i,0)*syz(i,0))*Ay(i,j)+(sx(i,0)*sxz(i,0)+sy(i,0)*syz(i,0)+sz(i,0)*szz(i,0))*Az(i,j);
+ }
+ }
+ AOper.invertInPlace();
+ Alap.mult(AOper,Temp);
+ for (int i=0;i<numNodes ; ++i){
+ for (int j=0;j<numNodes ; ++j){
+ Oper(i,j) = Temp(i,j);
+ }
+ }
+}
void GRbf::RbfLapSurface_global_CPM_low(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
@@ -719,7 +668,6 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
std::vector<MVertex*> ordered,
std::vector<double> coords,
std::map<MVertex*, SPoint3> &rbf_param){
-
int nb = Oper.size2();
UV.resize(nb,2);
@@ -945,9 +893,9 @@ bool GRbf::UVStoXYZ(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,xyz_local,nodes_eval,u_vec,u_temp);
@@ -970,3 +918,40 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
return converged;
}
+
+void GRbf::test(const fullMatrix<double> &cntrs) {
+
+ fullMatrix<double> theta(100,1),f_theta(100,1),f_prime_theta(100,1),f_prime_approx(100,1),rbfInvA(100,100),rbfMatB(100,100),D(100,100);
+ for (int j = 0; j< 100;j++){
+ theta(j,0) = 2*3.14*j/100;
+ f_theta(j,0) = cos(theta(j,0));
+ f_prime_theta(j,0) = -sin(theta(j,0));
+ }
+
+ for (int i = 0; i < 100; i++) {
+ for (int j = 0; j < 100; j++) {
+ double dx = theta(i,0)-theta(j,0);
+ double dy = 0;
+ double dz = 0;
+ rbfInvA(i,j) = evalRadialFnDer(0,dx,dy,dz,2);
+ }
+ }
+
+ rbfInvA.invertInPlace();
+
+ for (int i = 0; i < 100; i++) {
+ for (int j = 0; j < 100; j++) {
+ double dx = theta(i,0)-theta(j,0);
+ double dy = 0;
+ double dz = 0;
+ rbfMatB(i,j) = evalRadialFnDer(1,dx,dy,dz,2);
+ }
+ }
+
+ D.gemm(rbfMatB, rbfInvA, 1.0, 0.0);
+ f_prime_approx.gemm(D,f_theta, 1.0, 0.0);
+
+ f_prime_approx.print("f_prime_approx");
+ f_prime_theta.print("f_prime_theta");
+
+}
diff --git a/Geo/GRbf.h b/Geo/GRbf.h
index 39a0d57205..24087590aa 100644
--- a/Geo/GRbf.h
+++ b/Geo/GRbf.h
@@ -124,6 +124,11 @@ class GRbf {
void RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
fullMatrix<double> &D);
+
+ // Second method that Finds global surface differentiation matrix (method CPMH)
+ void RbfLapSurface_global_CPM_high_2(const fullMatrix<double> &cntrs,
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &D);
// Calculates the curvature of a surface at a certain node
@@ -147,6 +152,8 @@ class GRbf {
inline const fullMatrix<double> getXYZ() {return centers;};
inline const fullMatrix<double> getN() {return normals;};
+ void test(const fullMatrix<double> &cntrs);
+
};
#endif
--
GitLab