diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 07ffc98953805acd2e117a9dfe048dc2b0236beb..f703d10c682c5f7701bbfb0eeacbe409db1f72cc 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -665,7 +665,7 @@ bool GFaceCompound::parametrize() const
int variableEps = 0;
int radFunInd = 1; //MQ RBF
double delta = 0.33*getDistMin();
- double epsilon = 0.15*(allNodes.size()/150.0)/delta; //max(2.5, 0.5*(nbNodes/150.0)/dist_min);
+ double epsilon = 0.5/getDistMin(); //0.15*(allNodes.size()/150.0)/delta; //max(2.5, 0.5*(nbNodes/150.0)/dist_min);
double radius= 3.*getSizeH()/sqrt(allNodes.size());
Msg::Info("*****************************************");
@@ -685,8 +685,7 @@ bool GFaceCompound::parametrize() const
_rbf->solveHarmonicMap(Oper, _ordered, _coords, coordinates);
printStuff();
- exit(1);
-
+
}
buildOct();
@@ -780,18 +779,18 @@ double GFaceCompound::getSizeH() const
double GFaceCompound::getDistMin() const
{
- double dist_min = 1.e6;
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; itv++){
- for(std::set<MVertex *>::iterator itv2 = allNodes.begin(); itv2 !=allNodes.end() ; itv2++){
- MVertex *v1 = *itv;
- MVertex *v2 = *itv2;
- double dist = sqrt((v1->x()-v2->x())*(v1->x()-v2->x())+(v1->y()-v2->y())*(v1->y()-v2->y())
- +(v1->z()-v2->z())*(v1->z()-v2->z()));
- if (dist<dist_min && dist != 0.0) dist_min = dist;
- }
+ double dist_min = 1.e6;
+ double tol = 1.e-8;
+ for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; itv++){
+ for(std::set<MVertex *>::iterator itv2 = allNodes.begin(); itv2 !=allNodes.end() ; itv2++){
+ MVertex *v1 = *itv;
+ MVertex *v2 = *itv2;
+ double dist = sqrt((v1->x()-v2->x())*(v1->x()-v2->x())+(v1->y()-v2->y())*(v1->y()-v2->y())
+ +(v1->z()-v2->z())*(v1->z()-v2->z()));
+ if (dist<dist_min && dist > tol) dist_min = dist;
}
-
- return dist_min;
+ }
+ return dist_min;
}
double GFaceCompound::getSizeBB(const std::list<GEdge* > &elist) const
{
@@ -1446,11 +1445,15 @@ GPoint GFaceCompound::point(double par1, double par2) const
if(!oct) parametrize();
if (_mapping == RBF){
+ if (fabs(par1) > 1 || fabs(par2) > 1){
+ GPoint gp (3,3,0,this);
+ gp.setNoSuccess();
+ return gp;
+ }
double x, y, z;
SVector3 dXdu, dXdv;
- _rbf->UVStoXYZ(par1, par2,x,y,z, dXdu, dXdv);
- //printf("XYZ =%g %g %g \n", x,y,z);
- //exit(1);
+ bool conv = _rbf->UVStoXYZ(par1, par2,x,y,z, dXdu, dXdv);
+ if (!conv) printf("UV=%g %g \n", par1, par2);
return GPoint(x,y,z);
}
@@ -1543,7 +1546,7 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
if (_mapping == RBF){
double x, y, z;
SVector3 dXdu, dXdv ;
- _rbf->UVStoXYZ(param.x(), param.y(), x,y,z, dXdu, dXdv);
+ bool conv = _rbf->UVStoXYZ(param.x(), param.y(), x,y,z, dXdu, dXdv);
return Pair<SVector3, SVector3>(dXdu,dXdv);
}
@@ -1892,6 +1895,8 @@ double GFaceCompound::checkAspectRatio() const
void GFaceCompound::coherencePatches() const
{
+
+ if (_mapping == RBF) return;
Msg::Info("Re-orient all %d compound patches normals coherently", _compound.size());
std::map<MEdge, std::set<MElement*>, Less_Edge > edge2elems;
@@ -2033,7 +2038,7 @@ int GFaceCompound::genusGeom() const
void GFaceCompound::printStuff(int iNewton) const
{
- if( !CTX::instance()->mesh.saveAll) return;
+ //if( !CTX::instance()->mesh.saveAll) return;
std::list<GFace*>::const_iterator it = _compound.begin();
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index 4c6da4899ba470fe3d9c929774c76215eb796e96..0f55aa5ba079151b59162ff01f46ffe9293bab72 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -45,7 +45,7 @@ static int SphereInEle(void *a, double*c){
GRbf::GRbf (double eps, double del, double rad, int variableEps, int rbfFun,
std::map<MVertex*, SVector3> _normals, std::set<MVertex *> allNodes)
- : ep_scalar(eps), delta(del), radius (rad), variableShapeParam(variableEps),
+ : epsilonXYZ(eps), epsilonUV(eps), delta(del), radius (rad), variableShapeParam(variableEps),
radialFunctionIndex (rbfFun), XYZkdtree(0)
{
@@ -230,20 +230,23 @@ double GRbf::evalRadialFnDer (int p, double dx, double dy, double dz, double ep)
}
fullMatrix<double> GRbf::generateRbfMat(int p,
- const fullMatrix<double> &nodes1,
- const fullMatrix<double> &nodes2) {
+ const fullMatrix<double> &nodes1,
+ const fullMatrix<double> &nodes2,
+ int inUV) {
int m = nodes2.size1();
int n = nodes1.size1();
fullMatrix<double> rbfMat(m,n);
- fullVector<double > epsilon;
- computeEpsilon(nodes1, epsilon);
+ //fullVector<double > epsilon;
+ //computeEpsilon(nodes1, epsilon, inUN);
+ double eps = epsilonXYZ;
+ 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);
- rbfMat(i,j) = evalRadialFnDer(p,dx,dy,dz,epsilon(j));
+ rbfMat(i,j) = evalRadialFnDer(p,dx,dy,dz,eps);
}
}
@@ -254,28 +257,28 @@ fullMatrix<double> GRbf::generateRbfMat(int p,
void GRbf::RbfOp(int p,
const fullMatrix<double> &cntrs,
const fullMatrix<double> &nodes,
- fullMatrix<double> &D) {
+ fullMatrix<double> &D, int inUV) {
fullMatrix<double> rbfInvA, rbfMatB;
D.resize(nodes.size1(), cntrs.size1());
if (isLocal){
- rbfInvA = generateRbfMat(0,cntrs,cntrs);
+ rbfInvA = generateRbfMat(0,cntrs,cntrs, inUV);
rbfInvA.invertInPlace();
}
else{
- if (cntrs.size1() == nbNodes )
+ if (cntrs.size1() == nbNodes && !inUV)
rbfInvA = matAInv;
- else if (cntrs.size1() == 3*nbNodes )
+ else if (cntrs.size1() == 3*nbNodes && !inUV)
rbfInvA = matAInv_nn;
else{
- rbfInvA = generateRbfMat(0,cntrs,cntrs);
+ rbfInvA = generateRbfMat(0,cntrs,cntrs, inUV);
rbfInvA.invertInPlace();
}
}
- rbfMatB = generateRbfMat(p,cntrs,nodes);
+ rbfMatB = generateRbfMat(p,cntrs,nodes, inUV);
D.gemm(rbfMatB, rbfInvA, 1.0, 0.0);
}
@@ -284,19 +287,22 @@ void GRbf::evalRbfDer(int p,
const fullMatrix<double> &cntrs,
const fullMatrix<double> &nodes,
const fullMatrix<double> &fValues,
- fullMatrix<double> &fApprox) {
+ fullMatrix<double> &fApprox,
+ int inUV) {
fApprox.resize(nodes.size1(),fValues.size2());
fullMatrix<double> D;
- RbfOp(p,cntrs,nodes,D);
+ RbfOp(p,cntrs,nodes,D, inUV);
fApprox.gemm(D,fValues, 1.0, 0.0);
}
-void GRbf::computeEpsilon(const fullMatrix<double> &cntrs, fullVector<double> &epsilon){
+void GRbf::computeEpsilon(const fullMatrix<double> &cntrs, fullVector<double> &epsilon,
+ int inUV){
epsilon.resize(nbNodes*3);
- epsilon.setAll(ep_scalar);
+ if (inUV) epsilon.setAll(epsilonUV);
+ epsilon.setAll(epsilonXYZ);
if (variableShapeParam) {
@@ -342,7 +348,7 @@ void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
for (int j=0;j<3 ; ++j){
level_set_nodes(i,j) = cntrs(i,j);
level_set_nodes(i+numNodes,j) = cntrs(i,j)-delta*normals(i,j);
- level_set_nodes(i+2*numNodes,j) = cntrs(i,j)+delta*normals(i,j);
+ level_set_nodes(i+2*numNodes,j) = cntrs(i,j)+delta*normals(i,j);
}
level_set_funvals(i,0) = 0.0;
level_set_funvals(i+numNodes,0) = -1;
@@ -603,20 +609,12 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
int nb = Oper.size2();
UV.resize(nb,2);
- fullMatrix<double> dirichletBC(ordered.size(),2);
- for(unsigned int i = 0; i < ordered.size(); i++){
- MVertex *v = ordered[i];
- std::map<MVertex*, int>::iterator itm = _mapV.find(v);
- const double theta = 2 * M_PI * coords[i];
- dirichletBC(i,0) = itm->second;
- dirichletBC(i,1) = theta;
- }
-
fullMatrix<double> F(nb,2);
F.scale(0.0);
for (int i=0; i < ordered.size(); i++){
- int iFix = (int)dirichletBC(i,0);
- double theta = dirichletBC(i,1);
+ std::map<MVertex*, int>::iterator itm = _mapV.find(ordered[i]);
+ double theta = 2 * M_PI * coords[i];
+ int iFix = itm->second;
for (int j=0;j<nb ; ++j) Oper(iFix,j) = 0.0;
Oper(iFix,iFix) = 1.0;
F(iFix,0) = cos(theta);
@@ -626,17 +624,26 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
Oper.invertInPlace();
Oper.mult(F,UV);
+ //ANN UVtree + dist_min
+ double dist_min = 1.e6;
#if defined (HAVE_ANN)
UVnodes = annAllocPts(nbNodes, 3);
for(int i = 0; i < nbNodes; i++){
UVnodes[i][0] = UV(i,0);
UVnodes[i][1] = UV(i,1);
UVnodes[i][2] = 0.0;
+ for(int j = i+1; j < nbNodes; j++){
+ double dist = sqrt((UV(i,0)-UV(j,0))*(UV(i,0)-UV(j,0))+
+ (UV(i,1)-UV(j,1))*(UV(i,1)-UV(j,1)));
+ if (dist<dist_min) dist_min = dist;
+ }
}
UVkdtree = new ANNkd_tree(UVnodes, nbNodes, 3);
index = new ANNidx[num_neighbours];
dist = new ANNdist[num_neighbours];
#endif
+ epsilonUV = 0.5/dist_min;
+ deltaUV = 0.6*dist_min;
//fill rbf_param
std::map<MVertex*, int>::iterator itm = _mapV.begin();
@@ -728,7 +735,7 @@ void GRbf::UVStoXYZ_global(const double u_eval, const double v_eval,
}
-void GRbf::UVStoXYZ(const double u_eval, const double v_eval,
+bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
double &XX, double &YY, double &ZZ,
SVector3 &dXdu, SVector3& dXdv){
@@ -736,6 +743,7 @@ void GRbf::UVStoXYZ(const double u_eval, const double 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
//TO FILL IN
+ bool converged = true;
#if defined (HAVE_ANN)
double uvw[3] = { u_eval, v_eval, 0.0 };
@@ -753,15 +761,16 @@ void GRbf::UVStoXYZ(const double u_eval, const double v_eval,
//THE LOCAL UVS
u_vec(i,0) = UV(index[i],0);
u_vec(i,1) = UV(index[i],1);
- u_vec(i,2) = surfInterp(index[i],0);
+ u_vec(i,2) = 0.0;
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);
+ 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);
+ u_vec(i+2*num_neighbours,2) = surfInterp(index[i]+2*nbNodes,0)*deltaUV;
+
//THE LOCAL XYZ
xyz_local(i,0) = extendedX(index[i],0);
@@ -784,7 +793,13 @@ void 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);
+ 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);
u_temp(0,0) = u_eval;
u_temp(0,1) = v_eval;
@@ -793,7 +808,7 @@ void GRbf::UVStoXYZ(const double u_eval, const double v_eval,
int incr = 0;
double norm_s = 0.0;
fullMatrix<double> Jac(3,3);
- while(norm_s <5 && incr < 5){
+ while(norm_s <5 && incr < 10){
// Find the entries of the m Jacobians
evalRbfDer(1,xyz_local,nodes_eval,u_vec,ux);
@@ -822,14 +837,16 @@ void GRbf::UVStoXYZ(const double u_eval, const double v_eval,
incr++;
}
- if (norm_s < 5 )
- printf("Newton not converged for point (uv)=(%g,%g -> norm_s =%g )\n", u_eval, v_eval, norm_s);
-
+ 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;
+ }
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;
}
diff --git a/Geo/GRbf.h b/Geo/GRbf.h
index dc5ecf776a31586d7f488746916aa50f14c12344..cad135e2b427ede1cd9945311537157929cf6220 100644
--- a/Geo/GRbf.h
+++ b/Geo/GRbf.h
@@ -32,8 +32,10 @@ class GRbf {
int nn;
int num_neighbours;
- double ep_scalar; // Shape parameter
+ double epsilonXYZ; // Shape parameter
+ double epsilonUV; // Shape parameter
double delta; //offset level set
+ double deltaUV; //offset level set
double radius;
int variableShapeParam; // 1 if one chooses epsilon to vary spatially, 0 if one chooses it to be constant
int radialFunctionIndex; // Index corresponding to the radial function used (0 - GA,1 - MQ, ... )
@@ -79,22 +81,23 @@ class GRbf {
//(p)th derivative of the radial function w.r.t. the (q)th variable
fullMatrix<double> generateRbfMat(int p,
const fullMatrix<double> &nodes1,
- const fullMatrix<double> &nodes2);
+ const fullMatrix<double> &nodes2,
+ int inUV=0);
// Computes the interpolation(p==0) or the derivative (p!=0) operator(mxn) (n:number of centers, m: number of evaluation nodes)
void RbfOp(int p, // (p)th derivatives
const fullMatrix<double> &cntrs,
const fullMatrix<double> &nodes,
- fullMatrix<double> &D);
+ fullMatrix<double> &D, int inUV=0);
// Computes the interpolant(p==0) or the derivative (p!=0) of the function values entered and evaluates it at the new nodes
void evalRbfDer(int p, // (p)th derivatives
const fullMatrix<double> &cntrs,
const fullMatrix<double> &nodes,
const fullMatrix<double> &fValues,
- fullMatrix<double> &fApprox);
+ fullMatrix<double> &fApprox, int inUV=0);
- void computeEpsilon(const fullMatrix<double> &cntrs, fullVector<double> &epsilon);
+ void computeEpsilon(const fullMatrix<double> &cntrs, fullVector<double> &epsilon, int inUV=0);
// Finds surface differentiation matrix using the LOCAL projection method
void RbfLapSurface_local_projection(const fullMatrix<double> &cntrs,
@@ -129,12 +132,11 @@ class GRbf {
const fullMatrix<double> &node,
fullMatrix<double> &curvature);
- virtual void UVStoXYZ_global(const double u_eval,
- const double v_eval,
- double &XX, double &YY, double &ZZ,
- SVector3 &dXdu, SVector3& dxdv);
- virtual void UVStoXYZ(const double u_eval,
- const double v_eval,
+ void UVStoXYZ_global(const double u_eval, const double v_eval,
+ double &XX, double &YY, double &ZZ,
+ SVector3 &dXdu, SVector3& dxdv);
+
+ bool UVStoXYZ(const double u_eval, const double v_eval,
double &XX, double &YY, double &ZZ,
SVector3 &dXdu, SVector3& dxdv);