From 13eacc55a0429c14c8f31099ee8783bbc484181e Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 26 Sep 2008 13:38:52 +0000
Subject: [PATCH] added embedded points in 2D mesh
---
Fltk/GUI.cpp | 2 +-
Geo/GFace.h | 3 ++
Geo/MElement.cpp | 13 +++++-
Geo/MElement.h | 6 ++-
Mesh/BDS.cpp | 5 +-
Mesh/HighOrder.cpp | 90 +++++++++++++++++++++++++-----------
Mesh/meshGFace.cpp | 12 ++++-
Mesh/meshGFaceBDS.cpp | 27 +++++++++--
Mesh/meshGFaceOptimize.cpp | 11 +++++
Mesh/qualityMeasures.cpp | 13 ++++--
benchmarks/2d/Square-Emb.geo | 5 ++
benchmarks/2d/naca12_2d.geo | 10 ++--
12 files changed, 151 insertions(+), 46 deletions(-)
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index dab66cf2d4..d855f08e56 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -2455,7 +2455,7 @@ void GUI::create_option_window()
mesh_value[3]->align(FL_ALIGN_RIGHT);
mesh_value[3]->callback(mesh_options_ok_cb);
- mesh_butt[4] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 9 * BH, BW, BH, "Use incomplete second order elements");
+ mesh_butt[4] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 9 * BH, BW, BH, "Use incomplete high order elements");
mesh_butt[4]->type(FL_TOGGLE_BUTTON);
mesh_butt[4]->callback(mesh_options_ok_cb);
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 3f7948bea1..0c569cdc5b 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -68,6 +68,9 @@ class GFace : public GEntity
// edges that are embedded in the face
virtual std::list<GEdge*> embeddedEdges() const { return embedded_edges; }
+ // edges that are embedded in the face
+ virtual std::list<GVertex*> embeddedVertices() const { return embedded_vertices; }
+
// vertices that bound the face
virtual std::list<GVertex*> vertices() const;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index eda15b5d7c..b3d5a05eba 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -160,6 +160,17 @@ double MTetrahedron::distoShapeMeasure()
#endif
}
+double MTetrahedronN::distoShapeMeasure()
+{
+#if defined(HAVE_GMSH_EMBEDDED)
+ return 1.;
+#else
+ if (_disto < -1.e21)_disto = qmDistorsionOfMapping(this);
+ return _disto;
+#endif
+}
+
+
double MTriangle::gammaShapeMeasure()
{
#if defined(HAVE_GMSH_EMBEDDED)
@@ -893,7 +904,7 @@ const gmshFunctionSpace* MLine::getFunctionSpace(int o) const {
return NULL;
}
-const int numSubEdges = 6;
+const int numSubEdges = 12;
int MTriangleN::getNumFacesRep(){ return numSubEdges * numSubEdges; }
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 3ca34e483c..162bf45198 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -1403,20 +1403,22 @@ class MTetrahedronN : public MTetrahedron {
protected:
std::vector<MVertex *> _vs;
const short _order;
+ double _disto;
public:
MTetrahedronN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
std::vector<MVertex*> &v, int order, int num=0, int part=0)
- : MTetrahedron(v0, v1, v2, v3, num, part) , _vs (v), _order(order)
+ : MTetrahedron(v0, v1, v2, v3, num, part) , _vs (v), _order(order),_disto(-1.e22)
{
for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
}
MTetrahedronN(std::vector<MVertex*> &v, int order, int num=0, int part=0)
- : MTetrahedron(v[0], v[1], v[2], v[3], num, part) , _order(order)
+ : MTetrahedron(v[0], v[1], v[2], v[3], num, part) , _order(order),_disto(-1.e22)
{
for(unsigned int i = 4; i < v.size(); i++) _vs.push_back(v[i]);
for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
}
~MTetrahedronN(){}
+ virtual double distoShapeMeasure();
virtual int getPolynomialOrder() const { return _order; }
virtual int getNumVertices() const { return 4 + _vs.size(); }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 1f730084fa..4104574066 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1331,7 +1331,10 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
{
if(!p->config_modified) return false;
if(p->g && p->g->classif_degree <= 1) return false;
-
+ if(p->g && p->g->classif_tag < 0) {
+ p->config_modified = true;
+ return true;
+ }
std::list<BDS_Edge*>::iterator eit = p->edges.begin();
while(eit != p->edges.end()) {
if((*eit)->numfaces() == 1) return false;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 54793e1b86..5790393e32 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -538,7 +538,8 @@ bool reparamOnEdge(MVertex *v, GEdge *ge, double ¶m)
}
void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
- edgeContainer &edgeVertices, bool linear, int nPts = 1)
+ edgeContainer &edgeVertices, bool linear, int nPts = 1,
+ std::map<MVertex*,SVector3> *disp = 0)
{
for(int i = 0; i < ele->getNumEdges(); i++){
MEdge edge = ele->getEdge(i);
@@ -588,7 +589,11 @@ void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
}
else {
GPoint pc = ge->point(US[j+1]);
- v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge, US[j+1]);
+ v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge, US[j+1]);
+ if (disp){
+ SPoint3 pc2 = edge.interpolate(t);
+ (*disp)[v] = SVector3(pc2.x(),pc2.y(),pc2.z());
+ }
}
temp.push_back(v);
ge->mesh_vertices.push_back(v);
@@ -603,7 +608,8 @@ void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
}
void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
- edgeContainer &edgeVertices, bool linear, int nPts = 1)
+ edgeContainer &edgeVertices, bool linear, int nPts = 1,
+ std::map<MVertex*,SVector3> *disp = 0)
{
for(int i = 0; i < ele->getNumEdges(); i++){
MEdge edge = ele->getEdge(i);
@@ -638,7 +644,11 @@ void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
}
else{
GPoint pc = gf->point(US[j+1], VS[j+1]);
- v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, US[j+1], VS[j+1]);
+ v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, US[j+1], VS[j+1]);
+ if (disp){
+ SPoint3 pc2 = edge.interpolate(t);
+ (*disp)[v] = SVector3(pc2.x(),pc2.y(),pc2.z());
+ }
}
temp.push_back(v);
gf->mesh_vertices.push_back(v);
@@ -652,11 +662,13 @@ void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
}
}
-void getEdgeVertices(GRegion *gr, MElement *ele,
- std::vector<MVertex*> &ve,
+void getEdgeVertices(GRegion *gr, MElement *ele,
+ std::vector<MVertex*> &ve,
std::set<MVertex*>& blocked,
- edgeContainer &edgeVertices, bool linear,
- int nPts = 1)
+ edgeContainer &edgeVertices,
+ bool linear,
+ int nPts = 1,
+ std::map<MVertex*,SVector3> *displ = 0)
{
for(int i = 0; i < ele->getNumEdges(); i++){
MEdge edge = ele->getEdge(i);
@@ -693,11 +705,13 @@ void getFaceVertices(GFace *gf,
MElement *ele,
std::vector<MVertex*> &vf,
faceContainer &faceVertices,
- bool linear, int nPts = 1) {
+ bool linear, int nPts = 1,
+ std::map<MVertex*,SVector3> *displ = 0) {
Double_Matrix points;
int start = 0;
+
switch (nPts){
case 2 :
points = gmshFunctionSpaces::find(MSH_TRI_10).points;
@@ -777,14 +791,20 @@ void getFaceVertices(GFace *gf,
}
if (reparamOK){
GPoint gp = gf->closestPoint(SPoint3(X,Y,Z),GUESS);
- if (gp.g())
+ if (gp.g()){
v = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
- else
+ }
+ else{
v = new MVertex(X,Y,Z, gf);
+ }
}
else{
v = new MVertex(X,Y,Z, gf);
}
+ if (displ){
+ SPoint3 pc2 = face.interpolate(t1, t2);
+ (*displ)[v] = SVector3(pc2.x(),pc2.y(),pc2.z());
+ }
}
// should be expensive -> induces a new search each time
// faceVertices[p].push_back(v);
@@ -1200,13 +1220,13 @@ void getRegionVertices(GRegion *gr,
void setHighOrder(GEdge *ge, edgeContainer &edgeVertices, bool linear,
- int nbPts = 1)
+ int nbPts = 1, std::map<MVertex*,SVector3> *displ = 0)
{
std::vector<MLine*> lines2;
for(unsigned int i = 0; i < ge->lines.size(); i++){
MLine *l = ge->lines[i];
std::vector<MVertex*> ve;
- getEdgeVertices(ge, l, ve, edgeVertices, linear, nbPts);
+ getEdgeVertices(ge, l, ve, edgeVertices, linear, nbPts, displ);
if(nbPts == 1)
lines2.push_back(new MLine3(l->getVertex(0), l->getVertex(1), ve[0]));
else
@@ -1219,13 +1239,13 @@ void setHighOrder(GEdge *ge, edgeContainer &edgeVertices, bool linear,
void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
faceContainer &faceVertices, bool linear, bool incomplete,
- int nPts = 1)
+ int nPts = 1, std::map<MVertex*,SVector3> *displ = 0)
{
std::vector<MTriangle*> triangles2;
for(unsigned int i = 0; i < gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i];
std::vector<MVertex*> ve, vf;
- getEdgeVertices(gf, t, ve, edgeVertices, linear, nPts);
+ getEdgeVertices(gf, t, ve, edgeVertices, linear, nPts,displ);
if(nPts == 1){
triangles2.push_back
(new MTriangle6(t->getVertex(0), t->getVertex(1), t->getVertex(2),
@@ -1244,7 +1264,7 @@ void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
else{
MTriangleN incpl (t->getVertex(0), t->getVertex(1), t->getVertex(2),
ve, nPts + 1);
- getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts);
+ getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts,displ);
}
ve.insert(ve.end(), vf.begin(), vf.end());
triangles2.push_back
@@ -1260,7 +1280,7 @@ void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
MQuadrangle *q = gf->quadrangles[i];
std::vector<MVertex*> ve, vf;
- getEdgeVertices(gf, q, ve, edgeVertices, linear, nPts);
+ getEdgeVertices(gf, q, ve, edgeVertices, linear, nPts,displ);
if(incomplete){
quadrangles2.push_back
(new MQuadrangle8(q->getVertex(0), q->getVertex(1), q->getVertex(2),
@@ -1280,7 +1300,7 @@ void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
faceContainer &faceVertices, bool linear, bool incomplete,
- int nPts = 1)
+ int nPts = 1, std::map<MVertex*,SVector3> *displ = 0)
{
int nbCorr = 0;
@@ -1291,7 +1311,7 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
std::set<MVertex*> blocked;
std::vector<MVertex*> ve,vf,vr;
- getEdgeVertices(gr, t, ve, blocked, edgeVertices, linear, nPts);
+ getEdgeVertices(gr, t, ve, blocked,edgeVertices, linear, nPts,displ);
if (nPts == 1){
tetrahedra2.push_back
(new MTetrahedron10(t->getVertex(0), t->getVertex(1), t->getVertex(2),
@@ -1343,7 +1363,7 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
MHexahedron *h = gr->hexahedra[i];
std::vector<MVertex*> ve, vf;
std::set<MVertex*> blocked;
- getEdgeVertices(gr, h, ve, blocked, edgeVertices, linear, nPts);
+ getEdgeVertices(gr, h, ve, blocked, edgeVertices, linear, nPts,displ);
if(incomplete){
hexahedra2.push_back
(new MHexahedron20(h->getVertex(0), h->getVertex(1), h->getVertex(2),
@@ -1373,7 +1393,7 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
MPrism *p = gr->prisms[i];
std::vector<MVertex*> ve, vf;
std::set<MVertex*> blocked;
- getEdgeVertices(gr, p, ve, blocked,edgeVertices, linear, nPts);
+ getEdgeVertices(gr, p, ve, blocked, edgeVertices, linear, nPts,displ);
if(incomplete){
prisms2.push_back
(new MPrism15(p->getVertex(0), p->getVertex(1), p->getVertex(2),
@@ -1397,7 +1417,7 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
MPyramid *p = gr->pyramids[i];
std::vector<MVertex*> ve, vf;
std::set<MVertex*> blocked;
- getEdgeVertices(gr, p, ve, blocked, edgeVertices, linear, nPts);
+ getEdgeVertices(gr, p, ve, blocked, edgeVertices, linear, nPts,displ);
if(incomplete){
pyramids2.push_back
(new MPyramid13(p->getVertex(0), p->getVertex(1), p->getVertex(2),
@@ -2047,21 +2067,39 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
double t1 = Cpu();
// first, make sure to remove any existsing second order vertices/elements
- SetOrder1(m);
+ SetOrder1(m);
+
+ std::map<MVertex*,SVector3> *displ = 0;
+ if(CTX.mesh.smooth_internal_edges){
+ displ = new std::map<MVertex*,SVector3>;
+ }
// then create new second order vertices/elements
edgeContainer edgeVertices;
faceContainer faceVertices;
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
- setHighOrder(*it, edgeVertices, linear, nPts);
+ setHighOrder(*it, edgeVertices, linear, nPts,displ);
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
- setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
+ setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts,displ);
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
- setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
+ setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts,displ);
printJacobians(m, "detjIni.pos");
if(CTX.mesh.smooth_internal_edges){
+ FILE *fl = fopen("displ.dat","w");
+ std::map<MVertex*,SVector3>::iterator it = displ->begin();
+ fprintf(fl,"%d\n",displ->size());
+ for ( ; it!=displ->end();++it){
+ fprintf(fl,"%6d %22.15E %22.15E %22.15E %22.15E %22.15E %22.15E\n",(it->first)->getNum(),(it->first)->x(),(it->first)->y(),(it->first)->z(),
+ it->second.x(),it->second.y(),it->second.z());
+ }
+ fclose(fl);
+ delete displ;
+ }
+
+
+ if(0 && CTX.mesh.smooth_internal_edges){
checkHighOrderTriangles(m);
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
Msg::Info("Smoothing internal Edges in Surface %d",(*it)->tag());
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 206f3f4c06..670f307173 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -322,13 +322,16 @@ static bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
{
+ debug = true;
BDS_GeomEntity CLASS_F (1, 2);
typedef std::set<MVertex*> v_container;
v_container all_vertices;
std::map<int, MVertex*>numbered_vertices;
std::list<GEdge*> edges = gf->edges();
std::list<GEdge*> emb_edges = gf->embeddedEdges();
+ std::list<GVertex*> emb_vertx = gf->embeddedVertices();
std::list<GEdge*>::iterator it = edges.begin();
+ std::list<GVertex*>::iterator itvx = emb_vertx.begin();
// build a set with all points of the boundaries
it = edges.begin();
@@ -357,7 +360,14 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
}
++it;
}
-
+
+ // add embedded vertices
+ while(itvx != emb_vertx.end()){
+ all_vertices.insert((*itvx)->mesh_vertices.begin(),
+ (*itvx)->mesh_vertices.end() );
+ ++itvx;
+ }
+
if (all_vertices.size() < 3){
Msg::Warning("Mesh Generation of Model Face %d Skipped: "
"Only %d Mesh Vertices on The Contours",
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index ee315ea877..ac6008b8d2 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -519,8 +519,25 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
int MAXNP = m.MAXPOINTNUMBER;
+
+ // classify correctly the embedded vertices
+ // use a negative model face number to avoid
+ // mesh motion
+ std::list<GVertex*> emb_vertx = gf->embeddedVertices();
+ std::list<GVertex*>::iterator itvx = emb_vertx.begin();
+ while(itvx != emb_vertx.end()){
+ MVertex *v = *((*itvx)->mesh_vertices.begin());
+ BDS_Point *p = m.find_point(v->getIndex());
+ m.add_geom (-1, 2);
+ p->g = m.get_geom(-1,2);
+ p->lc() = (*itvx)->prescribedMeshSizeAtVertex();
+ p->lcBGM() = (*itvx)->prescribedMeshSizeAtVertex();
+ ++itvx;
+ }
+
// IF ASKED , compute nodal size field using 1D Mesh
if (computeNodalSizeField){
+
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while (itp != m.points.end()){
std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
@@ -535,10 +552,12 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
}
++it;
}
- if (!ne) L = 1.e22;
- if(CTX.mesh.lc_from_points)
- (*itp)->lc() = L;
- (*itp)->lcBGM() = L;
+ if ((*itp)->g && (*itp)->g->classif_tag > 0){
+ if (!ne) L = 1.e22;
+ if(CTX.mesh.lc_from_points)
+ (*itp)->lc() = L;
+ (*itp)->lcBGM() = L;
+ }
++itp;
}
}
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 376dd967a4..240502049f 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -81,6 +81,17 @@ void buidMeshGenerationDataStructures(GFace *gf, std::set<MTri3*, compareTri3Ptr
for (unsigned int i = 0;i < gf->triangles.size(); i++)
setLcs(gf->triangles[i], vSizesMap);
+ // take care of embedded vertices
+ {
+ std::list<GVertex*> emb_vertx = gf->embeddedVertices();
+ std::list<GVertex*>::iterator itvx = emb_vertx.begin();
+ while(itvx != emb_vertx.end()){
+ MVertex *v = *((*itvx)->mesh_vertices.begin());
+ vSizesMap[v] = std::min(vSizesMap[v],(*itvx)->prescribedMeshSizeAtVertex());
+ ++itvx;
+ }
+ }
+
int NUM = 0;
for (std::map<MVertex*, double>::iterator it = vSizesMap.begin();
it != vSizesMap.end(); ++it){
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 5b38f8bfad..7eb0bf020b 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -209,6 +209,7 @@ static double mesh_functional_distorsion(MTriangle *t, double u, double v)
double qmDistorsionOfMapping (MTriangle *e)
{
+ return 1.0;
if (e->getPolynomialOrder() == 1)return 1.0;
IntPt *pts;
int npts;
@@ -264,11 +265,13 @@ double qmDistorsionOfMapping (MTetrahedron *e)
const double di = mesh_functional_distorsion (e,u,v,w);
dmin = (i==0)? di : std::min(dmin,di);
}
- double p[4][3] = {{0,0,0},{0,1,0},{1,0,0},{0,0,1}};
- for (int i=0;i<4;i++){
- const double u = p[i][0];
- const double v = p[i][1];
- const double w = p[i][2];
+
+ const Double_Matrix& points = e->getFunctionSpace()->points;
+
+ for (int i=0;i<e->getNumPrimaryVertices();i++) {
+ const double u = points(i,0);
+ const double v = points(i,1);
+ const double w = points(i,2);
const double di = mesh_functional_distorsion (e,u,v,w);
dmin = std::min(dmin,di);
}
diff --git a/benchmarks/2d/Square-Emb.geo b/benchmarks/2d/Square-Emb.geo
index d6d70c7138..73cd80b8af 100644
--- a/benchmarks/2d/Square-Emb.geo
+++ b/benchmarks/2d/Square-Emb.geo
@@ -22,3 +22,8 @@ Spline(5) = {11,22,33,44};
Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
Line {5} In Surface {6};
+
+Point(111) = {0.5,0.5,0,lc2/10};
+Point(112) = {0.53,0.5,0,lc2/10};
+Point(113) = {0.5,0.51,0,lc2/10};
+Point {111,112,113} In Surface {6};
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index dbc2c9f090..3b2bd214b3 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -1,6 +1,6 @@
lc = 0.2;
-lc2 = 10;
-lc3 = 0.02;
+lc2 = 1;
+lc3 = 0.1;
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};
@@ -205,12 +205,12 @@ Point(200) = {9.997533e-01,0.000000e+00,3.498543e-05,lc};
Spline(1) = { 1 ... 50};
Spline(2) = { 50 ... 101};
Spline(3) = { 101 ... 151};
-Spline(4) = { 151 ... 200, 1};
+Spline(4) = { 151 ... 200,1};
Rotate { {1,0,0},{0,0,0},Pi/2 } { Line{1,2,3,4}; }
Translate {-0.5,0,0} { Line{1,2,3,4}; }
-d=10;
+d=4;
Point(1000) = {d,d,0,lc2};
Point(1001) = {-d,d,0,lc2};
Point(1002) = {-d,-d,0,lc2};
@@ -224,7 +224,7 @@ Line Loop(9) = {6,7,8,5};
Line Loop(10) = {2,3,4,1};
Plane Surface(11) = {9,10};
-Point(9999) = {0.6,0,0,1};
+//Point(9999) = {0.6,0,0,1};
//Recombine Surface {11};
--
GitLab