diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 4b7c68809b496d0bd4b528f1715ee2a82598afd5..d99711f7a3e7f5e54623c3d63b91bc8f2c7470c3 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -45,7 +45,7 @@ void GEdgeCompound::orderEdges()
edges.push_back(_compound[i]);
}
- // find a lonly edge
+ // find a lonely edge
std::map<GVertex*,GEdge*> tempv;
for (std::list<GEdge*>::iterator it = edges.begin() ; it != edges.end() ; ++it){
GVertex *v1 = (*it)->getBeginVertex();
@@ -120,7 +120,7 @@ void GEdgeCompound::orderEdges()
first = last;
last = temp;
_orientation[0] = 0;
- }
+ }
else {
Msg::Error("Compound Edge %d is wrong",tag());
return;
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 4613e00f6328bd5512486096471182fb72f8327d..46a59c6feb8676cebd6d9ffbcf90383078c49c5e 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -159,6 +159,8 @@ class GFace : public GEntity
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const = 0;
// return the curvature computed as the divergence of the normal
+ inline double curvature(const SPoint2 ¶m) const
+ {return curvatureMax(param);}
virtual double curvatureDiv(const SPoint2 ¶m) const;
// return the maximum curvature at a point
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index c3209ba38838e5a8977cb4975c752dfe52d1b6bf..309e6b2e7eb9697aac79e542bec3130aedd3768b 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -526,32 +526,24 @@ void GFaceCompound::computeNormals () const
for ( ; itn != _normals.end() ; ++itn) itn->second.normalize();
}
-double GFaceCompound::curvature(const SPoint2 ¶m) const
+double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
{
parametrize();
double U,V;
GFaceCompoundTriangle *lt;
getTriangle (param.x(),param.y(), <, U,V);
if (!lt){
+ printf("oops\n");
return 0.0;
}
-
- // if (lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
- // SPoint2 pv1, pv2, pv3;
- // bool ok = reparamMeshVertexOnFace(lt->t->getVertex(0), lt->gf, pv1);
- // ok |= reparamMeshVertexOnFace(lt->t->getVertex(1), lt->gf, pv2);
- // ok |= reparamMeshVertexOnFace(lt->t->getVertex(2), lt->gf, pv3);
- // if (ok){
- // SPoint2 pv = pv1*(1.-U-V) + pv2*U + pv3*V;
- // return lt->gf->curvature(pv));
- // }
- // }
-
- // return curvature(lt->t);
+ if (lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
+ SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
+ return lt->gf->curvatureMax(pv);
+ }
return 0.;
}
-double GFaceCompound::curvature(MTriangle *t) const
+double GFaceCompound::curvature(MTriangle *t, double u, double v) const
{
SVector3 n1 = _normals[t->getVertex(0)];
SVector3 n2 = _normals[t->getVertex(1)];
@@ -559,7 +551,7 @@ double GFaceCompound::curvature(MTriangle *t) const
double val[9] = {n1.x(),n2.x(),n3.x(),
n1.y(),n2.y(),n3.y(),
n1.z(),n2.z(),n3.z()};
- // return fabs(t->interpolateDiv (val,0.,0.,0.));
+ return fabs(t->interpolateDiv (val,u,v,0.0));
return 0.;
}
@@ -569,11 +561,12 @@ GPoint GFaceCompound::point(double par1, double par2) const
double U,V;
GFaceCompoundTriangle *lt;
getTriangle (par1, par2, <, U,V);
- SPoint3 p(0, 0, 0);
+ SPoint3 p(3, 3, 0);
if (!lt){
- // Msg::Warning("Re-Parametrized face %d --> point (%g %g) lies outside the domain", tag(), par1,par2);
-
- return GPoint(p.x(),p.y(),p.z(),this);
+ // Msg::Warning("Re-Parametrized face %d --> point (%g %g) lies outside the domain", tag(), par1,par2);
+ GPoint gp (p.x(),p.y(),p.z(),this);
+ gp.setNoSuccess();
+ return gp;
}
if (0 && lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
@@ -626,10 +619,11 @@ static void GFaceCompoundBB(void *a, double*mmin, double*mmax)
const double dx = mmax[0] - mmin[0];
const double dy = mmax[1] - mmin[1];
- mmin[0] -= .02*dx;
- mmin[1] -= .02*dy;
- mmax[0] += .02*dx;
- mmax[1] += .02*dy;
+ const double eps = 0.0;//1.e-12;
+ mmin[0] -= eps*dx;
+ mmin[1] -= eps*dy;
+ mmax[0] += eps*dx;
+ mmax[1] += eps*dy;
}
static void GFaceCompoundCentroid(void *a, double*c)
@@ -669,7 +663,17 @@ void GFaceCompound::getTriangle(double u, double v,
double uv[3] = {u, v, 0};
*lt = (GFaceCompoundTriangle*)Octree_Search(uv, oct);
- if (!(*lt)) return;
+ // if (!(*lt)) {
+// for (int i=0;i<nbT;i++){
+// if (GFaceCompoundInEle (&_gfct[i],uv)){
+// *lt = &_gfct[i];
+// break;
+// }
+// }
+// }
+ if (!(*lt)){
+ return;
+ }
double M[2][2],X[2],R[2];
const SPoint3 p0 = (*lt)->p1;
@@ -686,6 +690,8 @@ void GFaceCompound::getTriangle(double u, double v,
_v = X[1];
}
+
+
void GFaceCompound::buildOct() const
{
SBoundingBox3d bb;
@@ -714,6 +720,40 @@ void GFaceCompound::buildOct() const
it = _compound.begin();
count = 0;
+ for ( ; it != _compound.end() ; ++it){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ MTriangle *t = (*it)->triangles[i];
+ std::map<MVertex*,SPoint3>::const_iterator it0 =
+ coordinates.find(t->getVertex(0));
+ std::map<MVertex*,SPoint3>::const_iterator it1 =
+ coordinates.find(t->getVertex(1));
+ std::map<MVertex*,SPoint3>::const_iterator it2 =
+ coordinates.find(t->getVertex(2));
+ _gfct[count].p1 = it0->second;
+ _gfct[count].p2 = it1->second;
+ _gfct[count].p3 = it2->second;
+ if ((*it)->geomType() != GEntity::DiscreteSurface){
+ reparamMeshVertexOnFace(t->getVertex(0), *it, _gfct[count].gfp1);
+ reparamMeshVertexOnFace(t->getVertex(1), *it, _gfct[count].gfp2);
+ reparamMeshVertexOnFace(t->getVertex(2), *it, _gfct[count].gfp3);
+ }
+ _gfct[count].v1 = SPoint3(t->getVertex(0)->x(),t->getVertex(0)->y(),t->getVertex(0)->z());
+ _gfct[count].v2 = SPoint3(t->getVertex(1)->x(),t->getVertex(1)->y(),t->getVertex(1)->z());
+ _gfct[count].v3 = SPoint3(t->getVertex(2)->x(),t->getVertex(2)->y(),t->getVertex(2)->z());
+ _gfct[count].gf = *it;
+ Octree_Insert(&_gfct[count], oct);
+ count ++;
+ }
+ }
+ nbT = count;
+ Octree_Arrange(oct);
+ printStuff();
+}
+
+void GFaceCompound::printStuff() const
+{
+ std::list<GFace*> :: const_iterator it = _compound.begin();
+
FILE * uvx = fopen("UVX.pos","w");
FILE * uvy = fopen("UVY.pos","w");
FILE * uvz = fopen("UVZ.pos","w");
@@ -737,34 +777,32 @@ void GFaceCompound::buildOct() const
coordinates.find(t->getVertex(1));
std::map<MVertex*,SPoint3>::const_iterator it2 =
coordinates.find(t->getVertex(2));
- _gfct[count].p1 = it0->second;
- _gfct[count].p2 = it1->second;
- _gfct[count].p3 = it2->second;
- if ((*it)->geomType() != GEntity::DiscreteSurface){
- reparamMeshVertexOnFace(t->getVertex(0), *it, _gfct[count].gfp1);
- reparamMeshVertexOnFace(t->getVertex(1), *it, _gfct[count].gfp2);
- reparamMeshVertexOnFace(t->getVertex(2), *it, _gfct[count].gfp3);
- }
- _gfct[count].v1 = SPoint3(t->getVertex(0)->x(),t->getVertex(0)->y(),t->getVertex(0)->z());
- _gfct[count].v2 = SPoint3(t->getVertex(1)->x(),t->getVertex(1)->y(),t->getVertex(1)->z());
- _gfct[count].v3 = SPoint3(t->getVertex(2)->x(),t->getVertex(2)->y(),t->getVertex(2)->z());
- _gfct[count].gf = *it;
- fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(xyzu,"VT(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
- it0->second.x(),it1->second.x(),it2->second.x());
+ (35.*it0->second.x()-t->getVertex(0)->x()), -t->getVertex(0)->y(), (35.*it0->second.y()-t->getVertex(0)->z()),
+ (35.*it1->second.x()-t->getVertex(1)->x()), -t->getVertex(1)->y(), (35.*it1->second.y()-t->getVertex(1)->z()),
+ (35.*it2->second.x()-t->getVertex(2)->x()), -t->getVertex(2)->y(), (35.*it2->second.y()-t->getVertex(2)->z()));
fprintf(xyzv,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
it0->second.y(),it1->second.y(),it2->second.y());
- const double K = fabs(curvature (t));
+ /* fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
+ t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
+ t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
+ it0->second.x(),it1->second.x(),it2->second.x());*/
+ double K1 = curvature(t,it0->second.x(),it0->second.y());
+ double K2 = curvature(t,it1->second.x(),it1->second.y());
+ double K3 = curvature(t,it2->second.x(),it2->second.y());
+ // const double K = fabs(curvature (t));
fprintf(xyzc,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
- K, K, K);
+ K1, K2, K3);
fprintf(uvx,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
it0->second.x(), it0->second.y(), 0.0,
@@ -781,9 +819,6 @@ void GFaceCompound::buildOct() const
it1->second.x(), it1->second.y(), 0.0,
it2->second.x(), it2->second.y(), 0.0,
t->getVertex(0)->z(), t->getVertex(1)->z(), t->getVertex(2)->z());
-
- Octree_Insert(&_gfct[count], oct);
- count ++;
}
}
fprintf(uvx,"};\n");
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index e1587a3de0451acdb084d63e8ba58450af837b40..f4b126d17d09d1aaa1dcea6dc4e8fccb45d57fa4 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -47,6 +47,7 @@ class GFaceCompound : public GFace {
protected:
std::list<GFace*> _compound;
std::list<GEdge*> _U0, _U1, _V0, _V1;
+ mutable int nbT;
mutable GFaceCompoundTriangle *_gfct;
mutable Octree *oct;
mutable std::map<MVertex*,SPoint3> coordinates;
@@ -58,7 +59,8 @@ class GFaceCompound : public GFace {
void getBoundingEdges();
void getTriangle(double u, double v, GFaceCompoundTriangle **lt,
double &_u, double &_v) const;
- virtual double curvature(MTriangle *t) const;
+ virtual double curvature(MTriangle *t, double u, double v) const;
+ void printStuff() const;
public:
typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag,
@@ -79,7 +81,7 @@ public:
virtual SPoint2 getCoordinates(MVertex *v) const;
virtual bool buildRepresentationCross(){ return false; }
- virtual double curvature(const SPoint2 ¶m) const;
+ virtual double curvatureMax(const SPoint2 ¶m) const;
private:
typeOfIsomorphism _type;
};
diff --git a/Geo/GPoint.h b/Geo/GPoint.h
index 4248b457360b9fab5c3c643ad2cc3d0c891bdb47..e0d06f8cc347c0eac67d37cd71039c936eb08f09 100644
--- a/Geo/GPoint.h
+++ b/Geo/GPoint.h
@@ -16,6 +16,7 @@ class GPoint
double X, Y, Z;
const GEntity *e;
double par[2];
+ bool success;
public:
inline double x() const { return X; }
inline double y() const { return Y; }
@@ -27,16 +28,16 @@ class GPoint
inline double v() const { return par[1]; }
inline const GEntity* g() const { return e; }
GPoint (double _x=0, double _y=0, double _z=0, const GEntity *onwhat=0)
- : X(_x), Y(_y), Z(_z), e(onwhat)
+ : X(_x), Y(_y), Z(_z), e(onwhat) ,success(true)
{
}
GPoint (double _x, double _y, double _z, const GEntity *onwhat, double p)
- : X(_x), Y(_y), Z(_z), e(onwhat)
+ : X(_x), Y(_y), Z(_z), e(onwhat),success(true)
{
par[0] = p;
}
GPoint (double _x, double _y, double _z, const GEntity *onwhat, double p[2])
- : X(_x), Y(_y), Z(_z), e(onwhat)
+ : X(_x), Y(_y), Z(_z), e(onwhat),success(true)
{
par[0] = p[0];
par[1] = p[1];
@@ -48,6 +49,12 @@ class GPoint
double dz = Z - p.z();
return sqrt(dx * dx + dy * dy + dz * dz);
}
+ bool succeeded () const{
+ return success;
+ }
+ bool setNoSuccess (){
+ success = false;
+ }
};
#endif
diff --git a/Geo/Makefile b/Geo/Makefile
index 2e30b5be569cfd51a79e8efba5df10d61d30e5fa..f3b4d07faf6ab82efa9871d02d618ffb8452168d 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -15,7 +15,7 @@ INC = ${DASH}I../Common ${DASH}I../Geo ${DASH}I../Mesh\
CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE}
-SRC = GEntity.cpp\
+SRC = GEntity.cpp STensor3.cpp\
GVertex.cpp GEdge.cpp GFace.cpp GRegion.cpp\
GEdgeLoop.cpp GEdgeCompound.cpp GFaceCompound.cpp GRegionCompound.cpp\
gmshVertex.cpp gmshEdge.cpp gmshFace.cpp gmshRegion.cpp gmshSurface.cpp\
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 5183002a463138f0dd4b88073c5f70e058b5cbee..efc401792ca2c428dc5933844e78e30baf0d5dca 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1196,7 +1196,11 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
LC /= (sTot);
GPoint gp = gf->point(U * scalingU, V * scalingV);
-
+
+ if (!gp.succeeded()){
+ printf ("iha\n");
+ return false;
+ }
const double oldX = p->X;
const double oldY = p->Y;
const double oldZ = p->Z;
@@ -1307,6 +1311,8 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point *p, GFace *gf)
}
GPoint gp = gf->point(U * scalingU, V * scalingV);
+ if (!gp.succeeded())return false;
+
p->u = U;
p->v = V;
p->lc() = LC;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 6e68114106d52f3247775c1f7a18b9f29a1e2bd4..061f171330302aed8e94d08426c68afb4de24b38 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -40,9 +40,11 @@ static double max_surf_curvature(const GEdge *ge, double u)
std::list<GFace *> faces = ge->faces();
std::list<GFace *>::iterator it = faces.begin();
while(it != faces.end()){
- SPoint2 par = ge->reparamOnFace((*it), u, 1);
- double cc = (*it)->curvatureDiv(par);
- val = std::max(cc, val);
+ if ((*it)->geomType() != GEntity::CompoundSurface){
+ SPoint2 par = ge->reparamOnFace((*it), u, 1);
+ double cc = (*it)->curvature(par);
+ val = std::max(cc, val);
+ }
++it;
}
return val;
@@ -90,7 +92,7 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
case 2:
{
GFace *gf = (GFace *)ge;
- Crv = gf->curvatureDiv(SPoint2(U, V));
+ Crv = gf->curvature(SPoint2(U, V));
}
break;
}
@@ -150,7 +152,6 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
double l3 = MAX_LC;
if(CTX::instance()->mesh.lcFromCurvature && ge->dim() < 3)
l3 = LC_MVertex_CURV(ge, U, V);
-
// lc from fields
double l4 = MAX_LC;
FieldManager *fields = GModel::current()->getFields();
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index ef86a5d36e0e6d8cddc613f74f1e06dd862144c6..bba90346b896023330797649b6a0d01a7b143747 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1010,7 +1010,7 @@ extern double mesh_functional_distorsion(MTriangle *t, double u, double v);
static void printJacobians(GModel *m, const char *nm)
{
- const int n = 15;
+ const int n = 100;
double D[n][n], X[n][n], Y[n][n], Z[n][n];
FILE *f = fopen(nm,"w");
@@ -1025,9 +1025,12 @@ static void printJacobians(GModel *m, const char *nm)
double v = (double)k / (n - 1);
t->pnt(u, v, 0, pt);
D[i][k] = mesh_functional_distorsion(t, u, v);
- X[i][k] = pt.x();
- Y[i][k] = pt.y();
- Z[i][k] = pt.z();
+ //X[i][k] = u;
+ //Y[i][k] = v;
+ //Z[i][k] = 0.0;
+ X[i][k] = pt.x();
+ Y[i][k] = pt.y();
+ Z[i][k] = pt.z();
}
}
for(int i= 0; i < n -1; i++){
@@ -1128,7 +1131,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
if(displ2D) delete displ2D;
if(displ3D) delete displ3D;
- //printJacobians(m, "smoothness.pos");
+ // printJacobians(m, "smoothness.pos");
double t2 = Cpu();
Msg::Info("Meshing order %d complete (%g s)", order, t2 - t1);
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index 4ec085e0db7837c07effa0857297ef96ba223076..eba8d32bffee66c13155e30d2fb4926505c130c3 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -511,11 +511,11 @@ void gmshHighOrderSmoother::smooth_metric ( std::vector<MElement*> & all, GFace
double dx0 = smooth_metric_ ( v, gf, myAssembler, verticesToMove,El);
double dx = dx0;
- // printf(" dx0 = %12.5E\n",dx0);
+ printf(" dx0 = %12.5E\n",dx0);
int iter = 0;
while(1){
double dx2 = smooth_metric_ ( v,gf, myAssembler, verticesToMove,El);
- // printf(" dx2 = %12.5E\n",dx2);
+ printf(" dx2 = %12.5E\n",dx2);
if (fabs(dx2-dx) < 1.e-4 * dx0)break;
if (iter++ > 10)break;
dx = dx2;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index e282b100240deb71b0ee2296dc090eb25d68a79a..fd9b358a954450e9aeefee1a092f30661dd91518 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1370,7 +1370,7 @@ void deMeshGFace::operator() (GFace *gf)
gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
}
-const int debugSurface = -1;
+const int debugSurface = -1;//-100;
void meshGFace::operator() (GFace *gf)
{
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 491ef1912f50bcc9dd20c91cc3c3e80950a3d058..b9d0a36bb67eb67e19870185204f4fbbbfae25fd 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -38,9 +38,13 @@ inline double computeEdgeLinearLength(BDS_Point *p1,
double SCALINGU,
double SCALINGV)
{
+
GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
0.5 * (p1->v + p2->v) * SCALINGV));
+ if (!GP.succeeded())
+ return computeEdgeLinearLength(p1,p2);
+
const double dx1 = p1->X - GP.x();
const double dy1 = p1->Y - GP.y();
const double dz1 = p1->Z - GP.z();
@@ -410,17 +414,22 @@ void splitEdgePassUnsorted(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
const double coord = computeEdgeMiddleCoord((*it)->p1, (*it)->p2, gf,
m.scalingU, m.scalingV);
BDS_Point *mid;
- mid = m.add_point(++m.MAXPOINTNUMBER,
- coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
- coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v, gf);
- mid->lcBGM() = BGM_MeshSize
- (gf,
- (coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u)*m.scalingU,
- (coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
- mid->X,mid->Y,mid->Z);
- mid->lc() = 0.5 * ((*it)->p1->lc() + (*it)->p2->lc());
- if(!m.split_edge(*it, mid)) m.del_point(mid);
- else nb_split++;
+
+ GPoint gpp = gf->point(coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
+ coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v);
+ if (gpp.succeeded()){
+ mid = m.add_point(++m.MAXPOINTNUMBER,
+ coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
+ coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v, gf);
+ mid->lcBGM() = BGM_MeshSize
+ (gf,
+ (coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u)*m.scalingU,
+ (coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
+ mid->X,mid->Y,mid->Z);
+ mid->lc() = 0.5 * ((*it)->p1->lc() + (*it)->p2->lc());
+ if(!m.split_edge(*it, mid)) m.del_point(mid);
+ else nb_split++;
+ }
}
}
++it;
@@ -451,15 +460,21 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
BDS_Point *mid ;
double U = coord * e->p1->u + (1 - coord) * e->p2->u;
double V = coord * e->p1->v + (1 - coord) * e->p2->v;
- mid = m.add_point(++m.MAXPOINTNUMBER, U , V, gf);
- mid->lcBGM() = BGM_MeshSize
- (gf,
- (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());
- if(!m.split_edge(e, mid)) m.del_point(mid);
- else nb_split++;
+
+ GPoint gpp = gf->point(U,V);
+ if (gpp.succeeded()){
+ mid = m.add_point(++m.MAXPOINTNUMBER, gpp.x(),gpp.y(),gpp.z());
+ mid->u = U;
+ mid->v = V;
+ mid->lcBGM() = BGM_MeshSize
+ (gf,
+ (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());
+ if(!m.split_edge(e, mid)) m.del_point(mid);
+ else nb_split++;
+ }
}
}
}
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 8e65c3f09eaaac176b425311e48bd19f992dd5fb..8de0c48b2a06a58538d795d869b8bcf02e9e8dd6 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -18,6 +18,8 @@
const double LIMIT_ = 0.5 * sqrt(2.0);
+static const bool _experimental_anisotropic_blues_band_ = false;
+
// This function controls the frontal algorithm
static bool isActive(MTri3 *t, double limit_, int &active)
{
@@ -30,6 +32,7 @@ static bool isActive(MTri3 *t, double limit_, int &active)
return false;
}
+
bool circumCenterMetricInTriangle(MTriangle *base,
const double *metric,
const std::vector<double> &Us,
@@ -74,6 +77,43 @@ void circumCenterMetric(double *pa, double *pb, double *pc,
2. * (x[0] - pa[0]) * (x[1] - pa[1]) * b;
}
+void circumCenterMetricXYZ(double *p1, double *p2, double *p3, SMetric3 & metric,
+ double *res, double *uv, double &radius)
+{
+ double v1[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
+ double v2[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
+ double vx[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
+ double vy[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
+ double vz[3]; prodve(vx, vy, vz); prodve(vz, vx, vy);
+ norme(vx); norme(vy); norme(vz);
+ double p1P[2] = {0.0, 0.0};
+ double p2P[2]; prosca(v1, vx, &p2P[0]); prosca(v1, vy, &p2P[1]);
+ double p3P[2]; prosca(v2, vx, &p3P[0]); prosca(v2, vy, &p3P[1]);
+ double resP[2];
+
+ gmshMatrix<double> T(3,3);
+ for (int i=0;i<3;i++)T(0,i) = vx[i];
+ for (int i=0;i<3;i++)T(1,i) = vy[i];
+ for (int i=0;i<3;i++)T(2,i) = vz[i];
+ SMetric3 tra = metric.transform(T);
+ double mm[3] = {tra(0,0),tra(0,1),tra(1,1)};
+
+ circumCenterMetric(p1P, p2P, p3P, mm, resP, radius);
+
+ if(uv){
+ double mat[2][2] = {{p2P[0] - p1P[0], p3P[0] - p1P[0]},
+ {p2P[1] - p1P[1], p3P[1] - p1P[1]}};
+ double rhs[2] = {resP[0] - p1P[0], resP[1] - p1P[1]};
+ sys2x2(mat, rhs, uv);
+ }
+
+ res[0] = p1[0] + resP[0] * vx[0] + resP[1] * vy[0];
+ res[1] = p1[1] + resP[0] * vx[1] + resP[1] * vy[1];
+ res[2] = p1[2] + resP[0] * vx[2] + resP[1] * vy[2];
+}
+
+
+
void circumCenterMetric(MTriangle *base,
const double *metric,
const std::vector<double> &Us,
@@ -98,6 +138,39 @@ void buildMetric(GFace *gf, double *uv, double *metric)
metric[2] = dot(der.second(), der.second());
}
+// m 3x3
+// d 3x2
+// M = d^T m d
+
+void buildMetric(GFace *gf, double *uv, SMetric3 & m, double *metric)
+{
+ Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(uv[0], uv[1]));
+
+ SVector3 x1 (m(0,0) * der.first().x() +
+ m(1,0) * der.first().y() +
+ m(2,0) * der.first().z(),
+ m(0,1) * der.first().x() +
+ m(1,1) * der.first().y() +
+ m(2,1) * der.first().z(),
+ m(0,2) * der.first().x() +
+ m(1,2) * der.first().y() +
+ m(2,2) * der.first().z());
+ SVector3 x2 (m(0,0) * der.second().x() +
+ m(1,0) * der.second().y() +
+ m(2,0) * der.second().z(),
+ m(0,1) * der.second().x() +
+ m(1,1) * der.second().y() +
+ m(2,1) * der.second().z(),
+ m(0,2) * der.second().x() +
+ m(1,2) * der.second().y() +
+ m(2,2) * der.second().z());
+
+ metric[0] = dot(x1, der.first());
+ metric[1] = dot(x2, der.first());
+ metric[2] = dot(x2, der.second());
+}
+
+
int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3,
double *uv, double *metric)
{
@@ -137,22 +210,28 @@ int inCircumCircleAniso(GFace *gf, MTriangle *base,
return d2 < Radius2;
}
-MTri3::MTri3(MTriangle *t, double lc) : deleted(false), base(t)//,done(0)
+MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric) : deleted(false), base(t)//,done(0)
{
neigh[0] = neigh[1] = neigh[2] = 0;
-
- // compute the circumradius of the triangle
+ double center[3];
double pa[3] = {base->getVertex(0)->x(), base->getVertex(0)->y(), base->getVertex(0)->z()};
double pb[3] = {base->getVertex(1)->x(), base->getVertex(1)->y(), base->getVertex(1)->z()};
double pc[3] = {base->getVertex(2)->x(), base->getVertex(2)->y(), base->getVertex(2)->z()};
- double center[3];
- circumCenterXYZ(pa, pb, pc, center);
- const double dx = base->getVertex(0)->x() - center[0];
- const double dy = base->getVertex(0)->y() - center[1];
- const double dz = base->getVertex(0)->z() - center[2];
+
+ // compute the circumradius of the triangle
- circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
- circum_radius /= lc;
+ if (!metric){
+ circumCenterXYZ(pa, pb, pc, center);
+ const double dx = base->getVertex(0)->x() - center[0];
+ const double dy = base->getVertex(0)->y() - center[1];
+ const double dz = base->getVertex(0)->z() - center[2];
+ circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
+ circum_radius /= lc;
+ }
+ else{
+ double uv[2];
+ circumCenterMetricXYZ(pa, pb, pc, *metric, center, uv, circum_radius);
+ }
}
int MTri3::inCircumCircle(const double *p) const
@@ -369,24 +448,37 @@ bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
bool onePointIsTooClose = false;
while (it != shell.end()){
MTriangle *t = new MTriangle(it->v[0], it->v[1], v);
- double lc = 0.3333333333 * (vSizes[t->getVertex(0)->getNum()] +
- vSizes[t->getVertex(1)->getNum()] +
- vSizes[t->getVertex(2)->getNum()]);
- double lcBGM = 0.3333333333 * (vSizesBGM[t->getVertex(0)->getNum()] +
- vSizesBGM[t->getVertex(1)->getNum()] +
- vSizesBGM[t->getVertex(2)->getNum()]);
+ const double ONE_THIRD = 1./3.;
+ double lc = ONE_THIRD * (vSizes[t->getVertex(0)->getNum()] +
+ vSizes[t->getVertex(1)->getNum()] +
+ vSizes[t->getVertex(2)->getNum()]);
+ double lcBGM = ONE_THIRD * (vSizesBGM[t->getVertex(0)->getNum()] +
+ vSizesBGM[t->getVertex(1)->getNum()] +
+ vSizesBGM[t->getVertex(2)->getNum()]);
double LL = Extend1dMeshIn2dSurfaces() ? std::min(lc, lcBGM) : lcBGM;
- MTri3 *t4 = new MTri3(t, LL);
- double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
+ MTri3 *t4;
+ if (_experimental_anisotropic_blues_band_){
+ SMetric3 metrBGM = interpolation ( vMetricsBGM[t->getVertex(0)->getNum()],
+ vMetricsBGM[t->getVertex(1)->getNum()],
+ vMetricsBGM[t->getVertex(2)->getNum()],
+ ONE_THIRD,ONE_THIRD);
+
+ t4 = new MTri3(t, LL,&metrBGM);
+ }
+ else{
+ t4 = new MTri3(t, LL);
+ }
+
+ double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
(it->v[0]->y() - v->y()) * (it->v[0]->y() - v->y()) +
(it->v[0]->z() - v->z()) * (it->v[0]->z() - v->z()));
- double d2 = sqrt((it->v[1]->x() - v->x()) * (it->v[1]->x() - v->x()) +
+ double d2 = sqrt((it->v[1]->x() - v->x()) * (it->v[1]->x() - v->x()) +
(it->v[1]->y() - v->y()) * (it->v[1]->y() - v->y()) +
(it->v[1]->z() - v->z()) * (it->v[1]->z() - v->z()));
- if (d1 < LL * .25 || d2 < LL * .25) onePointIsTooClose = true;
-
- // if (t4->getRadius () < LIMIT_ / 2) onePointIsTooClose = true;
+ if (d1 < LL * .25 || d2 < LL * .25) onePointIsTooClose = true;
+
+ // if (t4->getRadius () < LIMIT_ / 2) onePointIsTooClose = true;
newTris[k++] = t4;
// all new triangles are pushed front in order to be able to
@@ -521,25 +613,21 @@ static void insertAPoint(GFace *gf,
GPoint p = gf->point(center[0], center[1]);
// 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]);
-
- // printf("%g %g -> %g %g %g\n",center[0], center[1],v->x(), v->y(), v->z());
-
v->setNum(Us.size());
double lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getNum()] +
uv[0] * vSizes [ptin->tri()->getVertex(1)->getNum()] +
uv[1] * vSizes [ptin->tri()->getVertex(2)->getNum()]);
- // double eigMetricSurface = gf->getMetricEigenvalue(SPoint2(center[0],center[1]));
double 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);
+ // vMetricsBGM.push_back(metr);
vSizesBGM.push_back(lc);
vSizes.push_back(lc1);
Us.push_back(center[0]);
Vs.push_back(center[1]);
if (!insertVertex(gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
- Us, Vs, metric)) {
- Msg::Debug("2D Delaunay : a cavity is not star shaped");
+ Us, Vs, metric) || !p.succeeded()) {
+ // Msg::Debug("2D Delaunay : a cavity is not star shaped");
AllTris.erase(it);
worst->forceRadius(-1);
AllTris.insert(worst);
@@ -550,7 +638,7 @@ static void insertAPoint(GFace *gf,
}
}
else {
- Msg::Debug("Point %g %g is outside (%g %g , %g %g , %g %g) (metric %g %g %g)",
+ /* Msg::Debug("Point %g %g is outside (%g %g , %g %g , %g %g) (metric %g %g %g)",
center[0], center[1],
Us[base->getVertex(0)->getNum()],
Vs[base->getVertex(0)->getNum()],
@@ -558,7 +646,7 @@ static void insertAPoint(GFace *gf,
Vs[base->getVertex(1)->getNum()],
Us[base->getVertex(2)->getNum()],
Vs[base->getVertex(2)->getNum()],
- metric[0], metric[1], metric[2]);
+ metric[0], metric[1], metric[2]);*/
AllTris.erase(it);
worst->forceRadius(0);
AllTris.insert(worst);
@@ -571,7 +659,7 @@ void gmshBowyerWatson(GFace *gf)
std::vector<double> vSizes, vSizesBGM, Us, Vs;
std::vector<SMetric3> vMetricsBGM;
- buildMeshGenerationDataStructures(gf, AllTris, vSizes, vSizesBGM, Us, Vs);
+ buildMeshGenerationDataStructures(gf, AllTris, vSizes, vSizesBGM, vMetricsBGM,Us, Vs);
// _printTris ("before.pos", AllTris, Us,Vs);
int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
@@ -600,7 +688,12 @@ void gmshBowyerWatson(GFace *gf)
(Vs[base->getVertex(0)->getNum()] +
Vs[base->getVertex(1)->getNum()] +
Vs[base->getVertex(2)->getNum()]) / 3.};
- buildMetric(gf, pa, metric);
+ SMetric3 m = interpolation (vMetricsBGM[base->getVertex(0)->getNum()],
+ vMetricsBGM[base->getVertex(1)->getNum()],
+ vMetricsBGM[base->getVertex(2)->getNum()],
+ pa[0],pa[1]);
+ buildMetric(gf, pa, m, metric);
+ //buildMetric(gf, pa, metric);
circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
insertAPoint(gf, AllTris.begin(), center, metric, Us, Vs, vSizes, vSizesBGM, vMetricsBGM,
AllTris);
@@ -677,7 +770,7 @@ void gmshBowyerWatsonFrontal(GFace *gf)
//testTensor();
- buildMeshGenerationDataStructures(gf, AllTris, vSizes, vSizesBGM, Us, Vs);
+ buildMeshGenerationDataStructures(gf, AllTris, vSizes, vSizesBGM, vMetricsBGM,Us, Vs);
// delaunise the initial mesh
int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 1a86bb49744f2f2b9a96636b31bb2b8ce56f0cd2..f3f748b5f449ad4c4c20250ab59db30488c52bc7 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -7,6 +7,7 @@
#define _MESH_GFACE_DELAUNAY_INSERTIONFACE_H_
#include "MTriangle.h"
+#include "STensor3.h"
#include <list>
#include <set>
#include <map>
@@ -50,7 +51,7 @@ class MTri3
void forceRadius (double r){ circum_radius = r; }
double getRadius () const { return circum_radius; }
- MTri3(MTriangle *t, double lc);
+ MTri3(MTriangle *t, double lc, SMetric3 *m = 0);
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]; }
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 0f0b657d58e7943d277e891342eb906f57a10549..ad0f32730cf9d9910c62eac8f56b890265baf1ac 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -50,6 +50,7 @@ void buildMeshGenerationDataStructures(GFace *gf,
std::set<MTri3*, compareTri3Ptr> &AllTris,
std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
+ std::vector<SMetric3> &vMetricsBGM,
std::vector<double> &Us,
std::vector<double> &Vs)
{
@@ -79,6 +80,7 @@ void buildMeshGenerationDataStructures(GFace *gf,
it->first->setNum(NUM++);
vSizes.push_back(it->second);
vSizesBGM.push_back(it->second);
+ vMetricsBGM.push_back(SMetric3(it->second));
SPoint2 param;
reparamMeshVertexOnFace(it->first, gf, param);
Us.push_back(param[0]);
@@ -235,9 +237,11 @@ void laplaceSmoothing(GFace *gf)
ver->setParameter(0, cu / fact);
ver->setParameter(1, cv / fact);
GPoint pt = gf->point(SPoint2(cu / fact, cv / fact));
- ver->x() = pt.x();
- ver->y() = pt.y();
- ver->z() = pt.z();
+ if (pt.succeeded()){
+ ver->x() = pt.x();
+ ver->y() = pt.y();
+ ver->z() = pt.z();
+ }
}
}
++it;
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index 729740391d7a3c260ff27da7cb8dd0f882ab39c9..67f7fa2f00a46565057e79789e9d0ba527926594 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -11,6 +11,7 @@
#include "MElement.h"
#include "MEdge.h"
#include "meshGFaceDelaunayInsertion.h"
+#include "STensor3.h"
class GFace;
class MVertex;
@@ -50,6 +51,7 @@ void buildMeshGenerationDataStructures(GFace *gf,
std::set<MTri3*, compareTri3Ptr> &AllTris,
std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
+ std::vector<SMetric3> &vMetricsBGM,
std::vector<double> &Us,
std::vector<double> &Vs);
void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 92fb0015baf81c9eaf22db4aa838f76129c3aff5..432e7cd86d62f842366c457d286e42484922f341 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -202,6 +202,7 @@ double mesh_functional_distorsion(MTriangle *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);
diff --git a/Numeric/GmshMatrix.cpp b/Numeric/GmshMatrix.cpp
index 07a4e792e3d254e7e0021583dcb2a7a3e3bf33fa..f8a21ad5ab50dbc9a40ff640198b7d416bf68d67 100644
--- a/Numeric/GmshMatrix.cpp
+++ b/Numeric/GmshMatrix.cpp
@@ -132,7 +132,7 @@ bool gmshMatrix<double>::invertInPlace()
if(info == 0) return true;
if(info > 0)
- Msg::Error("U(%d,%d)=0 in matric inversion", info, info);
+ Msg::Error("U(%d,%d)=0 in matrix inversion", info, info);
else
Msg::Error("Wrong %d-th argument in matrix inversion", -info);
return false;
diff --git a/Numeric/gmshElasticity.cpp b/Numeric/gmshElasticity.cpp
index 3d323ca272303b28d002025cf4b6d45afff410ee..da5a790994e7623ecddad123617bfa65e0e0fa43 100644
--- a/Numeric/gmshElasticity.cpp
+++ b/Numeric/gmshElasticity.cpp
@@ -71,6 +71,7 @@ void gmshElasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
}
BTH.set_all(0.);
BTH.gemm(BT, H);
+ // printf("detJ = %22.15E\n",detJ);
m.gemm(BTH, B, weight * detJ, 1.);
}
}
diff --git a/benchmarks/2d/circle2_GEO.geo b/benchmarks/2d/circle2_GEO.geo
index 4922514fe0ca3506df3070bc14e0109ce0002128..2d493a82fd0275143b127b2bc084df7814ffd4d2 100644
--- a/benchmarks/2d/circle2_GEO.geo
+++ b/benchmarks/2d/circle2_GEO.geo
@@ -6,4 +6,4 @@ CreateTopology;
Compound Line (100) = {1,2,3,4};
Compound Line (200) = {11,22,33,44};
-//Compound Surface(1000)={47} Boundary {{1,2,3,4}, {11,22,33,44}};
+Compound Surface(1000)={47} Boundary {{1,2,3,4}, {11,22,33,44}};
diff --git a/benchmarks/2d/square.geo b/benchmarks/2d/square.geo
index 6fa7083a5cafea225eb732e7eb90f59156a74f1b..ec29fb7bdc8f39fcfa8fd9ff2fd4d09415cde094 100644
--- a/benchmarks/2d/square.geo
+++ b/benchmarks/2d/square.geo
@@ -11,7 +11,3 @@ Line(4) = {1, 2};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(10) = {5};
-//Compound Line(150)={1,2};
-//Compound Line(160)={3,4};
-
-//Compound Surface(170)={10} Boundary {{}};
diff --git a/benchmarks/2d/square_CLASS_GEO.geo b/benchmarks/2d/square_CLASS_GEO.geo
index c584ddcd10174ea338150ec30000144211832c7d..d68ebf5567ffcb2f778dea1f347c12e186de9e79 100644
--- a/benchmarks/2d/square_CLASS_GEO.geo
+++ b/benchmarks/2d/square_CLASS_GEO.geo
@@ -1,6 +1,6 @@
Mesh.CharacteristicLengthFactor=6;
-Merge "square_CLASS.msh";
+Merge "square.msh";
CreateTopology;
Compound Line(150)={1,2};