diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index ef9a2f23046f9f090744c2c546709653b1cd5289..fddcdda40b4d1fba779d7b1b9afcafc20e8896b6 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -19,6 +19,7 @@
#include "Context.h"
#include "meshGFaceOptimize.h"
#include "meshGFaceLloyd.h"
+#include "BackgroundMesh.h"
#if defined(HAVE_BFGS)
#include "stdafx.h"
@@ -588,6 +589,47 @@ void GFace::getMeanPlaneData(double plan[3][3]) const
plan[i][j] = meanPlane.plan[i][j];
}
+
+void GFace::computeMeshSizeFieldAccuracy(double &avg,double &max_e, double &min_e,
+ int &nE, int &GS)
+{
+
+ std::set<MEdge, Less_Edge> es;
+ for(unsigned int i = 0; i < getNumMeshElements(); i++){
+ MElement *e = getMeshElement(i);
+ for(int j = 0; j < e->getNumEdges(); j++) es.insert(e->getEdge(j));
+ }
+
+ avg = 0.0;
+ min_e = 1.e22;
+ max_e = 0;
+ nE = es.size();
+ GS = 0;
+ double oneoversqr2 = 1. / sqrt(2.);
+ double sqr2 = sqrt(2.);
+ for (std::set<MEdge, Less_Edge>::const_iterator it = es.begin();
+ it!=es.end();++it){
+ double u1,v1,u2,v2;
+ MVertex *vert1 = it->getVertex(0);
+ vert1->getParameter(0, u1);
+ vert1->getParameter(1, v1);
+ MVertex *vert2 = it->getVertex(1);
+ vert2->getParameter(0, u2);
+ vert2->getParameter(1, v2);
+ double l1 = BGM_MeshSize(this,u1,v1,vert1->x(), vert1->y(),vert1->z());
+ double l2 = BGM_MeshSize(this,u2,v2,vert2->x(), vert2->y(),vert2->z());
+ double correctLC = 0.5*(l1+l2);
+ double lone = it->length()/correctLC;
+ if (lone > oneoversqr2 && lone < sqr2) GS++;
+ double add = lone >1 ? (1. / lone) - 1. : lone - 1.;
+ avg += lone >1 ? (1. / lone) - 1. : lone - 1.;
+ max_e = std::max(max_e, lone);
+ min_e = std::min(min_e, lone);
+ }
+ avg = 100*exp(1./nE*avg);
+
+}
+
double GFace::curvatureDiv(const SPoint2 ¶m) const
{
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 9aa758b0a1d67d9c2ca37b763199fe9f36523b91..bc2470cadbc5434ab1e3f30e4a3284ea25b606e2 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -244,6 +244,10 @@ class GFace : public GEntity
// for that face
void moveToValidRange(SPoint2 &pt) const;
+ //compute mesh statistics
+ void computeMeshSizeFieldAccuracy(double &avg,double &max_e, double &min_e,
+ int &nE, int &GS);
+
// compound
void setCompound(GFaceCompound *gfc) { compound = gfc; }
GFaceCompound *getCompound() const { return compound; }
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 6eb737e7fd84c37677da51cef8b241815370e86d..a4a57fc38520bfc82a8eaf741349901e741e4ff0 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -502,7 +502,7 @@ double discreteEdge::curvature(double par) const
Curvature& curvature = Curvature::getInstance();
if( !Curvature::valueAlreadyComputed() ) {
std::cout << "Need to compute discrete curvature (in discreteEdge)" << std::endl;
- Curvature::typeOfCurvature type = Curvature::RBF; //RUSIN; //RBF
+ Curvature::typeOfCurvature type = Curvature::RUSIN; //RUSIN; //RBF
curvature.computeCurvature(model(), type);
}
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 302e9b2e787d69cbded6d27a33dd88b0e0fd4f54..ede2ecc2826a964a3a5320ca955162be9a3001e3 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -385,6 +385,7 @@ static void Mesh1D(GModel *m)
static void PrintMesh2dStatistics(GModel *m)
{
+
FILE *statreport = 0;
if(CTX::instance()->createAppendMeshStatReport == 1)
statreport = fopen(CTX::instance()->meshStatReportFileName.c_str(), "w");
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 378c5deacac796deb5db04e815f7a46b3a99a7b1..e31cb2459782a6c175f2a684caba7d3a6d537dc7 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1092,29 +1092,21 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
Msg::Debug("Starting to add internal points");
// start mesh generation
if(!algoDelaunay2D(gf) && !onlyInitialMesh){
- // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) {
- // printf("coucou here !!!\n");
- // backgroundMesh::unset();
- // buildBackGroundMesh (gf);
- // }
+ // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) {
+ // printf("coucou here !!!\n");
+ // backgroundMesh::unset();
+ // buildBackGroundMesh (gf);
+ // }
refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true,
&recoverMapInv);
optimizeMeshBDS(gf, *m, 2);
refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, false,
&recoverMapInv);
optimizeMeshBDS(gf, *m, 2);
- // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) {
- // backgroundMesh::unset();
- // }
+ // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) {
+ // backgroundMesh::unset();
+ // }
}
- /*
- computeMeshSizeFieldAccuracy(gf, *m, gf->meshStatistics.efficiency_index,
- gf->meshStatistics.longest_edge_length,
- gf->meshStatistics.smallest_edge_length,
- gf->meshStatistics.nbEdge,
- gf->meshStatistics.nbGoodLength);
- */
- gf->meshStatistics.status = GFace::DONE;
// fill the small gmsh structures
BDS2GMSH(m, gf, recoverMap);
@@ -1156,6 +1148,15 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
if(CTX::instance()->mesh.remove4triangles)
removeFourTrianglesNodes(gf,false);
+ //Emi print efficiency index
+ gf->computeMeshSizeFieldAccuracy(gf->meshStatistics.efficiency_index,
+ gf->meshStatistics.longest_edge_length,
+ gf->meshStatistics.smallest_edge_length,
+ gf->meshStatistics.nbEdge,
+ gf->meshStatistics.nbGoodLength);
+ printf("----- Efficiency index is tau=%g\n", gf->meshStatistics.efficiency_index);
+ gf->meshStatistics.status = GFace::DONE;
+
// delete the mesh
delete m;
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index aae39273833383baad3b4a5fb18bb33cf89eca3a..f0d41df92dc6e2aac060fea66892f0acc757c5fb 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -39,13 +39,12 @@ inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
double SCALINGU, double SCALINGV)
{
// FIXME SUPER HACK
-
- // if (CTX::instance()->mesh.recombineAll || f->meshAttributes.recombine || 1){
- // double quadAngle = backgroundMesh::current()->getAngle (0.5 * (p1->u + p2->u) * SCALINGU, 0.5 * (p1->v + p2->v) * SCALINGV,0);
- // const double a [2] = {p1->u,p1->v};
- // const double b [2] = {p2->u,p2->v};
- // return lengthInfniteNorm(a,b, quadAngle);
- // }
+ // if (CTX::instance()->mesh.recombineAll || f->meshAttributes.recombine || 1){
+ // double quadAngle = backgroundMesh::current()->getAngle (0.5 * (p1->u + p2->u) * SCALINGU, 0.5 * (p1->v + p2->v) * SCALINGV,0);
+ // const double a [2] = {p1->u,p1->v};
+ // const double b [2] = {p2->u,p2->v};
+ // return lengthInfniteNorm(a,b, quadAngle);
+ // }
GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
0.5 * (p1->v + p2->v) * SCALINGV));
@@ -53,7 +52,6 @@ inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
if (!GP.succeeded())
return computeEdgeLinearLength(p1,p2);
-
const double dx1 = p1->X - GP.x();
const double dy1 = p1->Y - GP.y();
const double dz1 = p1->Z - GP.z();
@@ -181,7 +179,7 @@ void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
double &max_e, double &min_e, int &nE, int &GS)
{
std::list<BDS_Edge*>::iterator it = m.edges.begin();
- avg = 0;
+ avg = 0.0;
min_e = 1.e22;
max_e = 0;
nE = 0;
@@ -199,6 +197,7 @@ void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
}
++it;
}
+ avg = 100*exp(1./nE*avg);
}
// SWAP TESTS i.e. tell if swap should be done
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index db458fe917a77b7b1b1160528555c5a19f426465..a2e88efa35f1a1877a0fbb69e95fd87a0eef7845 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -884,7 +884,7 @@ void bowyerWatson(GFace *gf, int MAXPNT)
*/
double lengthInfniteNorm(const double p[2], const double q[2],
- const double quadAngle){
+ const double quadAngle){
double xp = p[0] * cos(quadAngle) + p[1] * sin(quadAngle);
double yp = -p[0] * sin(quadAngle) + p[1] * cos(quadAngle);
double xq = q[0] * cos(quadAngle) + q[1] * sin(quadAngle);
@@ -1437,7 +1437,6 @@ void buildBackGroundMesh (GFace *gf)
gf->triangles = TR;
// gf->quadrangles = QR;
}
- // printf("end build bak mesh\n");
}
diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h
index cb64fb86c6f8039930e5aff619e465b872446c07..ecf7864bb0a58d01c9ea8d53b439ed54d7a587a8 100644
--- a/Solver/convexCombinationTerm.h
+++ b/Solver/convexCombinationTerm.h
@@ -58,11 +58,12 @@ class convexCombinationTerm : public femTerm<double> {
tan = sin(angle*0.5)/cos(angle*0.5);
}
double length = vj->distance(vk);
- m(j, k) = -tan/length;
- //m(j, k) = -1;
+ m(j, k) = -tan/length; // mean value convex map
+ //m(j, k) = -1; //convex map of Tutte
}
+ //m(j,j) = (nbSF - 1);// convex map of Tutte
for (int k = 0; k < nbSF; k++){
- if (k != j) diag += (-m(j,k));
+ if (k != j) diag += (-m(j,k));
}
m(j, j) = diag;
}
diff --git a/benchmarks/curvature/Bifurcation/BifurcationRemeshCurvature.geo b/benchmarks/curvature/Bifurcation/BifurcationRemeshCurvature.geo
index 7943fe5eb7d8b4fd86a357a922c6731b61f04f44..67e5dba1eb371c9fec4739998635b4265147c09a 100644
--- a/benchmarks/curvature/Bifurcation/BifurcationRemeshCurvature.geo
+++ b/benchmarks/curvature/Bifurcation/BifurcationRemeshCurvature.geo
@@ -1,9 +1,13 @@
Mesh.Algorithm=7; //1=MeshAdapt, 5=Delaunay, 6=Frontal, 7=BAMG
-Mesh.CharacteristicLengthFactor=0.35;
+
+Mesh.RemeshParametrization=7; //0=harmonic_circle, 1=conformal_spectral, 2=rbf, 3=harmonic_plane, 4=convex_circle, 5=convex_plane, 6=harmonic square, 7=conformal_fe
+Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split metis
+
+Mesh.CharacteristicLengthFactor=0.25;
Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
Mesh.CharacteristicLengthMin = 0.05; //-clmin 0.05
-Mesh.CharacteristicLengthMax = 3.00; //-clmax 4.0
-Mesh.LcIntegrationPrecision=1.e-5; //-epslc1d
+Mesh.CharacteristicLengthMax = 6.00; //-clmax 4.0
+Mesh.LcIntegrationPrecision=1.e-3; //-epslc1d
Mesh.MinimumCirclePoints=30; //default=7
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthExtendFromBoundary=0;
@@ -13,8 +17,6 @@ Merge "BifurcationInput.msh";
CreateTopology;
Compound Surface(20) = {1};
-Mesh.RemeshParametrization=1; //(0) harmonic (1) conformal
-Mesh.RemeshAlgorithm=0; //(0) nosplit (1) automatic (2) split only with metis ///Default: 1
Physical Surface("Wall") = {20};