diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index c16ad6482a21c8567a3a8a58fd3a868432f8da17..a38d80711e8ed7b841b78a246dba1916a823f162 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.87 2006-12-16 15:44:27 geuzaine Exp $
+// $Id: CommandLine.cpp,v 1.88 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -406,8 +406,10 @@ void Get_Options(int argc, char *argv[])
else if(!strcmp(argv[i] + 1, "algo")) {
i++;
if(argv[i] != NULL) {
- if(!strncmp(argv[i], "iso", 3))
- CTX.mesh.algo2d = DELAUNAY_ISO;
+ if(!strncmp(argv[i], "bds", 3))
+ CTX.mesh.algo2d = MESHADAPT;
+ else if(!strncmp(argv[i], "del", 3))
+ CTX.mesh.algo2d = DELAUNAY2D;
else if(!strncmp(argv[i], "tri", 3))
CTX.mesh.algo2d = DELAUNAY_TRIANGLE;
else if(!strncmp(argv[i], "netgen", 6))
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 9625829a94a28f572f051db44ccd3b0a65468a13..08f9c12ece52727bdb3b14126b969c6127b60b82 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -836,8 +836,8 @@ StringXNumber GeometryOptions_Number[] = {
} ;
StringXNumber MeshOptions_Number[] = {
- { F|O, "Algorithm" , opt_mesh_algo2d , DELAUNAY_ISO ,
- "2D mesh algorithm (1=isotropic, 2=anisotropic, 3=triangle)" },
+ { F|O, "Algorithm" , opt_mesh_algo2d , MESHADAPT ,
+ "2D mesh algorithm (1=meshadapt, 2=delaunay)" },
{ F|O, "Algorithm3D" , opt_mesh_algo3d ,
#if defined(HAVE_TETGEN)
DELAUNAY_ISO,
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index d2e48de5c555bbb907e4dc40d438e9ad7d7b0445..a35873417a857a65b4be29061784fa8e747f2e22 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -66,6 +66,8 @@
#define DELAUNAY_TRIANGLE 3
#define FRONTAL_NETGEN 4
#define DELAUNAY_TETGEN 5
+#define MESHADAPT 6
+#define DELAUNAY2D 7
#define TRANSFINI 1
#define LIBRE 2
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index c868fa05010b45a31ea73dcdba2442de7d82dd8c..922c46c399863a96ef053b008916d489bbaaded1 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.499 2006-12-18 19:47:37 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.500 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -1092,7 +1092,8 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
opt_mesh_point_type(0, GMSH_SET, WID->mesh_choice[0]->value());
opt_mesh_algo2d(0, GMSH_SET,
- (WID->mesh_choice[2]->value() == 0) ? DELAUNAY_ISO :
+ (WID->mesh_choice[2]->value() == 0) ? MESHADAPT :
+ (WID->mesh_choice[2]->value() == 1) ? DELAUNAY2D :
DELAUNAY_TRIANGLE);
opt_mesh_recombine_algo(0, GMSH_SET,
(WID->mesh_choice[5]->value() == 0) ? 1 : 2);
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index bd7b1dde0e267eeaea11d762b2c2ab7fd6453c96..80bce128b52f37e6044ea90eb64e568ecbeb2e44 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.588 2006-12-21 21:47:05 geuzaine Exp $
+// $Id: GUI.cpp,v 1.589 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -2229,7 +2229,8 @@ void GUI::create_option_window()
o->hide();
static Fl_Menu_Item menu_2d_algo[] = {
- {"Isotropic", 0, 0, 0},
+ {"MeshAdapt", 0, 0, 0},
+ {"Delaunay", 0, 0, 0},
{"Triangle", 0, 0, 0},
{0}
};
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 9ea6da85471a4969e2671b40aaff4279add25200..0d3f234024717052afe07aea505308f4de8c91a0 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.26 2006-12-17 22:23:17 geuzaine Exp $
+// $Id: GModel.cpp,v 1.27 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -301,6 +301,7 @@ double GModel::getMeshSize ()
if(meshSize < 0){
SBoundingBox3d b = bounds();
return 0.1 * norm(SVector3(b.max(), b.min()));
+
}
else
return meshSize;
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 77a3860e0c0340f8f46721b29247ce058e95f884..c4d4cfcd23c966e4309d83b220697fc7f3bb03b9 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Geo.cpp,v 1.5 2007-01-12 08:10:32 geuzaine Exp $
+// $Id: GModelIO_Geo.cpp,v 1.6 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -128,6 +128,8 @@ int GModel::importTHEM()
ge->physicals.push_back(pnum);
}
}
+
+ meshSize = 1.e22;
Msg(DEBUG, "Gmsh model imported:");
Msg(DEBUG, "%d Vertices", vertices.size());
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index bbc667e492a9f193b5d8fa3d464c261c49b76a78..2c8ce8ed8915957dfe885cbd85b2f0ccf1669040 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.25 2006-12-20 15:50:57 remacle Exp $
+// $Id: MElement.cpp,v 1.26 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -412,6 +412,7 @@ bool MTriangle::invertmappingXY(double *p, double *uv, double tol)
b[0] = p[0] - getVertex(0)->x();
b[1] = p[1] - getVertex(0)->y();
sys2x2(mat, b, uv);
+
if(uv[0] >= -tol &&
uv[1] >= -tol &&
uv[0] <= 1. + tol &&
@@ -464,6 +465,8 @@ void MTriangle::circumcenterXY(double *res) const
res[0] = (double)((a1 * (y3 - y2) + a2 * (y1 - y3) + a3 * (y2 - y1)) / d);
res[1] = (double)((a1 * (x2 - x3) + a2 * (x3 - x1) + a3 * (x1 - x2)) / d);
+ // printf("%g %g - %g %g - %g %g cc %g %g\n",x1,y1,x2,y2,x3,y3,res[0],res[1]);
+
return ;
}
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index e78d4f19701d0b022fa2af4a8b7d246cfc1685e4..aeb38945894c44cc80d35626d5e996af3be85493 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -1,4 +1,4 @@
-// $Id: BackgroundMesh.cpp,v 1.10 2006-12-05 14:22:05 remacle Exp $
+// $Id: BackgroundMesh.cpp,v 1.11 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -173,6 +173,7 @@ double BGM_MeshSize ( GEntity *ge, double U, double V , double X, double Y, doub
double l3 = ge->model()->getMeshSize();
double l4 = LC_MVertex_BGM ( ge, X, Y , Z );
if (ge->dim() < 2) l2 = LC_MVertex_PNTS ( ge, U, V );
+ // printf("ge->dim() %d l2 = %g l4 = %g l3 = %g\n",ge->dim(),l2,l3,l4);
double l = std::min(std::min(l2,l4),l3);
l *= CTX.mesh.lc_factor ;
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 3e62ac0a87b54a467770b76a6474a89a5fc62edc..6bc744d029610775147303d9a4dc82fd618f648d 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.158 2006-12-21 17:10:15 geuzaine Exp $
+# $Id: Makefile,v 1.159 2007-01-12 13:16:59 remacle Exp $
#
# Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
#
@@ -36,6 +36,7 @@ SRC = Generator.cpp \
meshGFace.cpp \
meshGFaceTransfinite.cpp \
meshGFaceExtruded.cpp \
+ meshGFaceDelaunayInsertion.cpp \
meshGRegion.cpp \
meshGRegionDelaunayInsertion.cpp \
meshGRegionTransfinite.cpp \
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index ccb2282aaefdceb675e058f651e169e4f86c2a65..c103dd28000499d2bf020e2eb4cee99815653105 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.25 2006-12-16 15:48:51 geuzaine Exp $
+// $Id: meshGEdge.cpp,v 1.26 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -56,6 +56,7 @@ double F_Lc(double t)
else
lc_here = BGM_MeshSize(_myGEdge, t , 0 , p.x(),p.y(),p.z());
SVector3 der = _myGEdge -> firstDer(t) ;
+
const double d = norm(der);
return d/lc_here;
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 2cdc22bdee303c899aed43a62f81b52a1fa606cd..340e0c82bf51580dfdebbd13fd77247c3fe94f48 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.45 2006-12-21 17:10:15 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.46 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -20,6 +20,7 @@
// Please report all bugs and problems to <gmsh@geuz.org>.
#include "meshGFace.h"
+#include "meshGFaceDelaunayInsertion.h"
#include "DivideAndConquer.h"
#include "BackgroundMesh.h"
#include "GVertex.h"
@@ -807,16 +808,14 @@ bool gmsh2DMeshGenerator ( GFace *gf )
m->del_point(m->find_point(-3));
m->del_point(m->find_point(-4));
- // goto hhh;
- // start mesh generation
-
- RefineMesh (gf,*m,100);
- // OptimizeMesh (gf,*m,2);
-
-
- if (gf->meshAttributes.recombine)
+ // start mesh generation
+ // if (CTX.mesh.algo2d == MESHADAPT)
{
- m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
+ RefineMesh (gf,*m,100);
+ if (gf->meshAttributes.recombine)
+ {
+ m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
+ }
}
// fill the small gmsh structures
@@ -829,7 +828,6 @@ bool gmsh2DMeshGenerator ( GFace *gf )
if (numbered_vertices.find(p->iD) == numbered_vertices.end())
{
MVertex *v = new MFaceVertex (p->X,p->Y,p->Z,gf,p->u,p->v);
- //MVertex *v = new MFaceVertex (p->u,p->v,0,gf,p->u,p->v);
numbered_vertices[p->iD]=v;
gf->mesh_vertices.push_back(v);
}
@@ -861,11 +859,19 @@ bool gmsh2DMeshGenerator ( GFace *gf )
}
}
+ // the delaunay algo is based directly on internal gmsh structures
+ // BDS mesh is passed in order not to recompute local coordinates
+ // of vertices
+// if (CTX.mesh.algo2d == DELAUNAY2D ||CTX.mesh.algo2d == DELAUNAY_ISO)
+// {
+// insertVerticesInFace (gf,m) ;
+// }
// delete the mesh
delete m;
+
delete [] U_;
delete [] V_;
@@ -1260,9 +1266,6 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
m->cleanup();
-
-
-
{
std::list<BDS_Edge*>::iterator ite = m->edges.begin();
while (ite != m->edges.end())
@@ -1291,16 +1294,14 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
// goto hhh;
// start mesh generation
- RefineMesh (gf,*m,100);
- // OptimizeMesh (gf,*m,2);
-
-
-
- if (gf->meshAttributes.recombine)
+ // if (CTX.mesh.algo2d == MESHADAPT)
{
- m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
+ RefineMesh (gf,*m,100);
+ if (gf->meshAttributes.recombine)
+ {
+ m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
+ }
}
-
// fill the small gmsh structures
{
@@ -1344,11 +1345,17 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
// delete the mesh
- char name[245];
- sprintf(name,"param%d.pos",gf->tag());
+ // char name[245];
+ // sprintf(name,"param%d.pos",gf->tag());
//outputScalarField(m->triangles, name,1);
- sprintf(name,"real%d.pos",gf->tag());
+ // sprintf(name,"real%d.pos",gf->tag());
//outputScalarField(m->triangles, name,0);
+
+// if (CTX.mesh.algo2d == DELAUNAY2D ||CTX.mesh.algo2d == DELAUNAY_ISO)
+// {
+// insertVerticesInFace (gf,m) ;
+// }
+
delete m;
return true;
@@ -1390,8 +1397,11 @@ void meshGFace :: operator() (GFace *gf)
Msg(DEBUG1, "Generating the mesh");
// temp fix until we create MEdgeLoops in gmshFace:
- if(gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty())
- gmsh2DMeshGenerator ( gf ) ;
+ if(gf->getNativeType() == GEntity::GmshModel ||
+ gf->edgeLoops.empty() )
+ {
+ gmsh2DMeshGenerator ( gf );
+ }
else
{
if (!gmsh2DMeshGeneratorPeriodic ( gf ))
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a08cd3ceb1fcb983378f0e365275836ff230d1c6
--- /dev/null
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -0,0 +1,569 @@
+// $Id: meshGFaceDelaunayInsertion.cpp,v 1.1 2007-01-12 13:16:59 remacle Exp $
+//
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "BDS.h"
+#include "BackgroundMesh.h"
+#include "meshGFaceDelaunayInsertion.h"
+#include "GFace.h"
+#include "Numeric.h"
+#include "Message.h"
+#include <set>
+#include <map>
+#include <algorithm>
+
+// computes the center of the circum circle in the tangent plane
+// the metric given by a b d is supposed to be constant
+void MTri3::Center_Circum_Aniso(double a, double b, double d, double &x, double &y, double &r) const
+{
+ double sys[2][2], X[2];
+ double rhs[2];
+ double x1, y1, x2, y2, x3, y3;
+
+ x1 = base->getVertex(0)->x();
+ y1 = base->getVertex(0)->y();
+ x2 = base->getVertex(1)->x();
+ y2 = base->getVertex(1)->y();
+ x3 = base->getVertex(2)->x();
+ y3 = base->getVertex(2)->y();
+
+ sys[0][0] = 2. * a * (x1 - x2) + 2. * b * (y1 - y2);
+ sys[0][1] = 2. * d * (y1 - y2) + 2. * b * (x1 - x2);
+ sys[1][0] = 2. * a * (x1 - x3) + 2. * b * (y1 - y3);
+ sys[1][1] = 2. * d * (y1 - y3) + 2. * b * (x1 - x3);
+
+ rhs[0] =
+ a * (x1 * x1 - x2 * x2) + d * (y1 * y1 - y2 * y2) + 2. * b * (x1 * y1 -
+ x2 * y2);
+ rhs[1] =
+ a * (x1 * x1 - x3 * x3) + d * (y1 * y1 - y3 * y3) + 2. * b * (x1 * y1 -
+ x3 * y3);
+ sys2x2(sys, rhs, X);
+
+ x = X[0];
+ y = X[1];
+
+ r = sqrt((X[0] - x1) * (X[0] - x1) * a
+ + (X[1] - y1) * (X[1] - y1) * d
+ + 2. * (X[0] - x1) * (X[1] - y1) * b);
+}
+
+
+// find a common basis for 2 metrics in 2D
+void simultaneousMetricReduction( const gmsh2dMetric &M1, const gmsh2dMetric &M2,
+ double & l1,double & l2, double V[2][2])
+{
+ double a11=M1.a11,a21=M1.a21,a22=M1.a22;
+ double b11=M2.a11,b21=M2.a21,b22=M2.a22;
+ const double /*c11 = a11*a11,*/ c21= a21*a21;
+ const double /*d11 = b11*b11,*/ d21= b21*b21;
+ const double a=b11*b22 - d21;
+ const double b=-a11*b22-a22*b11+2*a21*b21;
+ const double c=-c21+a11*a22;
+ const double bb = b*b,ac= a*c;
+ const double delta = bb - 4 * ac;
+ if (bb + fabs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) )
+ {
+ if (fabs(a) < 1.e-30 )
+ l1 = l2 = 0;
+ else
+ l1=l2=-b/(2*a);
+ V[0][0] = 1.;
+ V[1][0] = 0.;
+ V[0][1] = 0.;
+ V[1][1] = 1.;
+ }
+ else {
+ const double delta2 = sqrt(delta);
+ l1= (-b - delta2)/(2*a);
+ l2= (-b + delta2)/(2*a);
+ double v0 = a11-l1*b11, v1 = a21-l1*b21,v2 = a22 - l1*b22;
+ double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
+ double vp1x,vp1y,vp2x,vp2y;
+ if(s1 < s0)
+ s0=sqrt(s0),vp1x=v1/s0,vp1y=-v0/s0;
+ else
+ s1=sqrt(s1),vp1x=v2/s1,vp1y=-v1/s1;
+
+ v0 = a11-l2*b11, v1 = a21-l2*b21,v2 = a22 - l2*b22;
+ s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
+ if(s1 < s0)
+ s0=sqrt(s0),vp2x=v1/s0,vp2y=-v0/s0;
+ else
+ s1=sqrt(s1),vp2x=v2/s1,vp2y=-v1/s1;
+ V[0][0] = vp1x;
+ V[0][1] = vp2x;
+ V[1][0] = vp1y;
+ V[1][1] = vp2y;
+ }
+}
+
+gmsh2dMetric metricIntersection(const gmsh2dMetric &Ma,const gmsh2dMetric &Mb)
+{
+ double M[2][2],M1[2][2];
+ double l1,l2;
+ simultaneousMetricReduction(Ma,Mb,l1,l2,M);
+ inv2x2 ( M, M1 );
+ const double M1ev1 = Ma.eval(M[0][0],M[0][1]);
+ const double M1ev2 = Ma.eval(M[0][0],M[0][1]);
+ const double M2ev1 = Mb.eval(M[1][0],M[1][1]);
+ const double M2ev2 = Mb.eval(M[1][0],M[1][1]);
+
+ const double D0 = std::max( M1ev1,M1ev2);
+ const double D1 = std::max( M2ev1,M2ev2);
+
+ return gmsh2dMetric(M1[0][0] * D0 * M1[0][0] + M1[0][1] * D1 * M1[0][1] ,
+ 0.5 * ( M1[0][0] * D0 * M1[1][0] + M1[0][1] * D1 * M1[1][1] +
+ M1[1][0] * D0 * M1[0][0] + M1[1][1] * D1 * M1[0][1] ),
+ M1[1][0] * D0 * M1[1][0] + M1[1][1] * D1 * M1[1][1] );
+}
+
+gmsh2dMetric metricInterpolationTriangle(const gmsh2dMetric &M1,
+ const gmsh2dMetric &M2,
+ const gmsh2dMetric &M3,
+ const double u,
+ const double v)
+{
+ double m1[2][2] = {{M1.a11,0.5*M1.a21},{0.5*M1.a21,M1.a22}};
+ double m2[2][2] = {{M2.a11,0.5*M2.a21},{0.5*M2.a21,M2.a22}};
+ double m3[2][2] = {{M3.a11,0.5*M3.a21},{0.5*M3.a21,M3.a22}};
+ double invm1 [2][2]; inv2x2 ( m1, invm1 );
+ double invm2 [2][2]; inv2x2 ( m2, invm2 );
+ double invm3 [2][2]; inv2x2 ( m3, invm3 );
+ double invm [2][2] ;
+ for (int i=0;i<2;i++)
+ for (int j=0;j<2;j++)
+ invm[i][j] = (1.-u-v) * invm1[i][j] + u * invm2[i][j] + v * invm3[i][j];
+ double m [2][2]; inv2x2 ( invm, m );
+ return gmsh2dMetric ( m[0][0] , m[0][1] , m[1][1] );
+}
+
+gmsh2dMetric :: gmsh2dMetric ( double lc )
+ : a11 ( 1./(lc*lc) ) ,a21 ( 0.0 ) ,a22 ( 1./(lc*lc) )
+{
+}
+
+
+MTri3::MTri3 ( MTriangle * t, std::vector<gmsh2dMetric> & sizes) : deleted(false), base (t)
+{
+ neigh[0] = neigh[1] = neigh[2] = 0;
+ // compute the metric at the centroid
+ metric = metricInterpolationTriangle(sizes [base->getVertex(0)->getNum()] ,
+ sizes [base->getVertex(1)->getNum()] ,
+ sizes [base->getVertex(2)->getNum()] ,
+ 1./3., 1./3.);
+ Center_Circum_Aniso(metric.a11,metric.a21,metric.a22, xc, yc, circum_radius);
+ // printf("metric %g %g %g circum_radius = %g\n",metric.a11,metric.a21,metric.a22,circum_radius);
+}
+
+int inCircumCircle ( MTri3 *t , const double *p , const gmsh2dMetric &metric )
+{
+ double xc, yc;
+ double eps, d1, d2, x[2];
+ t->Center_Circum_Aniso(metric.a11,metric.a21,metric.a22,xc, yc, d1);
+
+ x[0] = xc - p[0];
+ x[1] = yc - p[1];
+
+ d2 = sqrt(fabs(metric.eval(x[0],x[1])));
+
+ eps = fabs(d1 - d2) / (d1 + d2);
+ if(eps < 1.e-12) {
+ return 0;
+ }
+ if(d2 < d1)
+ return 1;
+ return 0;
+
+}
+
+int MTri3::inCircumCircle ( const double *p ) const
+{
+
+// double pa[2] = {base->getVertex(0)->x(),base->getVertex(0)->y()};
+// double pb[2] = {base->getVertex(1)->x(),base->getVertex(1)->y()};
+// double pc[2] = {base->getVertex(2)->x(),base->getVertex(2)->y()};
+// double result = gmsh::incircle(pa, pb, pc, (double*)p) * gmsh::orient2d(pa, pb, pc);
+// return (result > 0) ? 1 : 0;
+
+
+ double eps, d1, d2, x[2];
+
+ x[0] = xc - p[0];
+ x[1] = yc - p[1];
+
+ d1 = circum_radius;
+ d2 = sqrt(fabs(metric.eval(x[0],x[1])));
+
+ eps = fabs(d1 - d2) / (d1 + d2);
+ if(eps < 1.e-12) {
+ return 0;
+ }
+ if(d2 < d1)
+ return 1;
+ return 0;
+}
+
+struct edgeXface
+{
+ MVertex *v[2];
+ MTri3 * t1;
+ int i1;
+ edgeXface ( MTri3 *_t, int iFac)
+ : t1(_t),i1(iFac)
+ {
+ v[0] = t1->tri()->getVertex ( iFac == 0 ? 2 : iFac-1 );
+ v[1] = t1->tri()->getVertex ( iFac );
+ std::sort ( v, v+2 );
+ }
+ inline bool operator < ( const edgeXface & other) const
+ {
+ if (v[0] < other.v[0])return true;
+ if (v[0] > other.v[0])return false;
+ if (v[1] < other.v[1])return true;
+ return false;
+ }
+};
+
+template <class ITER>
+void connectTris ( ITER beg, ITER end)
+{
+ std::set<edgeXface> conn;
+ {
+ while (beg != end)
+ {
+ if (!(*beg)->isDeleted())
+ {
+ for (int i=0;i<3;i++)
+ {
+ edgeXface fxt ( *beg , i );
+ std::set<edgeXface>::iterator found = conn.find (fxt);
+ if (found == conn.end())
+ conn.insert ( fxt );
+ else if (found->t1 != *beg)
+ {
+ found->t1->setNeigh(found->i1,*beg);
+ (*beg)->setNeigh ( i, found->t1);
+ }
+ }
+ }
+ ++beg;
+ }
+ }
+}
+
+
+void recurFindCavity (const gmsh2dMetric &metric,
+ std::list<edgeXface> & shell,
+ std::list<MTri3*> & cavity,
+ double *v,
+ MTri3 *t)
+{
+ t->setDeleted(true);
+ // the cavity that has to be removed
+ // because it violates delaunay criterion
+ cavity.push_back(t);
+
+ for (int i=0;i<3;i++)
+ {
+ MTri3 *neigh = t->getNeigh(i) ;
+ if (!neigh)
+ shell.push_back ( edgeXface ( t, i ) );
+ else if (!neigh->isDeleted())
+ {
+ int circ = inCircumCircle ( neigh, v, metric );
+ if (circ)
+ recurFindCavity ( metric,shell, cavity,v , neigh);
+ else
+ shell.push_back ( edgeXface ( t, i ) );
+ }
+ }
+}
+
+bool insertVertex (MVertex *v ,
+ MTri3 *t ,
+ std::set<MTri3*,compareTri3Ptr> &allTets,
+ std::vector<gmsh2dMetric> & vSizes,
+ const gmsh2dMetric &metric)
+{
+ std::list<edgeXface> shell;
+ std::list<MTri3*> cavity;
+ std::list<MTri3*> new_cavity;
+
+ double p[2] = {v->x(),v->y()};
+
+ recurFindCavity ( metric,shell, cavity, p , t);
+
+ // check that volume is conserved
+ double newVolume = 0;
+ double oldVolume = 0;
+
+ std::list<MTri3*>::iterator ittet = cavity.begin();
+ std::list<MTri3*>::iterator ittete = cavity.end();
+ while ( ittet != ittete )
+ {
+ oldVolume += fabs((*ittet)->getSurfaceXY());
+ ++ittet;
+ }
+
+// char name[245];
+// int III = 1;
+// sprintf(name,"test%d.pos",III);
+// FILE *ff = fopen (name,"w");
+// fprintf(ff,"View\"test\"{\n");
+
+
+ MTri3** newTris = new MTri3*[shell.size()];;
+ int k = 0;
+ std::list<edgeXface>::iterator it = shell.begin();
+
+ while (it != shell.end())
+ {
+ MTriangle *t = new MTriangle ( it->v[0],
+ it->v[1],
+ v);
+ MTri3 *t4 = new MTri3 ( t , vSizes);
+// fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
+// it->v[0]->x(),
+// it->v[0]->y(),
+// it->v[0]->z(),
+// it->v[1]->x(),
+// it->v[1]->y(),
+// it->v[1]->z(),
+// v->x(),
+// v->y(),
+// v->z());
+ newTris[k++]=t4;
+ // all new triangles are pushed front in order to
+ // ba able to destroy them if the cavity is not
+ // star shaped around the new vertex.
+ new_cavity.push_back (t4);
+ MTri3 *otherSide = it->t1->getNeigh(it->i1);
+
+ if (otherSide)
+ new_cavity.push_back (otherSide);
+ // if (!it->t1->isDeleted())throw;
+ double ss = fabs(t4->getSurfaceXY());
+ if (ss < 1.e-12)ss = 1256172121;
+ newVolume += ss;
+ ++it;
+ }
+// fprintf(ff,"};\n");
+// fclose (ff);
+// printf("%d %d newVolume %g oldVolume %g\n",cavity.size(),new_cavity.size(),newVolume,oldVolume);
+// getchar();
+
+
+ if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume )
+ {
+ connectTris ( new_cavity.begin(),new_cavity.end() );
+ allTets.insert(newTris,newTris+shell.size());
+ delete [] newTris;
+ return true;
+ }
+ // The cavity is NOT star shaped
+ else
+ {
+ for (unsigned int i=0;i<shell.size();i++)delete newTris[i];
+ delete [] newTris;
+ ittet = cavity.begin();
+ ittete = cavity.end();
+ while ( ittet != ittete )
+ {
+ (*ittet)->setDeleted(false);
+ ++ittet;
+ }
+ return false;
+ }
+}
+
+static gmsh2dMetric evalMetricTensor ( GFace *gf , MVertex *mv , double lc )
+{
+ Pair<SVector3,SVector3> der = gf->firstDer(SPoint2(mv->x(),mv->y())) ;
+ const double oneoverlc2 = 1./(lc*lc);
+ const double a11 = dot(der.first() ,der.first() ) * oneoverlc2;
+ const double a21 = dot(der.second(),der.first() ) * oneoverlc2;
+ const double a22 = dot(der.second(),der.second()) * oneoverlc2;
+
+ // printf("metric : %g %g %g - %g %g\n",a11,a21,a22,oneoverlc2,lc);
+
+ return gmsh2dMetric ( a11, a21 , a22);
+}
+
+static void setLcs ( MTriangle *t, std::map<MVertex*,double> &vSizes)
+{
+ for (int i=0;i<3;i++)
+ {
+ for (int j=i+1;j<3;j++)
+ {
+ MVertex *vi = t->getVertex(i);
+ MVertex *vj = t->getVertex(j);
+ double dx = vi->x()-vj->x();
+ double dy = vi->y()-vj->y();
+ double dz = vi->z()-vj->z();
+ double l = sqrt(dx*dx+dy*dy+dz*dz);
+ std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
+ if (iti==vSizes.end() || iti->second > l)vSizes[vi] = l;
+ if (itj==vSizes.end() || itj->second > l)vSizes[vj] = l;
+ }
+ }
+}
+
+void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
+{
+
+ std::set<MTri3*,compareTri3Ptr> AllTris;
+ std::map<MVertex*,double> vSizesMap;
+ std::map<const BDS_Point*,MVertex*> recoverMap;
+ std::vector<gmsh2dMetric> vMetrics;
+ std::vector<double> vSizes;
+
+
+ // compute edge sizes on the contour and propagate thos sizes inside the surface
+ for (unsigned int i=0;i<gf->triangles.size();i++)setLcs ( gf->triangles[i] , vSizesMap);
+
+ // set vertex positions to local coordinates of the face
+ // compute metric tensors of the surface
+ int NUM=0;
+ for (std::map<MVertex*,double>::iterator it = vSizesMap.begin();it!=vSizesMap.end();++it)
+ {
+ BDS_Point *p;
+ p = bds->find_point(it->first->getNum());
+ it->first->x() = p->u;
+ it->first->y() = p->v;
+ it->first->z() = 0;
+ it->first->setNum(NUM++);
+ vSizes.push_back(it->second);
+ recoverMap[p]= (it->first);
+ vMetrics.push_back(evalMetricTensor ( gf , it->first, it->second));
+ }
+
+ for (unsigned int i=0;i<gf->triangles.size();i++)
+ AllTris.insert ( new MTri3 ( gf->triangles[i] ,vMetrics ) );
+
+ gf->triangles.clear();
+ connectTris ( AllTris.begin(), AllTris.end() );
+
+ Msg(DEBUG,"All %d tris were connected",AllTris.size());
+
+ // here the classification should be done
+
+ int ITER = 0;
+
+ while (1)
+ {
+ MTri3 *worst = *AllTris.begin();
+
+ if (worst->isDeleted())
+ {
+ delete worst->tri();
+ delete worst;
+ AllTris.erase(AllTris.begin());
+ // Msg(INFO,"Worst tet is deleted");
+ }
+ else
+ {
+ if(ITER++%5000 ==0)
+ Msg(DEBUG,"%7d points created -- Worst tri radius is %8.3f",vSizes.size(),worst->getRadius());
+ if (worst->getRadius() < 0.5 * sqrt(2.0)) break;
+ double center[2],uv[2];
+ worst->getCenter(center);
+ // worst->tri()->circumcenterXY(center);
+ bool inside = worst->tri()->invertmappingXY(center,uv);
+ if (inside || 1)
+ {
+ // we use here local coordinates as real coordinates
+ // x,y and z will be computed hereafter
+ // Msg(INFO,"Point is inside");
+ MVertex *v = new MFaceVertex (center[0],center[1],0.0,gf,center[0],center[1]);
+ v->setNum(NUM++);
+ double lc1 = (1.-uv[0]-uv[1]) * vSizes [worst->tri()->getVertex(0)->getNum()] +
+ uv[0] * vSizes [worst->tri()->getVertex(1)->getNum()] +
+ uv[1] * vSizes [worst->tri()->getVertex(2)->getNum()] ;
+ GPoint p = gf->point (center[0],center[1]);
+ double lc = std::min(lc1,BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z()));
+ vSizes.push_back(lc);
+ gmsh2dMetric metr = evalMetricTensor ( gf , v, lc);
+ vMetrics.push_back(metr);
+
+ // compute mesh spacing there
+ if (!insertVertex ( v , worst, AllTris,vMetrics,metr))
+ {
+ AllTris.erase(AllTris.begin());
+ worst->forceRadius(0);
+ AllTris.insert(worst);
+ delete v;
+ }
+ else
+ {
+ gf->mesh_vertices.push_back(v);
+ }
+ }
+ else
+ {
+ // Msg(INFO,"Point is outside");
+ AllTris.erase(AllTris.begin());
+ worst->forceRadius(0);
+ AllTris.insert(worst);
+ }
+ }
+ }
+
+ // restore real coordinates of vertices
+ if (1){
+ std::map<const BDS_Point*,MVertex*>::const_iterator itx = recoverMap.begin();
+ while (itx != recoverMap.end())
+ {
+ MVertex *v = (itx->second);
+ const BDS_Point *p = (itx->first);
+ v->x() = p->X;
+ v->y() = p->Y;
+ v->z() = p->Z;
+ ++itx;
+ }
+ for (unsigned int i=0;i<gf->mesh_vertices.size();i++)
+ {
+ double u = gf->mesh_vertices[i]->x();
+ double v = gf->mesh_vertices[i]->y();
+ GPoint p = gf->point(u,v);
+ gf->mesh_vertices[i]->x() = p.x();
+ gf->mesh_vertices[i]->y() = p.y();
+ gf->mesh_vertices[i]->z() = p.z();
+ }
+ }
+
+ // fill new gmsh structures with triangles
+ while (1)
+ {
+ if (AllTris.begin() == AllTris.end() ) break;
+ MTri3 *worst = *AllTris.begin();
+ if (worst->isDeleted())
+ {
+ delete worst->tri();
+ }
+ else
+ {
+ gf->triangles.push_back(worst->tri());
+ }
+ delete worst;
+ AllTris.erase(AllTris.begin());
+ }
+}
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
new file mode 100644
index 0000000000000000000000000000000000000000..778489ec07721ebc51606ced190dff72e57c6c1f
--- /dev/null
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -0,0 +1,114 @@
+#ifndef _DELAUNAYINSERTIONFACE_H_
+#define _DELAUNAYINSERTIONFACE_H_
+
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "MElement.h"
+#include <list>
+#include <set>
+#include <map>
+
+class GFace;
+class BDS_Mesh;
+class BDS_Point;
+
+struct gmsh2dMetric
+{
+ double a11,a21,a22;
+ gmsh2dMetric (double lc); // uniform metric
+ gmsh2dMetric (double _a11 = 1, double _a21= 0, double _a22=1)
+ : a11(_a11),a21(_a21),a22(_a22)
+ {}
+ inline double eval ( const double &x, const double &y ) const
+ {
+ return x * a11 * x + 2 * x * a21 * y + y * a22 * y;
+ }
+};
+
+
+class MTri3
+{
+ bool deleted;
+ gmsh2dMetric metric;
+ double xc,yc;
+ double circum_radius;
+ MTriangle *base;
+ MTri3 *neigh[3];
+ public :
+
+ bool isDeleted () const {return deleted;}
+ void forceRadius (double r){circum_radius=r;}
+ double getRadius ()const {return circum_radius;}
+ inline void getCenter (double c[2]) const {c[0]=xc;c[1]=yc;}
+
+ MTri3 ( MTriangle * t, std::vector<gmsh2dMetric> & sizes);
+ inline MTriangle * tri() const {return base;}
+ inline void setNeigh (int iN , MTri3 *n) {neigh[iN]=n;}
+ inline MTri3 *getNeigh (int iN ) const {return neigh[iN];}
+ int inCircumCircle ( const double *p ) const;
+ inline int inCircumCircle ( double x, double y ) const
+ {
+ const double p[2] = {x,y};
+ return inCircumCircle ( p );
+ }
+ inline int inCircumCircle ( const MVertex * v) const
+ {
+ return inCircumCircle ( v->x(), v->y() );
+ }
+
+ void Center_Circum_Aniso(double a, double b, double d, double &x, double &y, double &r) const ;
+
+ double getSurfaceXY () const { return base -> getSurfaceXY() ; };
+ inline void setDeleted (bool d)
+ {
+ deleted = d;
+ }
+ inline bool assertNeigh() const
+ {
+ if (deleted) return true;
+ for (int i=0;i<3;i++)
+ if (neigh[i] && (neigh[i]->isNeigh(this)==false))return false;
+ return true;
+ }
+
+ inline bool isNeigh (const MTri3 *t) const
+ {
+ for (int i=0;i<3;i++)
+ if (neigh[i]==t) return true;
+ return false;
+ }
+
+};
+
+void connectTriangles ( std::list<MTri3*> & );
+void insertVerticesInFace (GFace *gf, BDS_Mesh *);
+
+class compareTri3Ptr
+{
+ public:
+ inline bool operator () ( const MTri3 *a, const MTri3 *b )
+ {
+ if (a->getRadius() > b->getRadius())return true;
+ if (a->getRadius() < b->getRadius())return false;
+ return a<b;
+ }
+};
+
+#endif
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 07d1204a156360b5c15ee2e8dfeaadb240435c3a..aec76e24f8adb9bbce9d317346cf2697ec842a09 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.9 2006-12-01 16:16:50 remacle Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.10 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -319,7 +319,7 @@ void insertVerticesInRegion (GRegion *gr)
gr->tetrahedra.clear();
connectTets ( allTets.begin(), allTets.end() );
- Msg(INFO,"All %d tets were connected",allTets.size());
+ Msg(DEBUG,"All %d tets were connected",allTets.size());
// here the classification should be done
@@ -339,7 +339,7 @@ void insertVerticesInRegion (GRegion *gr)
else
{
if(ITER++%5000 ==0)
- Msg(INFO,"%d points created -- Worst tet radius is %g",vSizes.size(),worst->getRadius());
+ Msg(DEBUG,"%d points created -- Worst tet radius is %g",vSizes.size(),worst->getRadius());
if (worst->getRadius() < 1) break;
double center[3];
worst->tet()->circumcenter(center);
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 2316304da01137324c88bd23e001ff7ac4ec4359..7274d758340e529902b12f19ed20bddf2bdfcf3f 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -1,4 +1,4 @@
-// $Id: Numeric.cpp,v 1.30 2006-11-27 22:22:17 geuzaine Exp $
+// $Id: Numeric.cpp,v 1.31 2007-01-12 13:16:59 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -257,6 +257,30 @@ int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det)
return out;
}
+double det2x2(double mat[2][2])
+{
+ return mat[0][0] *mat[1][1] -mat[1][0] *mat[0][1];
+}
+
+double inv2x2(double mat[2][2], double inv[2][2])
+{
+ const double det = det2x2 ( mat );
+ if(det){
+ double ud = 1. / det;
+ inv[0][0] = ( mat[1][1] ) * ud ;
+ inv[1][0] = -( mat[1][0] ) * ud ;
+ inv[0][1] = -( mat[0][1] ) * ud ;
+ inv[1][1] = ( mat[0][0] ) * ud ;
+ }
+ else{
+ Msg(GERROR, "Singular matrix");
+ for(int i = 0; i < 2; i++)
+ for(int j = 0; j < 2; j++)
+ inv[i][j] = 0.;
+ }
+ return det;
+}
+
double inv3x3(double mat[3][3], double inv[3][3])
{
double det = det3x3(mat);
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 675a0eea5e6783235587ec2d86a002f9d51b9804..2829edb7dabb2c0e388ea66bcdb333c3cd442f16 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -70,10 +70,12 @@ void normal3points(double x0, double y0, double z0,
int sys2x2(double mat[2][2], double b[2], double res[2]);
int sys3x3(double mat[3][3], double b[3], double res[3], double *det);
int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det);
+double det2x2(double mat[2][2]);
double det3x3(double mat[3][3]);
double trace3x3(double mat[3][3]);
double trace2 (double mat[3][3]);
double inv3x3(double mat[3][3], double inv[3][3]);
+double inv2x2(double mat[2][2], double inv[2][2]);
void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
double angle_02pi(double A3);
double angle_plan(double V[3], double P1[3], double P2[3], double n[3]);
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index 65cfffe728984947a7fda7737214087223ac334b..99b5a7a7c81a318ebb3ad151fb658c1bfaf53c51 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $Id: OpenFile.cpp,v 1.137 2007-01-10 13:48:50 geuzaine Exp $
+// $Id: OpenFile.cpp,v 1.138 2007-01-12 13:17:00 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -386,8 +386,8 @@ int MergeFile(char *name, int warn_if_missing)
SBoundingBox3d bb;
bb = GMODEL->bounds();
-
- GMODEL->setMeshSize(0.1 * norm( SVector3 ( bb.max() , bb.min() )));
+ if (GMODEL->getMeshSize() < 0)
+ GMODEL->setMeshSize(0.1 * norm( SVector3 ( bb.max() , bb.min() )));
CTX.geom.draw = 1;
CTX.mesh.changed = ENT_ALL;
diff --git a/benchmarks/2d/Square-01.geo b/benchmarks/2d/Square-01.geo
index 0feef348ca6643e9ef529f2f6aaa9a52dcc41f2c..72c04246dbca630c5e0c642ebd48b6dabad34709 100644
--- a/benchmarks/2d/Square-01.geo
+++ b/benchmarks/2d/Square-01.geo
@@ -1,7 +1,7 @@
/******************************
Square uniformly meshed
******************************/
-lc = .49999;
+lc = .3;
Point(1) = {0.0,0.0,0,lc};
Point(2) = {1,0.0,0,lc};
Point(3) = {1,1,0,lc};
@@ -14,5 +14,5 @@ Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
//Attractor Point{2} = {0.05,0.05,2};
//Mesh.Algorithm = 2;
-Transfinite Line {4,-2} = 130 Using Bump 5;
-Transfinite Line {1,3} = 5 Using Progression 1;
+//Transfinite Line {4,-2} = 130 Using Bump 5;
+//Transfinite Line {1,3} = 5 Using Progression 1;