diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index b7c31552fa488724d38fe3b21c938cae42a55674..ba14c57cd44e5a51aac099b6775566a59a8cd5c5 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -42,6 +42,7 @@
#include "multiscaleLaplace.h"
#include "GRbf.h"
#include "Curvature.h"
+#include "MPoint.h"
static void fixEdgeToValue(GEdge *ed, double value, dofManager<double> &myAssembler)
{
@@ -670,9 +671,9 @@ bool GFaceCompound::parametrize() const
fullMatrix<double> Oper(3*allNodes.size(),3*allNodes.size());
_rbf = new GRbf(sizeBox, variableEps, radFunInd, _normals, allNodes, _ordered);
- //_rbf->RbfLapSurface_global_CPM_low(_rbf->getXYZ(), _rbf->getN(), Oper);
+ _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);
+ //_rbf->RbfLapSurface_global_CPM_high(_rbf->getXYZ(), _rbf->getN(), Oper);
//_rbf->RbfLapSurface_local_CPM(false, _rbf->getXYZ(), _rbf->getN(), Oper);
//_rbf->RbfLapSurface_global_projection(_rbf->getXYZ(), _rbf->getN(), Oper);
//_rbf->RbfLapSurface_local_projection(_rbf->getXYZ(), _rbf->getN(), Oper);
@@ -959,7 +960,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
getBoundingEdges();
_type = UNITCIRCLE;
-
+
nbSplit = 0;
fillTris.clear();
}
@@ -1065,6 +1066,7 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
dofManager<double> myAssembler(_lsys);
if(_type == UNITCIRCLE){
+ printf("unit circle \n");
for(unsigned int i = 0; i < _ordered.size(); i++){
MVertex *v = _ordered[i];
const double theta = 2 * M_PI * _coords[i];
@@ -1481,31 +1483,29 @@ 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;
- 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);
- }
+
double U,V;
double par[2] = {par1,par2};
GFaceCompoundTriangle *lt;
getTriangle (par1, par2, <, U,V);
SPoint3 p(3, 3, 0);
- if(!lt){
+ if(!lt && _mapping != RBF){
GPoint gp (p.x(),p.y(),p.z(),this);
gp.setNoSuccess();
return gp;
}
-
- const bool LINEARMESH = true; //false
+ else if (!lt && _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;
+ bool conv = _rbf->UVStoXYZ(par1, par2,x,y,z, dXdu, dXdv);
+ return GPoint(x,y,z);
+ }
if (lt->gf->geomType() != GEntity::DiscreteSurface){
SPoint2 pParam = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
@@ -1513,6 +1513,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
return GPoint(pp.x(),pp.y(),pp.z(),this,par);
}
+ const bool LINEARMESH = true; //false
if(LINEARMESH){
//linear Lagrange mesh
@@ -1580,24 +1581,29 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
if(trivial())
return (*(_compound.begin()))->firstDer(param);
- if (_mapping == RBF){
+ double U, V;
+ GFaceCompoundTriangle *lt;
+ getTriangle(param.x(), param.y(), <, U,V);
+ if(!lt && _mapping != RBF)
+ return Pair<SVector3, SVector3>(SVector3(1, 0, 0), SVector3(0, 1, 0));
+ else if (!lt && _mapping == RBF){
double x, y, z;
SVector3 dXdu, dXdv ;
bool conv = _rbf->UVStoXYZ(param.x(), param.y(), x,y,z, dXdu, dXdv);
return Pair<SVector3, SVector3>(dXdu,dXdv);
}
- double U, V;
- GFaceCompoundTriangle *lt;
- getTriangle(param.x(), param.y(), <, U,V);
- if(!lt)
- return Pair<SVector3, SVector3>(SVector3(1, 0, 0), SVector3(0, 1, 0));
-
double mat[2][2] = {{lt->p2.x() - lt->p1.x(), lt->p3.x() - lt->p1.x()},
{lt->p2.y() - lt->p1.y(), lt->p3.y() - lt->p1.y()}};
double inv[2][2];
- inv2x2(mat,inv);
-
+ double det = inv2x2(mat,inv);
+ if (!det && _mapping == RBF){
+ double x, y, z;
+ SVector3 dXdu, dXdv ;
+ bool conv = _rbf->UVStoXYZ(param.x(), param.y(), x,y,z, dXdu, dXdv);
+ return Pair<SVector3, SVector3>(dXdu,dXdv);
+ }
+
SVector3 dXdxi(lt->v2 - lt->v1);
SVector3 dXdeta(lt->v3 - lt->v1);
@@ -1716,14 +1722,14 @@ void GFaceCompound::getTriangle(double u, double v,
{
double uv[3] = {u, v, 0};
*lt = (GFaceCompoundTriangle*)Octree_Search(uv, oct);
- // if(!(*lt)) {
+ // if(!(*lt)) {
// for(int i=0;i<nbT;i++){
// if(GFaceCompoundInEle (&_gfct[i],uv)){
// *lt = &_gfct[i];
// break;
// }
// }
- // }
+ // }
if(!(*lt)){
return;
}
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 64aefa92e5df859cff75ff074a05df35b9283723..569c6661a6aa31c1be373f2caef47796a3d6e82a 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -195,7 +195,7 @@ double inv2x2(double mat[2][2], double inv[2][2])
inv[1][1] = mat[0][0] * ud;
}
else{
- Msg::Error("Singular matrix");
+ Msg::Error("Singular matrix 2x2");
for(int i = 0; i < 2; i++)
for(int j = 0; j < 2; j++)
inv[i][j] = 0.;
@@ -219,7 +219,7 @@ double inv3x3(double mat[3][3], double inv[3][3])
inv[2][2] = (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]) * ud;
}
else{
- Msg::Error("Singular matrix");
+ Msg::Error("Singular matrix 3x3");
for(int i = 0; i < 3; i++)
for(int j = 0; j < 3; j++)
inv[i][j] = 0.;
diff --git a/benchmarks/curvature/Torus.geo b/benchmarks/curvature/Torus.geo
index d74ba0e05b94c2c649474bbcf41a6fb670144012..5e6adbe6e22f8bf31787844eca473615773cda4b 100644
--- a/benchmarks/curvature/Torus.geo
+++ b/benchmarks/curvature/Torus.geo
@@ -1,3 +1,12 @@
+//Mesh.CharacteristicLengthFactor=0.3;
+Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
+Mesh.CharacteristicLengthMin = 0.1; //-clmin
+Mesh.CharacteristicLengthMax = 2.5; //-clmax
+Mesh.LcIntegrationPrecision=1.e-5; //-epslc1d
+
+Mesh.MinimumCirclePoints=15; //default=7
+Mesh.CharacteristicLengthExtendFromBoundary=0;
+
lc = 0.1;
Point(1) = {2.0,0.0,0.0,lc};
Point(2) = {2.0,1,0.0,lc};