diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 416e6bff5966845b035a52e0b09e3263b9d6b4ee..2863b339c3ce0449c3041f0abb057f9abafbfe00 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -85,6 +85,7 @@ void PrintUsage(const char *name)
Msg::Direct(" -rand float Set random perturbation factor");
Msg::Direct(" -bgm file Load background mesh from file");
Msg::Direct(" -check Perform various consistency checks on mesh");
+ Msg::Direct(" -mpass int Do several passes on the mesh for complex backround fields");
#if defined(HAVE_FLTK)
Msg::Direct("Post-processing options:");
Msg::Direct(" -link int Select link mode between views (0, 1, 2, 3, 4)");
@@ -438,6 +439,16 @@ void GetOptions(int argc, char *argv[])
else
Msg::Fatal("Missing number");
}
+ else if(!strcmp(argv[i] + 1, "mpass")) {
+ i++;
+ if(argv[i]) {
+ CTX::instance()->mesh.multiplePasses = atoi(argv[i++]);
+ if(CTX::instance()->mesh.multiplePasses <= 0)
+ Msg::Fatal("Number of Mesh Passes must be > 0");
+ }
+ else
+ Msg::Fatal("Missing number");
+ }
else if(!strcmp(argv[i] + 1, "edgelmin")) {
i++;
if(argv[i]) {
diff --git a/Common/Context.h b/Common/Context.h
index 466f44274c56d86067c16e466a5d3e22e56546f1..5dbe2af47c8414ce18e4ab3b6777b6822543e771 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -37,6 +37,7 @@ struct contextMeshOptions {
int saveElementTagType;
int switchElementTags;
int highOrderNoMetric;
+ int multiplePasses;
};
struct contextGeometryOptions {
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 25a2c2ac6435fb9b419ede1351c2384812d123d0..d1d4947561dd37a0578f31b279953d0da705c65c 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1152,6 +1152,8 @@ StringXNumber MeshOptions_Number[] = {
"Version of the MSH file format to use" },
{ F|O, "MshFilePartitioned" , opt_mesh_msh_file_partitioned , 0. ,
"Split MSH file by mesh partition" },
+ { F|O, "MultiplePassesMeshes" , opt_mesh_multiple_passes , 0. ,
+ "Do a first simple mesh and use it for complex background meshes (curvatures...)" },
{ F|O, "PartitionHexWeight" , opt_mesh_partition_hex_weight , 1 ,
"Weight of hexahedral element for METIS load balancing" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 931c14bf60cf88242b54ac795567063e59c931a5..9cd4ce697ce5b7d5b29acf5b09f705f663374a07 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5839,6 +5839,13 @@ double opt_mesh_second_order_experimental(OPT_ARGS_NUM)
return CTX::instance()->mesh.secondOrderExperimental;
}
+double opt_mesh_multiple_passes(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET)
+ CTX::instance()->mesh.multiplePasses = (int)val;
+ return CTX::instance()->mesh.multiplePasses;
+}
+
double opt_mesh_second_order_linear(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 1f13d6f411fe6365be60eb35cd7af16bb2fef2c1..a43ed93917bf50bb9c7db78b85a0ff3454655849 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -539,6 +539,7 @@ double opt_mesh_mesh_only_visible(OPT_ARGS_NUM);
double opt_mesh_min_circ_points(OPT_ARGS_NUM);
double opt_mesh_allow_swap_edge_angle(OPT_ARGS_NUM);
double opt_mesh_min_curv_points(OPT_ARGS_NUM);
+double opt_mesh_multiple_passes(OPT_ARGS_NUM);
double opt_mesh_order(OPT_ARGS_NUM);
double opt_mesh_smooth_internal_edges(OPT_ARGS_NUM);
double opt_mesh_second_order_experimental(OPT_ARGS_NUM);
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index eb2544dfdb600ea971c808b73b7b24f89cd2f9cf..b16f0f125f32c129240aa9bda499167fbbc21018 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -16,6 +16,7 @@
#include "MElementOctree.h"
#include "MLine.h"
#include "MTriangle.h"
+#include "MQuadrangle.h"
#include "MVertex.h"
#if defined(HAVE_SOLVER)
@@ -410,6 +411,10 @@ void backgroundMesh::unset()
backgroundMesh::backgroundMesh(GFace *_gf)
{
+ // create a bunch of triangles on the parametric space
+ // those triangles are local to the backgroundMesh so that
+ // they do not depend on the actual mesh that can be deleted
+
for (unsigned int i = 0; i < _gf->triangles.size(); i++){
MTriangle *e = _gf->triangles[i];
MVertex *news[3];
@@ -429,7 +434,11 @@ backgroundMesh::backgroundMesh(GFace *_gf)
}
_triangles.push_back(new MTriangle(news[0],news[1],news[2]));
}
+
+ // build a search structure
_octree = new MElementOctree(_triangles);
+
+ // compute the mesh sizes at nodes
if (CTX::instance()->mesh.lcFromPoints)
propagate1dMesh(_gf);
else {
@@ -438,7 +447,12 @@ backgroundMesh::backgroundMesh(GFace *_gf)
_sizes[itv2->first] = MAX_LC;
}
}
+ // ensure that other criteria are fullfilled
updateSizes(_gf);
+
+ // compute optimal mesh orientations
+ propagatecrossField(_gf);
+
_3Dto2D.clear();
_2Dto3D.clear();
}
@@ -450,13 +464,9 @@ backgroundMesh::~backgroundMesh()
delete _octree;
}
-void backgroundMesh::propagate1dMesh(GFace *_gf)
+static void propagateValuesOnFace (GFace *_gf,
+ std::map<MVertex*,double> &dirichlet)
{
-#if defined(HAVE_SOLVER)
- std::list<GEdge*> e = _gf->edges();
- std::list<GEdge*>::const_iterator it = e.begin();
- std::map<MVertex*,double> sizes;
-
linearSystem<double> *_lsys = 0;
#if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
_lsys = new linearSystemPETSc<double>;
@@ -469,8 +479,50 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
#else
_lsys = new linearSystemFull<double>;
#endif
-
+
dofManager<double> myAssembler(_lsys);
+
+ // fix boundary conditions
+ std::map<MVertex*, double>::iterator itv = dirichlet.begin();
+ for ( ; itv != dirichlet.end(); ++itv){
+ myAssembler.fixVertex(itv->first, 0, 1, itv->second);
+ }
+
+
+ // Number vertices
+ std::set<MVertex*> vs;
+ for (unsigned int k = 0; k < _gf->triangles.size(); k++)
+ for (int j=0;j<3;j++)vs.insert(_gf->triangles[k]->getVertex(j));
+ for (unsigned int k = 0; k < _gf->quadrangles.size(); k++)
+ for (int j=0;j<4;j++)vs.insert(_gf->quadrangles[k]->getVertex(j));
+ for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it)
+ myAssembler.numberVertex(*it, 0, 1);
+
+ // Assemble
+ simpleFunction<double> ONE(1.0);
+ laplaceTerm l(0, 1, &ONE);
+ for (unsigned int k = 0; k < _gf->triangles.size(); k++){
+ MTriangle *t = _gf->triangles[k];
+ SElement se(t);
+ l.addToMatrix(myAssembler, &se);
+ }
+
+ // Solve
+ _lsys->systemSolve();
+
+ // save solution
+ for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
+ myAssembler.getDofValue(*it, 0, 1, dirichlet[*it]);
+ }
+
+ delete _lsys;
+}
+
+void backgroundMesh::propagate1dMesh(GFace *_gf)
+{
+ std::list<GEdge*> e = _gf->edges();
+ std::list<GEdge*>::const_iterator it = e.begin();
+ std::map<MVertex*,double> sizes;
for( ; it != e.end(); ++it ){
if (!(*it)->isSeam(_gf)){
@@ -484,48 +536,93 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
MVertex *v = (*it)->lines[i]->getVertex(k);
std::map<MVertex*, double>::iterator itv = sizes.find(v);
if (itv == sizes.end())
- sizes[v] = d;
+ sizes[v] = log(d);
else
- itv->second = 0.5 * (itv->second + d);
+ itv->second = 0.5 * (itv->second + log(d));
}
}
}
}
- std::map<MVertex*, double>::iterator itv = sizes.begin();
- for ( ; itv != sizes.end(); ++itv){
- myAssembler.fixVertex(itv->first, 0, 1, itv->second);
+
+ propagateValuesOnFace(_gf,sizes);
+
+ std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
+ for ( ; itv2 != _2Dto3D.end(); ++itv2){
+ MVertex *v_2D = itv2->first;
+ MVertex *v_3D = itv2->second;
+ _sizes[v_2D] = exp(sizes[v_3D]);
}
+}
- simpleFunction<double> ONE(1.0);
- laplaceTerm l(0, 1, &ONE);
+// C R O S S F I E L D S
- for (unsigned int k = 0; k < _gf->triangles.size(); k++){
- MTriangle *t = _gf->triangles[k];
- myAssembler.numberVertex(t->getVertex(0), 0, 1);
- myAssembler.numberVertex(t->getVertex(1), 0, 1);
- myAssembler.numberVertex(t->getVertex(2), 0, 1);
+crossField2d :: crossField2d (MVertex* v, GEdge* ge){
+ double p;
+ bool success = reparamMeshVertexOnEdge(v, ge, p);
+ if (!success){
+ Msg::Warning("cannot reparametrize a point in crossField");
+ _angle = 0;
+ return;
}
+ SVector3 t = ge->firstDer (p);
+ t.normalize();
+ _angle = atan2 (t.y(),t.x());
+ crossField2d::normalizeAngle (_angle);
+}
+
+
+
+
+void backgroundMesh::propagatecrossField(GFace *_gf)
+{
+ std::map<MVertex*,double> _cosines4,_sines4;
+ std::list<GEdge*> e = _gf->edges();
+ std::list<GEdge*>::const_iterator it = e.begin();
- for (unsigned int k = 0; k < _gf->triangles.size(); k++){
- MTriangle *t = _gf->triangles[k];
- SElement se(t);
- l.addToMatrix(myAssembler, &se);
+ for( ; it != e.end(); ++it ){
+ if (!(*it)->isSeam(_gf)){
+ for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
+ MVertex *v[2];
+ v[0] = (*it)->lines[i]->getVertex(0);
+ v[1] = (*it)->lines[i]->getVertex(1);
+ SPoint2 p1,p2;
+ bool success = reparamMeshEdgeOnFace(v[0],v[1],_gf,p1,p2);
+ // double angle = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
+ double angle = atan2 ( v[0]->y()-v[1]->y() , v[0]->x()- v[1]->x() );
+ crossField2d::normalizeAngle (angle);
+ for (int i=0;i<2;i++){
+ std::map<MVertex*,double>::iterator itc = _cosines4.find(v[i]);
+ std::map<MVertex*,double>::iterator its = _sines4.find(v[i]);
+ if (itc != _cosines4.end()){
+ itc->second = 0.5*(itc->second + cos(4*angle));
+ its->second = 0.5*(its->second + sin(4*angle));
+ }
+ else {
+ _cosines4[v[i]] = cos(4*angle);
+ _sines4[v[i]] = sin(4*angle);
+ }
+ }
+ }
+ }
}
- _lsys->systemSolve();
+
+ propagateValuesOnFace(_gf,_cosines4);
+ propagateValuesOnFace(_gf,_sines4);
std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
for ( ; itv2 != _2Dto3D.end(); ++itv2){
MVertex *v_2D = itv2->first;
MVertex *v_3D = itv2->second;
- myAssembler.getDofValue(v_3D, 0, 1, _sizes[v_2D]);
+ double angle = atan2(_sines4[v_3D],_cosines4[v_3D]) / 4.0;
+ crossField2d::normalizeAngle (angle);
+ _angles[v_2D] = angle;
}
- delete _lsys;
-#endif
+
}
+
void backgroundMesh::updateSizes(GFace *_gf)
{
-#if defined(HAVE_SOLVER)
std::map<MVertex*,double>::iterator itv = _sizes.begin();
for ( ; itv != _sizes.end(); ++itv){
SPoint2 p;
@@ -548,80 +645,6 @@ void backgroundMesh::updateSizes(GFace *_gf)
itv->second = std::max(itv->second, CTX::instance()->mesh.lcMin);
itv->second = std::min(itv->second, CTX::instance()->mesh.lcMax);
}
- // return;
- // now do some diffusion
-
- std::list<GEdge*> e = _gf->edges();
- std::list<GEdge*>::const_iterator it = e.begin();
-
- linearSystem<double> *_lsys = 0;
-#if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
- _lsys = new linearSystemPETSc<double>;
-#elif defined(HAVE_GMM) && !defined(HAVE_TAUCS)
- linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
- _lsysb->setGmres(1);
- _lsys = _lsysb;
-#elif defined(HAVE_TAUCS)
- _lsys = new linearSystemCSRTaucs<double>;
-#else
- _lsys = new linearSystemFull<double>;
-#endif
-
- dofManager<double> myAssembler(_lsys);
-
- // M (U^{t+1} - U^{t})/DT + K U^{t+1} = 0
-
- for( ; it != e.end(); ++it ){
- if (!(*it)->isSeam(_gf)){
- for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
- MVertex *v1 = (*it)->lines[i]->getVertex(0);
- MVertex *v2 = (*it)->lines[i]->getVertex(1);
- MVertex *v1_2D = _3Dto2D[ v1 ];
- MVertex *v2_2D = _3Dto2D[ v2 ];
-
- myAssembler.fixVertex(v1, 0, 1, _sizes[v1_2D]);
- myAssembler.fixVertex(v2, 0, 1, _sizes[v2_2D]);
- }
- }
- }
-
- double DT = 0.01;
- simpleFunction<double> ONEOVERDT(1.0/DT);
- simpleFunction<double> ONE(1.0);
- helmholtzTerm<double> diffusionTerm(0,1, 1, &ONE,&ONEOVERDT);
- helmholtzTerm<double> massTerm(0,1, 1, 0,&ONEOVERDT);
-
- for (unsigned int k= 0; k< _gf->triangles.size(); k++){
- MTriangle *t = _gf->triangles[k];
- myAssembler.numberVertex(t->getVertex(0), 0, 1);
- myAssembler.numberVertex(t->getVertex(1), 0, 1);
- myAssembler.numberVertex(t->getVertex(2), 0, 1);
- }
- fullMatrix<double> mass_matrix(3,3);
- for (unsigned int k= 0; k< _gf->getNumMeshElements(); k++){
- MElement *e = _gf->getMeshElement(k);
- SElement se(e);
- diffusionTerm.addToMatrix(myAssembler, &se);
- massTerm.elementMatrix(&se,mass_matrix);
- for (int i = 0; i < e->getNumVertices();i++){
- double TERM = 0.0;
- for (int j = 0; j < e->getNumVertices();j++){
- MVertex *v_2D = _3Dto2D[ e->getVertex(j) ];
- TERM+=_sizes[v_2D] * mass_matrix(i,j);
- }
- myAssembler.assemble(e->getVertex(i),0,1,TERM);
- }
- }
- _lsys->systemSolve();
-
- std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
- for ( ; itv2 != _2Dto3D.end(); ++itv2){
- MVertex *v_2D = itv2->first;
- MVertex *v_3D = itv2->second;
- myAssembler.getDofValue(v_3D, 0, 1, _sizes[v_2D]);
- }
- delete _lsys;
-#endif
}
double backgroundMesh::operator() (double u, double v, double w) const
@@ -641,7 +664,35 @@ double backgroundMesh::operator() (double u, double v, double w) const
return itv1->second * (1-uv2[0]-uv2[1]) + itv2->second * uv2[0] + itv3->second * uv2[1];
}
-void backgroundMesh::print (const std::string &filename, GFace *gf) const
+double backgroundMesh::getAngle (double u, double v, double w) const
+{
+ double uv[3] = {u, v, w};
+ double uv2[3];
+ // return 1.0;
+ MElement *e = _octree->find(u, v, w);
+ if (!e) {
+ Msg::Error("cannot find %g %g", u, v);
+ return 0.0;
+ }
+ e->xyz2uvw(uv, uv2);
+ std::map<MVertex*,double>::const_iterator itv1 = _angles.find(e->getVertex(0));
+ std::map<MVertex*,double>::const_iterator itv2 = _angles.find(e->getVertex(1));
+ std::map<MVertex*,double>::const_iterator itv3 = _angles.find(e->getVertex(2));
+
+ double cos4 = cos (4*itv1->second) * (1-uv2[0]-uv2[1]) +
+ cos (4*itv2->second) * uv2[0] +
+ cos (4*itv3->second) * uv2[1] ;
+ double sin4 = sin (4*itv1->second) * (1-uv2[0]-uv2[1]) +
+ cos (4*itv2->second) * uv2[0] +
+ cos (4*itv3->second) * uv2[1] ;
+ double angle = atan2(sin4,cos4)/4.0;
+ crossField2d::normalizeAngle (angle);
+
+ return angle;
+}
+
+
+void backgroundMesh::print (const std::string &filename, GFace *gf, const std::map<MVertex*,double> &_whatToPrint) const
{
FILE *f = fopen (filename.c_str(),"w");
fprintf(f,"View \"Background Mesh\"{\n");
@@ -649,9 +700,9 @@ void backgroundMesh::print (const std::string &filename, GFace *gf) const
MVertex *v1 = _triangles[i]->getVertex(0);
MVertex *v2 = _triangles[i]->getVertex(1);
MVertex *v3 = _triangles[i]->getVertex(2);
- std::map<MVertex*,double>::const_iterator itv1 = _sizes.find(v1);
- std::map<MVertex*,double>::const_iterator itv2 = _sizes.find(v2);
- std::map<MVertex*,double>::const_iterator itv3 = _sizes.find(v3);
+ std::map<MVertex*,double>::const_iterator itv1 = _whatToPrint.find(v1);
+ std::map<MVertex*,double>::const_iterator itv2 = _whatToPrint.find(v2);
+ std::map<MVertex*,double>::const_iterator itv3 = _whatToPrint.find(v3);
if (!gf){
fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
v1->x(),v1->y(),v1->z(),
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index 8477fbaf31ae9a264a8af2cadb408cd7cbbbd1df..f98db7599d65ba80e70f18539a2f9123206b155b 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -18,6 +18,21 @@ class MElement;
class MVertex;
class GEntity;
+struct crossField2d
+{
+ double _angle;
+ static void normalizeAngle (double &angle) {
+ if (angle < 0)
+ while ( angle < 0 ) angle += (M_PI * .5);
+ else if (angle >= M_PI * .5)
+ while ( angle >= M_PI * .5 ) angle -= (M_PI * .5);
+ }
+ crossField2d (MVertex*, GEdge*);
+ crossField2d (double a) : _angle(a){}
+ crossField2d & operator += ( const crossField2d & );
+};
+
+
class backgroundMesh : public simpleFunction<double>
{
MElementOctree *_octree;
@@ -27,6 +42,7 @@ class backgroundMesh : public simpleFunction<double>
std::map<MVertex*,MVertex*> _3Dto2D;
std::map<MVertex*,MVertex*> _2Dto3D;
std::map<MVertex*,double> _distance;
+ std::map<MVertex*,double> _angles;
static backgroundMesh * _current;
backgroundMesh(GFace *);
~backgroundMesh();
@@ -35,9 +51,17 @@ class backgroundMesh : public simpleFunction<double>
static void unset();
static backgroundMesh *current () { return _current; }
void propagate1dMesh(GFace *);
+ void propagatecrossField(GFace *);
void updateSizes(GFace *);
double operator () (double u, double v, double w) const;
- void print (const std::string &filename, GFace *gf) const;
+ double getAngle(double u, double v, double w) const ;
+ void print (const std::string &filename, GFace *gf, const std::map<MVertex*,double>&) const;
+ void print (const std::string &filename, GFace *gf, int choice = 0) const {
+ switch(choice) {
+ case 0 : print(filename,gf,_sizes); return;
+ default : print(filename,gf,_angles); return;
+ }
+ }
};
double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double Z);
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index ad5d796fcf0e74e171a29b21bf07b5951ae62695..cc79e9695cc2899683e1709b50c21f034fb0813c 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -22,6 +22,7 @@
#include "GmshMessage.h"
#include "Numeric.h"
#include "mathEvaluator.h"
+#include "BackgroundMesh.h"
#if defined(HAVE_POST)
#include "OctreePost.h"
@@ -1710,10 +1711,16 @@ public:
// (dist/hwall)*(ratio-1) + 1 = ratio^{m}
// lc = dist*(ratio-1) + hwall
const double ll1 = dist*(ratio-1) + hwall_n;
- const double lc_n = std::min(ll1,hfar);
+ double lc_n = std::min(ll1,hfar);
const double ll2 = dist*(ratio-1) + hwall_t;
double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2,hfar));
+ if (backgroundMesh::current()){
+ const double lcBG = backgroundMesh::current()->operator() (x,y,z);
+ lc_n = std::min(lc_n, lcBG);
+ lc_t = std::min(lc_t, lcBG);
+ }
+
SVector3 t1,t2,t3;
double L1,L2,L3;
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 293861cb3ff911adffbcaff9fd5328d87c7e6e82..abbc9828ee7a606b3a3eb3475f5b30e954c0708d 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -448,16 +448,17 @@ static void Mesh2D(GModel *m)
int nIter = 0;
while(1){
- meshGFace mesher;
int nbPending = 0;
for(std::set<GFace*>::iterator it = f.begin(); it != f.end(); ++it){
if ((*it)->meshStatistics.status == GFace::PENDING){
+ meshGFace mesher (true, CTX::instance()->mesh.multiplePasses);
mesher(*it);
nbPending++;
}
}
for(std::set<GFace*>::iterator it = cf.begin(); it != cf.end(); ++it){
if ((*it)->meshStatistics.status == GFace::PENDING){
+ meshGFace mesher (true, CTX::instance()->mesh.multiplePasses);
mesher(*it);
nbPending++;
}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 59920098e7a39d33f764555230808e4170b74eb6..240f1173dfe8dd782d4e8ebd9fddd84f0941a5a3 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -29,32 +29,7 @@
static bool mappingIsInvertible(MTetrahedron *e)
{
if (e->getPolynomialOrder() == 1) return true;
-
- double mat[3][3];
- e->getPrimaryJacobian(0., 0., 0., mat);
- double det0 = det3x3(mat);
-
- fullMatrix<double> df;
- {
- const fullMatrix<double> &alldf =
- e->getGradShapeFunctionsAtIntegrationPoints(e->getPolynomialOrder());
- for (int i = 0; i < alldf.size2()/3; i++){
- df.setAsProxy(alldf, 3*i, 3);
- e->getJacobian(df, mat);
- if (det0 * det3x3(mat) <= 0.) return false;
- }
- }
- {
- const fullMatrix<double> &points = e->getFunctionSpace()->points;
- const fullMatrix<double> &alldf =
- e->getGradShapeFunctionsAtNodes(e->getPolynomialOrder());
- for (int i = 0; i < alldf.size2() / 3; i++){
- df.setAsProxy(alldf, 3*i, 3);
- e->getJacobian(df, mat);
- if (det0 * det3x3(mat) <= 0.) return false;
- }
- }
- return true;
+ if (e->getPolynomialOrder() == 1) return e->distoShapeMeasure() > 0;
}
// The aim here is to build a polynomial representation that consist
@@ -267,8 +242,8 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
if(relax < 1.e-2)
Msg::Warning
("Failed to compute equidistant parameters (relax = %g, value = %g) "
- "for edge %d-%d parametrized with %g %g on GEdge %d",
- relax, US[1], v0->getNum(), v1->getNum(),u0,u1,ge->tag());
+ "for edge %d-%d parametrized with %g %g on GEdge %d linear %d",
+ relax, US[1], v0->getNum(), v1->getNum(),u0,u1,ge->tag(), linear);
}
}
for(int j = 0; j < nPts; j++){
@@ -277,8 +252,9 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
double uc = (1. - t) * u0 + t * u1; // can be wrong, that's ok
MVertex *v;
if(linear || !reparamOK || uc < std::min(u0,u1) || uc > std::max(u0,u1)){
- Msg::Warning("We don't have a valid parameter on curve %d-%d",
- v0->getNum(), v1->getNum());
+ if (!linear)
+ Msg::Warning("We don't have a valid parameter on curve %d-%d",
+ v0->getNum(), v1->getNum());
// we don't have a (valid) parameter on the curve
SPoint3 pc = edge.interpolate(t);
v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
@@ -1256,10 +1232,10 @@ static void checkHighOrderTetrahedron(const char* cc, GModel *m,
}
}
if (minJGlob < 0)
- Msg::Info("%s : Worst Tetrahedron Smoothness %g Gamma %g NbFair = %d",
- cc, minJGlob, minGGlob, nbfair);
- else
- Msg::Warning("%s : Worst Tetrahedron Smoothness %g (%d negative jacobians) "
+ Msg::Warning("%s : Worst Tetrahedron Smoothness %g Gamma %g NbFair = %d NbBad = %d",
+ cc, minJGlob, minGGlob, nbfair, bad.size());
+ else
+ Msg::Info("%s : Worst Tetrahedron Smoothness %g (%d negative jacobians) "
"Worst Gamma %g Avg Smoothness %g", cc, minJGlob, bad.size(),
minGGlob, avg / (count ? count : 1));
}
@@ -1407,7 +1383,8 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
if (!linear){
hot.ensureMinimumDistorsion(0.1);
- checkHighOrderTriangles("final mesh", m, bad, worst);
+ checkHighOrderTriangles("Final mesh", m, bad, worst);
+ checkHighOrderTetrahedron("Final mesh", m, bad, worst);
}
// m->writeMSH("CORRECTED.msh");
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index 42f21fcbf2f8a44838d86468f41061089b8356d4..65a335a661fc889c85be69eaf69b8f152a46987b 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -402,7 +402,7 @@ void highOrderSmoother::optimize(GFace * gf,
// }
// then try to swap for better configurations
- smooth(gf, !CTX::instance()->mesh.highOrderNoMetric);
+ //smooth(gf, !CTX::instance()->mesh.highOrderNoMetric);
// for (int i=0;i<100;i++){
@@ -411,7 +411,7 @@ void highOrderSmoother::optimize(GFace * gf,
// printf("%d swaps\n",nbSwap);
// }
- // smooth(gf,true);
+ smooth(gf,true);
// smooth(gf,true);
// smooth(gf,true);
// smooth(gf,true);
@@ -838,10 +838,8 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
{
#if defined(HAVE_TAUCS)
linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
- printf("coucou\n");
#elif defined(HAVE_PETSC)
linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
- printf("coucou PETSC\n");
#else
linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
lsys->setNoisy(1);
@@ -854,9 +852,9 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
std::vector<MElement*> layer, v;
double minD;
- getDistordedElements(all, 0.85, v, minD);
+ getDistordedElements(all, 0.15, v, minD);
- Msg::Info("%d elements / %d distorted min Disto = %g\n",
+ Msg::Info("%d elements / %d distorted min Disto = %g",
all.size(), v.size(), minD);
if (!v.size());
diff --git a/Mesh/highOrderTools.cpp b/Mesh/highOrderTools.cpp
index c20e7ba8fdcc2aab7f201153527143562955a45e..ea884c357ece57e0a21819ff511f9ed732b1230c 100644
--- a/Mesh/highOrderTools.cpp
+++ b/Mesh/highOrderTools.cpp
@@ -481,6 +481,55 @@ void highOrderTools::computeStraightSidedPositions ()
}
}
+ for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
+ for (int i=0;i<(*it)->mesh_vertices.size();i++){
+ MVertex *v = (*it)->mesh_vertices[i];
+ _targetLocation[v] = SVector3(v->x(),v->y(),v->z());
+ }
+ for (int i=0;i<(*it)->tetrahedra.size();i++){
+ _dim = 3;
+ MTetrahedron *t = (*it)->tetrahedra[i];
+ MTetrahedron t_1 ((*it)->tetrahedra[i]->getVertex(0),
+ (*it)->tetrahedra[i]->getVertex(1),
+ (*it)->tetrahedra[i]->getVertex(2),
+ (*it)->tetrahedra[i]->getVertex(3));
+ const polynomialBasis* fs = t->getFunctionSpace();
+ for (int j=0;j<t->getNumVertices();j++){
+ if (t->getVertex(j)->onWhat() == *it){
+ double t1 = fs->points(j, 0);
+ double t2 = fs->points(j, 1);
+ double t3 = fs->points(j, 2);
+ SPoint3 pc; t_1.pnt(t1, t2, t3,pc);
+ _straightSidedLocation [t->getVertex(j)] =
+ SVector3(pc.x(),pc.y(),pc.z());
+ }
+ }
+ }
+ for (int i=0;i<(*it)->hexahedra.size();i++){
+ _dim = 3;
+ MHexahedron *h = (*it)->hexahedra[i];
+ MHexahedron h_1 ((*it)->hexahedra[i]->getVertex(0),
+ (*it)->hexahedra[i]->getVertex(1),
+ (*it)->hexahedra[i]->getVertex(2),
+ (*it)->hexahedra[i]->getVertex(3),
+ (*it)->hexahedra[i]->getVertex(4),
+ (*it)->hexahedra[i]->getVertex(5),
+ (*it)->hexahedra[i]->getVertex(6),
+ (*it)->hexahedra[i]->getVertex(7));
+ const polynomialBasis* fs = h->getFunctionSpace();
+ for (int j=0;j<h->getNumVertices();j++){
+ if (h->getVertex(j)->onWhat() == *it){
+ double t1 = fs->points(j, 0);
+ double t2 = fs->points(j, 1);
+ double t3 = fs->points(j, 2);
+ SPoint3 pc; h_1.pnt(t1, t2, t3,pc);
+ _straightSidedLocation [h->getVertex(j)] =
+ SVector3(pc.x(),pc.y(),pc.z());
+ }
+ }
+ }
+ }
+
Msg::Info("highOrderTools has been set up : %d nodes are considered",_straightSidedLocation.size());
}
@@ -621,10 +670,10 @@ void highOrderTools::ensureMinimumDistorsion (std::vector<MElement*> &all,
double minD;
std::vector<MElement*> disto;
getDistordedElements(all, threshold, disto, minD);
- if (num == disto.size())break;
+ // if (num == disto.size())break;
if (!disto.size())break;
num = disto.size();
- Msg::Info("fixing %d bad curved elements (worst disto %g)",disto.size(),minD);
+ Msg::Info("Fixing %d bad curved elements (worst disto %g)",disto.size(),minD);
for (int i=0;i<disto.size();i++){
ensureMinimumDistorsion(disto[i],threshold);
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 2fff56356f959c9d82d2d6b1dc4dd1f619eab0fe..828d00a090ba3fe71191203eb04344daf5662393 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1576,19 +1576,22 @@ void meshGFace::operator() (GFace *gf)
Msg::Debug("Type %d %d triangles generated, %d internal vertices",
gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size());
- // test : recompute the background mesh using a PDE
- /*
- if (backgroundMesh::current()){
+ // do the 2D mesh in several passes. For second and other passes,
+ // a background mesh is constructed with the previous mesh and
+ // nodal values of the metric are computed that take into account
+ // complex size fields that are tedious to evaluate on the fly
+ if (!twoPassesMesh)return;
+ twoPassesMesh--;
+ if (backgroundMesh::current()){
backgroundMesh::unset();
- }
- else{
- backgroundMesh::set(gf);
- char name[256];
- sprintf(name,"bgm-%d.pos",gf->tag());
- backgroundMesh::current()->print(name,gf);
- (*this)(gf);
- }
- */
+ }
+ backgroundMesh::set(gf);
+ char name[256];
+ sprintf(name,"bgm-%d.pos",gf->tag());
+ backgroundMesh::current()->print(name,gf);
+ sprintf(name,"cross-%d.pos",gf->tag());
+ backgroundMesh::current()->print(name,gf,1);
+ (*this)(gf);
}
bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges)
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index e2c870cb331ebaef53010147a0477648cce98f27..3f79ce359b3926e0b8321207d80b719eacb97906 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -17,15 +17,16 @@ class GFaceCompound;
// Create the mesh of the face
class meshGFace {
const bool repairSelfIntersecting1dMesh;
+ int twoPassesMesh;
public :
- meshGFace (bool r = true, bool s = false) : repairSelfIntersecting1dMesh(r){}
+ meshGFace (bool r = true, int t = 0) : repairSelfIntersecting1dMesh(r), twoPassesMesh(t){}
void operator () (GFace *);
};
// Destroy the mesh of the face
class deMeshGFace {
public :
- deMeshGFace (bool s = false){}
+ deMeshGFace (){}
void operator () (GFace *);
};
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 7478cde5149ca89521b5ea2fd94588b487fc8674..ef3e7c514c2efa1e292e2832551e3d0cf6fead38 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -558,7 +558,8 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
(coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
(coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
mid->X,mid->Y,mid->Z);
- mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
+ // mid->lc() = exp(0.5 * (log(e->p1->lc()) + log(e->p2->lc())));
+ mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
}
if(!m.split_edge(e, mid)) m.del_point(mid);
else nb_split++;
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index eb635f6c31373825cffeb21743042eb6ee35217f..d166913081869e9ab61ecbd61c9dd856f387e1bc 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -616,10 +616,17 @@ static void insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
// printf("the new point is %g %g\n",p.x(),p.y());
MVertex *v = new MFaceVertex(p.x(), p.y(), p.z(), gf, center[0], center[1]);
v->setIndex(Us.size());
- double lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getIndex()] +
- uv[0] * vSizes[ptin->tri()->getVertex(1)->getIndex()] +
- uv[1] * vSizes[ptin->tri()->getVertex(2)->getIndex()]);
- double lc = BGM_MeshSize(gf, center[0], center[1], p.x(), p.y(), p.z());
+ double lc1,lc;
+ if (backgroundMesh::current()){
+ lc1 = lc =
+ backgroundMesh::current()->operator()(center[0], center[1], 0.0);
+ }
+ else {
+ lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getIndex()] +
+ uv[0] * vSizes[ptin->tri()->getVertex(1)->getIndex()] +
+ uv[1] * vSizes[ptin->tri()->getVertex(2)->getIndex()]);
+ lc = BGM_MeshSize(gf, center[0], center[1], p.x(), p.y(), p.z());
+ }
//SMetric3 metr = BGM_MeshMetric(gf, center[0], center[1], p.x(), p.y(), p.z());
// vMetricsBGM.push_back(metr);
vSizesBGM.push_back(lc);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 3705a67305a9332c41ee51b94fa208eac07b4c90..911f522b669c0bfea356af8fe8ed2a35d891df68 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1764,7 +1764,7 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
for (int k=0;k<elist[0];k++){
int i1 = elist[1+3*k],i2 = elist[1+3*k+1],an=elist[1+3*k+2];
// FIXME !!
- if (an == 100000){// || an == 1000){
+ if (an == 100000 || an == 1000){
toProcess.push_back(std::make_pair(n2t[i1],n2t[i2]));
Msg::Debug("Extra edge found in blossom algorithm, optimization will be required");
}
diff --git a/Mesh/meshRefine.cpp b/Mesh/meshRefine.cpp
index 2ee8cca87607a0ebde0c3cb2099a7baef77bb009..1b4fa1e6f65bf0685710a5d0b0cd65eeb8f12219 100644
--- a/Mesh/meshRefine.cpp
+++ b/Mesh/meshRefine.cpp
@@ -377,11 +377,17 @@ void RefineMesh(GModel *m, bool linear, bool splitIntoQuads, bool splitIntoHexas
// mesh
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
Subdivide(*it);
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
Subdivide(*it, splitIntoQuads, splitIntoHexas, faceVertices);
+ if (splitIntoQuads){
+ recombineIntoQuads(*it,true,true);
+ }
+ }
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
Subdivide(*it, splitIntoHexas, faceVertices);
+
+
double t2 = Cpu();
Msg::StatusBar(2, true, "Done refining mesh (%g s)", t2 - t1);
}
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index c311c180b0fa43e82cf5d24a1141a016f88ddc16..18fe5f57aa5f5e87c7359fbe14ff4a97c12ae9b5 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -203,10 +203,10 @@ double mesh_functional_distorsion(MElement *t, double u, double v)
double sign = 1.0;
prosca(normal1, normal, &sign);
double det = norm3(normal) * (sign > 0 ? 1. : -1.) / nn;
- //double det = norm3(normal);
- // compute distorsion
- // double dist = std::min(1. / det, det);
+ // printf("%g %g : %g : %g n1 (%g,%g,%g)\n",u,v,sign,det, normal1[0], normal1[1], normal1[2]);
+ // printf("n (%g,%g,%g)\n", normal[0], normal[1], normal[2]);
+
return det;
}
@@ -305,17 +305,27 @@ double mesh_functional_distorsion_p2_exact(MTriangle *t)
double mesh_functional_distorsion_pN(MElement *t)
{
const bezierBasis *jac = t->getJacobianFuncSpace()->bezier;
- // jac->monomials.print("coucou");
fullVector<double>Ji(jac->points.size1());
+ // printf("%d points for bez \n",jac->points.size1());
for (int i=0;i<jac->points.size1();i++){
- const double u = jac->points(i,0);
- const double v = jac->points(i,1);
- Ji(i) = mesh_functional_distorsion(t,u,v);
+ double u = jac->points(i,0);
+ double v = jac->points(i,1);
+
+ // JF : bezier points are defined in the [0,1] x [0,1] quad
+ if (t->getType() == TYPE_QUA){
+ u = -1 + 2*u;
+ v = -1 + 2*v;
+ }
+
+ Ji(i) = mesh_functional_distorsion(t,u,v);
+ // printf("J(%g,%g) = %12.5E\n",u,v,Ji(i));
}
fullVector<double> Bi( jac->matrixLag2Bez.size1() );
jac->matrixLag2Bez.mult(Ji,Bi);
+ // jac->matrixLag2Bez.print("Lag2Bez");
+
return *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
}
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 5776a7109f5fae8e541ba7ca07b275ca8d2768d4..0ef87ad82679a4c6dbf1b8c368766eb14fdcd64d 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -108,6 +108,8 @@ static fullMatrix<double> generateExposantsQuad(int order)
}
}
+ // exposants.print("expq");
+
return exposants;
}
static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
@@ -760,19 +762,19 @@ static fullMatrix<double> generatePointsPrism(int order, bool serendip)
return point;
}
-static fullMatrix<double> generatePointsQuad(int order, bool serendip)
+static fullMatrix<double> generatePointsQuadRecur(int order, bool serendip)
{
int nbPoints = serendip ? order*4 : (order+1)*(order+1);
fullMatrix<double> point(nbPoints, 2);
if (order > 0) {
- point(0, 0) = 0;
- point(0, 1) = 0;
+ point(0, 0) = -1;
+ point(0, 1) = -1;
point(1, 0) = 1;
- point(1, 1) = 0;
+ point(1, 1) = -1;
point(2, 0) = 1;
point(2, 1) = 1;
- point(3, 0) = 0;
+ point(3, 0) = -1;
point(3, 1) = 1;
if (order > 1) {
@@ -786,8 +788,8 @@ static fullMatrix<double> generatePointsQuad(int order, bool serendip)
point(index, 1) = point(p0, 1) + i*(point(p1,1)-point(p0,1))/order;
}
}
- if (order > 2 && !serendip) {
- fullMatrix<double> inner = generatePointsQuad(order - 2, false);
+ if (order >= 2 && !serendip) {
+ fullMatrix<double> inner = generatePointsQuadRecur(order - 2, false);
inner.scale(1. - 2./order);
point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
}
@@ -797,9 +799,22 @@ static fullMatrix<double> generatePointsQuad(int order, bool serendip)
point(0, 0) = 0;
point(0, 1) = 0;
}
+
+
+ return point;
+}
+
+static fullMatrix<double> generatePointsQuad(int order, bool serendip){
+ fullMatrix<double> point = generatePointsQuadRecur(order, serendip);
+ // rescale to [0,1] x [0,1]
+ for (int i=0;i<point.size1();i++){
+ point(i,0) = (1.+point(i,0))*.5;
+ point(i,1) = (1.+point(i,1))*.5;
+ }
return point;
}
+
// Sub Control Points
static std::vector< fullMatrix<double> > generateSubPointsLine(int order)
{
@@ -1231,6 +1246,7 @@ const bezierBasis *bezierBasis::find(int tag)
B.points = generatePointsQuad(F->order,false);
dimSimplex = 0;
subPoints = generateSubPointsQuad(F->order);
+ // B.points.print("points");
break;
}
case TYPE_PRI : {
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index b90a654082cd5370004bf22f6baa16cfb761ba60..64aefa92e5df859cff75ff074a05df35b9283723 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -522,7 +522,7 @@ void FindCubicRoots(const double coef[4], double real[3], double imag[3])
double d = coef[0];
if(!a || !d){
- Msg::Error("Degenerate cubic: use a second degree solver!");
+ // Msg::Error("Degenerate cubic: use a second degree solver!");
return;
}
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index 70ba9b91ff4a8d860daabdb9c452731a6d502651..6d1b45d0d4653271f039e6b3dda161886ff8126d 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -1,6 +1,6 @@
-lc = 0.2 ;
-lc2 = 10 ;
-lc3 = 0.1 ;
+lc = .033 ;
+lc2 = 2.2 ;
+lc3 = .03 ;
Point(1) = {1.000000e+00,0.000000e+00,0.000000e+00,lc3};
Point(2) = {9.997533e-01,0.000000e+00,-3.498543e-05,lc};
Point(3) = {9.990134e-01,0.000000e+00,-1.398841e-04,lc};
diff --git a/benchmarks/2d/wing-splines.geo b/benchmarks/2d/wing-splines.geo
index 0aea8841a1e4b89961e10d51b4c16029446fc027..49d30ddc7b75270e0cac582bb4c32fcaee9c1a20 100644
--- a/benchmarks/2d/wing-splines.geo
+++ b/benchmarks/2d/wing-splines.geo
@@ -1,7 +1,7 @@
scale = .01 ;
-lc_wing = 0.001 * scale ;
+lc_wing = 0.05 * scale ;
lc_box = 10 * scale ;
Point(3895) = {1.177410e-02*scale,-2.768003e-03*scale,0,lc_wing};
@@ -520,32 +520,32 @@ Point(4) = {-4.447736e+01*scale,-4.504225e+01*scale,0,lc_box};
Point(1) = {4.552264e+01*scale,-4.504225e+01*scale,0,lc_box};
//one
-Spline(1) = {3846,3981,4001,3942,3940,3889,3892,3891,3943,3895,
+Spline(1) = {3846,3981,4001,3942,3940,3892,3943,3895,
3897,3896,3968,3995,4003,3857,3856,3860,3861,3863,
- 3864,3865,3866,3867,3868,3869,3870,3871,3872,3977,
+ 3864,3865,3866,3867,3869,3870,3871,3872,3977,
3877,3876,3878,3934,3873,3874,3875,3935,3880,3879,
3881,3936,3882,3883,3885,3884,1218,3933,3996,3989,
- 3990,3978,3979,3974,3973,3963,3947,3948,3904,3903,
+ 3990,3978,3979,3974,3973,3963,3948,3904,3903,
3946,3902,3901,3900,3908,3907,3951,3950,3847};
Spline(2) = {3847,3949,3952,3905,3906,3909,3969,3970,3997,3998,
4004,3959,3960,3972,3984,3988,3961,3962,3937,3938,
3886,3887,3888,3993,3994,3971,3918,3919,3920,3956,
- 3955,3965,3966,3910,3913,3912,3911,3914,3915,3916,
+ 3965,3910,3913,3912,3914,3915,3916,
3917,3953,3954,3964,3975,3992,3991,3999,3939,3967,
- 3987,3985,3986,4002,3922,4005,3851,3850,3853,3852,
- 3849,3848,3921,3983,3982,3976,3929,3958,3928,3957,
- 3927,3926,3925,3924,3932,3931,3930,3923,3862,3859,
+ 3987,3986,4002,3922,4005,3851,3850,3853,3852,
+ 3849,3848,3921,3983,3982,3976,3929,3958,3957,
+ 3927,3926,3925,3924,3931,3930,3923,3862,3859,
3858,3855,3854,3945,3944,3899,3898,3894,3890,3893,
3941,4000,3980,3846};
// two
Spline(3) = {4063,4062,4061,4114,4153,4178,4179,4191,4203,4109,
- 4202,4199,4201,4200,4100,4099,4045,4044,4043,4147,
+ 4202,4201,4200,4100,4099,4045,4044,4043,4147,
4146,4180,4167,4166,4134,4133,4098,4097,4042,4041,
4040,4095,4096,4132,4131,4165,4164,4163,4130,4129,
4094,4093,4039,4038,4037,4092,4091,4128,4127,4162,
4161,4183,4160,4159,4126,4125,4088,4087,4033,4032,
- 4032,4031,4086,4085,4124,4123,4158,4157,4156,4122,
+ 4031,4086,4085,4124,4123,4158,4157,4156,4122,
4121,4081,4080,4020,4019,4018,4143,4118,4022,4016,
4015,4011,4010,4009,4219,4177,4176,4149,4059,4058,
4055,4113,4110,4213,4210,4217,4006};
@@ -560,24 +560,24 @@ Spline(6) = {4008,4030,4084,4194,4195,4154,4155,4119,4120,4089,
4106,4052,4053,4054,4107,4108,4141,4142,4174,4175,
4187,4188,4145,4046,4047,4048,4196,4197,4198,4189,
4190,4214,4181,4182,4151,4152,4064,4065,4066,4205,
- 4206,4207,4072,4074,4073,4076,4075,4079,4078,4077,
+ 4206,4207,4072,4074,4076,4079,4078,4077,
4068,4067,4069,4070,4071,4115,4116,4204,4063};
//three
Spline(10) = {4298,4297,4222};
Spline(11) = {4222,4250,4252,4221};
-Spline(12) = {4221,4251,4299,4335,4337,4267,4294,4269,4233,4234,
- 4273,4293,4314,4317,4317,4318,4318,4274,4235,4236,
- 4237,4301,4319,4321,4320,4322,4347,4323,4324,4325,
- 4276,4238,4240,4239,4241,4277,4326,4327,4328,4350,
+Spline(12) = {4221,4251,4299,4335,4337,4267,4269,4233,4234,
+ 4273,4293,4314,4317,4318,4274,4235,4236,
+ 4237,4301,4319,4321,4322,4347,4323,4324,4325,
+ 4276,4238,4240,4241,4277,4326,4327,4328,4350,
4348,4349,4332,4333,4334,4278,4247,4248,4249,4282,
4283,4341,4342,4343,4289,4290,4286,4255,4257,4256,
- 4346,4288,4287,4260,4259,4258,4262,4261,4264,4264,
- 4263,4265,4266,4302,4304,4305,4306,4307,4308,
+ 4346,4288,4287,4260,4259,4258,4262,4261,4264,
+ 4263,4265,4266,4304,4307,
4309,4345,4344,4254,4253,4285,4284,4340,4339,4338,
4246,4245,4281,4280,4279,4244,4243,4242,4296,4295,
- 4300,4330,4331,4329,4232,4231,4230,4229,4228,4227,
- 4225,4226,4312,4223,4224,4311,4310,4275,4316,4315,
+ 4300,4331,4329,4232,4231,4230,4229,4228,4227,
+ 4225,4226,4312,4224,4311,4310,4275,4316,4315,
4313,4292,4272,4271,4270,4268,4291,4336,4298};
//box
diff --git a/benchmarks/2d/wing.geo b/benchmarks/2d/wing.geo
index febf083bd3a23fbcf2f8389d6d62681ccb3d8841..8cc01c722036375e4251d55e9ba345b0cb43c3fe 100644
--- a/benchmarks/2d/wing.geo
+++ b/benchmarks/2d/wing.geo
@@ -1,6 +1,6 @@
scale = 2 ;
-lc1 = 5.e-3 *scale ;
+lc1 = 15.e-3 *scale ;
lc2 = 1.e-2 *scale ;
lc3 = 10 *scale ;