diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 1780772f696ed9f95e4601609f2ba51c30c6f658..1f9d03dd6265f780c9a43c4a313a153a822a1941 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -802,9 +802,9 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax){
wt = _cornerareas[EIdx][j]/_pointareas[vj];
// wt = 1.0;
- _curv1[vj] += wt*c1;
+ _curv1[vj] += wt*c1;
_curv12[vj] += wt*c12;
- _curv2[vj] += wt*c2;
+ _curv2[vj] += wt*c2;
}
} //End of loop over the element (iElem)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index f63a2251a18d9306808c2a257e9be367e528d26b..0e34284b40d664b15ea03199154372aaacbd2511 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -671,7 +671,6 @@ bool GFaceCompound::parametrize() const
fullMatrix<double> Oper(3*allNodes.size(),3*allNodes.size());
_rbf = new GRbf(sizeBox, variableEps, radFunInd, _normals, allNodes, _ordered);
- double t0 = Cpu();
//_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);
@@ -680,12 +679,11 @@ bool GFaceCompound::parametrize() const
//_rbf->RbfLapSurface_local_projection(_rbf->getXYZ(), _rbf->getN(), Oper);
_rbf->solveHarmonicMap(Oper, _ordered, _coords, coordinates);
- double t1 = Cpu();
- printf("Time =%g \n", t1-t0);
+ //_rbf->computeCurvature(coordinates);
+ //exit(1);
//printStuff();
-
}
buildOct();
@@ -1400,7 +1398,7 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
std::cout << "Getting instance of curvature" << std::endl;
curvature.setGModel( model() );
- int computeMax = 1;
+ int computeMax = 0;
curvature.computeCurvature_Rusinkiewicz(computeMax);
curvature.writeToPosFile("curvature.pos");
curvature.writeToVtkFile("curvature.vtk");
@@ -1411,7 +1409,7 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
double c1;
double c2;
curvature.triangleNodalValues(lt->tri,c0, c1, c2, 1);
-
+
double cv = (1-U-V)*c0 + U*c1 + V*c2;
return cv;
@@ -1586,53 +1584,7 @@ void GFaceCompound::secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
{
- if(!oct) parametrize();
-
- //use central differences
-
- //EMI: TODO should take size of two or three triangles
- // double eps = 1e+2;
-
- // double u = param.x();
- // double v = param.y();
- // Pair<SVector3,SVector3> Der_u, Der_ueps, Der_v, Der_veps;
-
- // if(u - eps < 0.0) {
- // Der_u = firstDer(SPoint2(u,v));
- // Der_ueps = firstDer(SPoint2(u+eps,v));
- // }
- // else {
- // Der_u = firstDer(SPoint2(u-eps,v));
- // Der_ueps = firstDer(SPoint2(u,v));
- // }
-
- // if(v - eps < 0.0) {
- // Der_v = firstDer(SPoint2(u,v));
- // Der_veps = firstDer(SPoint2(u,v+eps));
- // }
- // else {
- // Der_v = firstDer(SPoint2(u,v-eps));
- // Der_veps = firstDer(SPoint2(u,v));
- // }
-
- // SVector3 dXdu_u = Der_u.first();
- // SVector3 dXdv_u = Der_u.second();
- // SVector3 dXdu_ueps = Der_ueps.first();
- // SVector3 dXdv_ueps = Der_ueps.second();
- // SVector3 dXdu_v = Der_v.first();
- // SVector3 dXdv_v = Der_v.second();
- // SVector3 dXdu_veps = Der_veps.first();
- // SVector3 dXdv_veps = Der_veps.second();
-
- // double inveps = 1./eps;
- // *dudu = inveps * (dXdu_u - dXdu_ueps) ;
- // *dvdv = inveps * (dXdv_v - dXdv_veps) ;
- // *dudv = inveps * (dXdu_v - dXdu_veps) ;
-
- //printf("der second dudu = %g %g %g \n", dudu->x(), dudu->y(), dudu->z());
- //printf("der second dvdv = %g %g %g \n", dvdv->x(), dvdv->y(), dvdv->z());
- //printf("der second dudv = %g %g %g \n", dudv->x(), dudv->y(), dudv->z());
-
+ if(!oct) parametrize();
Msg::Debug("Computation of the second derivatives is not implemented for compound faces");
}
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index c0b3b3369b2bd1566b7c66a1183c1a2e938e7172..f10e6ab5819f5b515f78be9d494a03405085acb6 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -241,7 +241,6 @@ double OCCFace::curvatures(const SPoint2 ¶m,
double *curvMax,
double *curvMin) const
{
-
const double eps = 1.e-12;
BRepAdaptor_Surface sf(s, Standard_True);
BRepLProp_SLProps prop(sf, 2, eps);
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 17f6682b44fe87c8e2e879ec664874235aaacf14..e484541ac03e23439f87578f2bf1ad6616fc8539 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -483,7 +483,7 @@ SVector3 discreteEdge::firstDer(double par) const
MVertex *vE = lines[iEdge]->getVertex(1);
double dx, dy, dz;
- double dt = 1.0;
+ double dt = 2.0;
dx = (vE->x() - vB->x()) / dt;
dy = (vE->y() - vB->y()) / dt;
dz = (vE->z() - vB->z()) / dt;
@@ -506,7 +506,7 @@ double discreteEdge::curvature(double par) const{
std::cout << "Getting instance of curvature" << std::endl;
curvature.setGModel( model() );
- int computeMax = 1;
+ int computeMax = 0;
curvature.computeCurvature_Rusinkiewicz(computeMax);
curvature.writeToPosFile("curvature.pos");
curvature.writeToVtkFile("curvature.vtk");
@@ -515,7 +515,6 @@ double discreteEdge::curvature(double par) const{
}
curvature.edgeNodalValues(lines[iEdge],c0, c1, 1);
-
double cv = (1-tLoc)*c0 + tLoc*c1;
//printf("curv edge =%g \n", cv);
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index c439a0a63d7b0338c416437666c33db5728c754f..f26124388f9444b9968baa026c6c09a2feb2c8b0 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -330,6 +330,8 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
lc = l1;
}
+ //printf("BGM X Y Z =%g %g %g L4=%g L3=%g L2=%g L1=%g LC=%g LFINAL=%g \n", X, Y, Z, l4, l3, l2, l1, lc , lc* CTX::instance()->mesh.lcFactor);
+
return lc * CTX::instance()->mesh.lcFactor;
}
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 20e0e3752dd7c4328ab26dc216cc1c89ef25edf6..b3506dc9596701a407142508749d0873ffe361d2 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -22,6 +22,7 @@ typedef struct {
static double F_Lc(GEdge *ge, double t)
{
+ //printf("FLC \n");
GPoint p = ge->point(t);
double lc_here;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 349b11ff10bad7630957f18cf8adf4c58e180a5b..f022f3e972e69c16c0cff387188ad358f149b629 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1561,8 +1561,7 @@ void meshGFace::operator() (GFace *gf)
break;
}
- Msg::Info("Meshing surface %d (%s, %s)", gf->tag(),
- gf->getTypeString().c_str(), algo);
+ Msg::Info("Meshing surface %d (%s, %s)", gf->tag(), gf->getTypeString().c_str(), algo);
// compute loops on the fly (indices indicate start and end points
// of a loop; loops are not yet oriented)
diff --git a/benchmarks/curvature/TorusRemesh.geo b/benchmarks/curvature/TorusRemesh.geo
index eb97cef0485425b884b83e19a4b4418ef42649b9..82272f3322c2198ab079bb333fe64863c0499b84 100644
--- a/benchmarks/curvature/TorusRemesh.geo
+++ b/benchmarks/curvature/TorusRemesh.geo
@@ -2,15 +2,15 @@
// The following settings were copied from benchmarks/step/linkrods.geo
// ======================================================================
-Mesh.RemeshAlgorithm=6; //1=MeshAdapt, 5=Delaunay, 6=Frontal
+Mesh.Algorithm=6; //1=MeshAdapt, 5=Delaunay, 6=Frontal, 7=bamg
Mesh.RemeshParametrization= 1; //(0) harmonic (1) conformal
Mesh.RemeshAlgorithm = 1;
-Mesh.CharacteristicLengthFactor=1.2;
+Mesh.CharacteristicLengthFactor=0.4;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
-Mesh.CharacteristicLengthMin = 0.001; //-clmin
-Mesh.CharacteristicLengthMax = 0.5; //-clmax
+Mesh.CharacteristicLengthMin = 0.01; //-clmin
+Mesh.CharacteristicLengthMax = 2.0; //-clmax
Mesh.LcIntegrationPrecision = 1.e-5; //-epslc1d
Mesh.MinimumCirclePoints=15; //default=7
Mesh.CharacteristicLengthExtendFromBoundary=0;