diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 39bc1bb08eaef7fc2e25138c1dab211452e93edf..afab1b5a77abdd8cec253f95c4ed2eab946a2165 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.95 2007-03-11 20:18:57 geuzaine Exp $
+// $Id: CommandLine.cpp,v 1.96 2007-03-13 16:11:42 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -86,7 +86,8 @@ void Print_Usage(char *name){
Msg(DIRECT, " -algo string Select mesh algorithm (iso, netgen, tetgen)");
Msg(DIRECT, " -smooth int Set number of mesh smoothing steps");
Msg(DIRECT, " -optimize Optimize quality of tetrahedral elements");
- Msg(DIRECT, " -order int Set mesh order (1, 2)");
+ Msg(DIRECT, " -order int Set mesh order (1, ..., 5)");
+ Msg(DIRECT, " -optimize_hom Optimize higher order meshes (in 2D)");
Msg(DIRECT, " -clscale float Set characteristic length scaling factor");
Msg(DIRECT, " -clcurv Compute characteristic lengths from curvatures");
Msg(DIRECT, " -rand float Set random perturbation factor");
@@ -373,6 +374,10 @@ void Get_Options(int argc, char *argv[])
exit(1);
}
}
+ else if(!strcmp(argv[i] + 1, "optimize_hom")) {
+ i++;
+ opt_mesh_smooth_internal_edges(0, GMSH_SET, 1);
+ }
else if(!strcmp(argv[i] + 1, "format") || !strcmp(argv[i] + 1, "f")) {
i++;
if(argv[i] != NULL) {
diff --git a/Common/Context.h b/Common/Context.h
index 0374e518fb3f8c2b6fdae3306ec491dc7473ed35..c385df687097aa139fa57a3309cb482e5504a073 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -173,7 +173,7 @@ public :
int dual;
int light, light_two_side, light_lines;
int format, nb_smoothing, algo2d, algo3d, algo_recombine;
- int order, second_order_linear, second_order_incomplete;
+ int order, second_order_linear, second_order_incomplete, smooth_internal_edges;
int min_circ_points;
int bgmesh_view_num, constrained_bgmesh, lc_from_curvature;
double normals, tangents, explode;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 2765c041e5e79f98a373b872acba3b97bad3a9b0..84a5d86d0e712bf852a0b3886a4b65323507aa50 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -907,7 +907,9 @@ StringXNumber MeshOptions_Number[] = {
"Display the dual mesh obtained by barycentric subdivision" },
{ F|O, "ElementOrder" , opt_mesh_order , 1. , // "Order" is a reserved token in the parser
- "Element order (1=linear elements, 2=quadratic elements)" },
+ "Element order (1=linear elements, N (<6) = elements of higher order)" },
+ { F|O, "SmoothInternalEdges" , opt_mesh_smooth_internal_edges , 0 , // "Order" is a reserved token in the parser
+ "Number of smoothing steps of internal edges for high order meshes" },
{ F|O, "Explode" , opt_mesh_explode , 1.0 ,
"Element shrinking factor (between 0 and 1)" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index f3e5a26d62662061f9e78cce6ab897329a7192aa..0c6aad7bd07325bd198d9374017ab7beaadbb813 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.336 2007-02-16 08:54:03 geuzaine Exp $
+// $Id: Options.cpp,v 1.337 2007-03-13 16:11:42 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -4804,11 +4804,22 @@ double opt_mesh_order(OPT_ARGS_NUM)
CTX.mesh.order = (int)val;
#if defined(HAVE_FLTK)
if(WID && (action & GMSH_GUI))
- WID->mesh_butt[3]->value(CTX.mesh.order == 2);
+ WID->mesh_value[3]->value(CTX.mesh.order);
#endif
return CTX.mesh.order;
}
+double opt_mesh_smooth_internal_edges(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET)
+ CTX.mesh.smooth_internal_edges = (int)val;
+#if defined(HAVE_FLTK)
+ if(WID && (action & GMSH_GUI))
+ WID->mesh_butt[3]->value(CTX.mesh.smooth_internal_edges);
+#endif
+ return CTX.mesh.smooth_internal_edges;
+}
+
double opt_mesh_second_order_linear(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 8020326e8241b962c9d5db381eb2321906903f9f..6d083a38a92f7d255919d0fdfd2de534c50837ae 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -468,6 +468,7 @@ double opt_mesh_recombine_algo(OPT_ARGS_NUM);
double opt_mesh_min_circ_points(OPT_ARGS_NUM);
double opt_mesh_constrained_bgmesh(OPT_ARGS_NUM);
double opt_mesh_order(OPT_ARGS_NUM);
+double opt_mesh_smooth_internal_edges(OPT_ARGS_NUM);
double opt_mesh_second_order_linear(OPT_ARGS_NUM);
double opt_mesh_second_order_incomplete(OPT_ARGS_NUM);
double opt_mesh_dual(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 0b27d00dd61c2b2efab6fa7fb2b8342752f43854..b02f8eff678f884760f41ecfc7aa636d021d521a 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.517 2007-03-11 20:18:57 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.518 2007-03-13 16:11:42 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -1085,7 +1085,8 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
opt_mesh_reverse_all_normals(0, GMSH_SET, WID->mesh_butt[0]->value());
opt_mesh_lc_from_curvature(0, GMSH_SET, WID->mesh_butt[1]->value());
opt_mesh_optimize(0, GMSH_SET, WID->mesh_butt[2]->value());
- opt_mesh_order(0, GMSH_SET, WID->mesh_butt[3]->value()? 2 : 1);
+ opt_mesh_order(0, GMSH_SET, WID->mesh_value[3]->value());
+ opt_mesh_smooth_internal_edges(0, GMSH_SET, WID->mesh_butt[3]->value());
opt_mesh_second_order_incomplete(0, GMSH_SET, WID->mesh_butt[4]->value());
opt_mesh_constrained_bgmesh(0, GMSH_SET, WID->mesh_butt[5]->value());
opt_mesh_points(0, GMSH_SET, WID->mesh_butt[6]->value());
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index d1c6bf721ce3f5601d7422f3f931f57c285d8a42..984e715e08ea73425505a03388fda6466231eefd 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.600 2007-02-21 08:17:16 geuzaine Exp $
+// $Id: GUI.cpp,v 1.601 2007-03-13 16:11:42 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -2253,33 +2253,40 @@ void GUI::create_option_window()
mesh_value[0]->align(FL_ALIGN_RIGHT);
mesh_value[0]->callback(mesh_options_ok_cb);
- mesh_value[2] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 5 * BH, IW, BH, "Characteristic length factor");
+ mesh_value[3] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 5 * BH, IW, BH, "Element orders (1-5)");
+ mesh_value[3]->minimum(1);
+ mesh_value[3]->maximum(5);
+ mesh_value[3]->step(1);
+ mesh_value[3]->align(FL_ALIGN_RIGHT);
+ mesh_value[3]->callback(mesh_options_ok_cb);
+
+ mesh_value[2] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 6 * BH, IW, BH, "Characteristic length factor");
mesh_value[2]->minimum(0.001);
mesh_value[2]->maximum(1000);
mesh_value[2]->step(0.01);
mesh_value[2]->align(FL_ALIGN_RIGHT);
mesh_value[2]->callback(mesh_options_ok_cb);
- mesh_butt[1] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 6 * BH, BW, BH, "Compute characteritic lenghts automatically from curvatures" );
+ mesh_butt[1] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 7 * BH, BW, BH, "Compute characteritic lenghts automatically from curvatures" );
mesh_butt[1]->type(FL_TOGGLE_BUTTON);
mesh_butt[1]->callback(mesh_options_ok_cb);
- mesh_butt[5] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 7 * BH, BW, BH, "Constrain background mesh with characteristic length field");
+ mesh_butt[5] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 8 * BH, BW, BH, "Constrain background mesh with characteristic length field");
mesh_butt[5]->type(FL_TOGGLE_BUTTON);
mesh_butt[5]->callback(mesh_options_ok_cb);
- mesh_butt[2] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 8 * BH, BW, BH, "Optimize quality of tetrahedra");
+ mesh_butt[2] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 9 * BH, BW, BH, "Optimize quality of tetrahedra");
mesh_butt[2]->type(FL_TOGGLE_BUTTON);
#if !defined(HAVE_NETGEN)
mesh_butt[2]->deactivate();
#endif
mesh_butt[2]->callback(mesh_options_ok_cb);
- mesh_butt[3] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 9 * BH, BW, BH, "Generate second order elements");
+ mesh_butt[3] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 10 * BH, BW, BH, "Optimize high order mesh (currently for 2D - plane)");
mesh_butt[3]->type(FL_TOGGLE_BUTTON);
mesh_butt[3]->callback(mesh_options_ok_cb);
- mesh_butt[4] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 10 * BH, BW, BH, "Use incomplete second order elements (8-node quads, etc.)");
+ mesh_butt[4] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 11 * BH, BW, BH, "Use incomplete second order elements (8-node quads, etc.)");
mesh_butt[4]->type(FL_TOGGLE_BUTTON);
mesh_butt[4]->callback(mesh_options_ok_cb);
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 5545a6ca0fb7db587991ce1d1568d073e2a20e1c..a640b6b6eb3bc5ff3f119d237fca5e1d6eda0b34 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -313,6 +313,12 @@ class MTriangle : public MElement {
bool invertmappingXY(double *p, double *uv, double tol = 1.e-8);
virtual int getNumVertices(){ return 3; }
virtual MVertex *getVertex(int num){ return _v[num]; }
+ virtual MVertex *getOtherVertex(MVertex *v1, MVertex *v2){
+ if( _v[0] != v1 && _v[0] != v2 ) return _v[0];
+ if( _v[1] != v1 && _v[1] != v2 ) return _v[1];
+ if( _v[2] != v1 && _v[2] != v2 ) return _v[2];
+ throw;
+ }
virtual int getNumEdges(){ return 3; }
virtual MEdge getEdge(int num)
{
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 7e5f0c44a957a09d24dda76b589946548323ab58..be39b4f68568a57aa5e9dbd9e6bd4a4a4e218aa6 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshEdge.cpp,v 1.32 2007-03-05 11:07:14 remacle Exp $
+// $Id: gmshEdge.cpp,v 1.33 2007-03-13 16:11:43 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -123,8 +123,6 @@ int gmshEdge::minimumMeshSegments () const
if(geomType() == Circle || geomType() == Ellipse)
return (int)(fabs(c->Circle.t1 - c->Circle.t2) *
(double)CTX.mesh.min_circ_points / Pi) - 1;
- else if(geomType() == Nurb)
- return 5;
else
return GEdge::minimumMeshSegments () ;
}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index a445588a113204a51affff01af8b577eff0731cb..67aea10db72f3ab6f4c0ae2cc4eb9cba9594cd32 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.7 2007-03-02 08:44:55 geuzaine Exp $
+// $Id: HighOrder.cpp,v 1.8 2007-03-13 16:11:43 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -601,6 +601,137 @@ void SetOrder1(GModel *m)
removeHighOrderVertices(*it);
}
+bool straightLine ( std::vector<MVertex *> &l , MVertex *n1, MVertex *n2)
+{
+ // x = a * t + b
+ // x1 = b
+ // x2 = a + b
+ for (int i=0;i<l.size();i++)
+ {
+ MVertex *v = l[i];
+ double b = n1->x();
+ double a = n2->x() - b;
+ double t = (v->x() - b)/a;
+ double by = n1->y();
+ double ay = n2->y() - by;
+ double y = ay * t + by;
+ if (fabs(y-v->y()) > 1.e-11 * CTX.lc)
+ {
+ // printf("coucou %g %g\n",y,v->y());
+ return false;
+ }
+ }
+
+ return true;
+
+}
+
+void smoothInternalEdges ( GFace *gf , edgeContainer &edgeVertices)
+{
+ typedef std::map<std::pair<MVertex*, MVertex*> , std::vector<MElement*> > edge2tris;
+ edge2tris e2t;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ MTriangle *t = gf->triangles[i];
+ for(int j = 0; j < t->getNumEdges(); j++){
+ MEdge edge = t->getEdge(j);
+ std::pair<MVertex*, MVertex*> p(edge.getMinVertex(), edge.getMaxVertex());
+ e2t[p].push_back(t);
+ }
+ }
+ for ( edge2tris::iterator it = e2t.begin() ; it != e2t.end() ; ++it)
+ {
+ std::pair<MVertex*, MVertex*> edge = it->first;
+ std::vector<MVertex*> e1,e2,e3,e4,e;
+ std::vector<MElement*> triangles = it->second;
+
+ if (triangles.size() == 2)
+ {
+ MVertex *n2 = edge.first;
+ MVertex *n4 = edge.second;
+ MTriangle *t1 = (MTriangle*)triangles[0];
+ MTriangle *t2 = (MTriangle*)triangles[1];
+ MVertex *n1 = t1->getOtherVertex (n2,n4);
+ MVertex *n3 = t2->getOtherVertex (n2,n4);
+ if (n1<n2)
+ e1 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n1,n2)];
+ else
+ e1 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2,n1)];
+ if (n2<n3)
+ e2 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2,n3)];
+ else
+ e2 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n3,n2)];
+ if (n3<n4)
+ e3 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n3,n4)];
+ else
+ e3 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4,n3)];
+ if (n4<n1)
+ e4 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4,n1)];
+ else
+ e4 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n1,n4)];
+ if (n2<n4)
+ e = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2,n4)];
+ else
+ e = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4,n2)];
+
+ if ((!straightLine ( e1,n1,n2) ||
+ !straightLine ( e2,n2,n3) ||
+ !straightLine ( e3,n3,n4) ||
+ !straightLine ( e4,n4,n1) ))
+ // && straightLine ( e,n2,n4))
+ {
+ for (int i=0 ; i < e.size() ; i++)
+ {
+ double v = (double)(i+1) / (e.size()+1);
+ double u = 1.-v;
+
+// printf("%g %g\n",u,v);
+
+ MVertex *vert = (n2<n4)?e[i]:e[e.size()-i-1];
+ MVertex *vert1 = (n1<n2)?e1[e1.size()-i-1]:e1[i];
+ MVertex *vert3 = (n3<n4)?e3[i]:e3[e3.size()-i-1];
+ MVertex *vert4 = (n4<n1)?e4[e4.size()-i-1]:e4[i];
+ MVertex *vert2 = (n2<n3)?e2[i]:e2[e2.size()-i-1];
+
+
+// printf("n1 %g %g\n",n1->x(),n1->y());
+// printf("n2 %g %g\n",n2->x(),n2->y());
+// printf("n3 %g %g\n",n3->x(),n3->y());
+// printf("n4 %g %g\n",n4->x(),n4->y());
+
+// printf("vert1 %g %g\n",vert1->x(),vert1->y());
+// printf("vert2 %g %g\n",vert2->x(),vert2->y());
+// printf("vert3 %g %g\n",vert3->x(),vert3->y());
+// printf("vert4 %g %g\n",vert4->x(),vert4->y());
+
+// printf("vert %g %g\n",vert->x(),vert->y());
+
+
+ vert->x() = vert->x() + 0.05 * (
+ (1.-u) * vert4->x() + u * vert2->x() +
+ (1.-v) * vert1->x() + v * vert3->x() -
+ ( (1.-u)*(1.-v) * n1->x()
+ + u * (1.-v) * n2->x()
+ + u*v*n3->x()
+ + (1.-u) * v * n4->x()) - vert->x());
+ vert->y() = vert->y() + 0.05 * (
+ (1.-u) * vert4->y() + u * vert2->y() +
+ (1.-v) * vert1->y() + v * vert3->y() -
+ ( (1.-u)*(1.-v) * n1->y()
+ + u * (1.-v) * n2->y()
+ + u*v*n3->y()
+ + (1.-u) * v * n4->y()) - vert->y());
+
+ // printf("vert %g %g\n\n",vert->x(),vert->y());
+ }
+
+ }
+ }
+ }
+
+
+}
+
+
void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
{
// replace all the elements in the mesh with second order elements
@@ -636,6 +767,13 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
+
+ if (CTX.mesh.smooth_internal_edges)
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
+ {
+ for (int i=0;i<10;i++)smoothInternalEdges ( *it , edgeVertices);
+ }
+
double t2 = Cpu();
Msg(INFO, "Mesh second order complete (%g s)", t2 - t1);
Msg(STATUS1, "Mesh");
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index 43358fce04fd05b4c2090d7337cd935d1d77b610..b754796ced51b6d0dcf60c9806078e29331f4706 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -1,4 +1,4 @@
-lc = 0.001;
+lc = 0.2;
lc2 = 0.2;
lc3 = 0.05;
Point(1) = {1.000000e+00,0.000000e+00,0.000000e+00,lc};
@@ -224,5 +224,5 @@ Line Loop(10) = {2,3,4,1};
Plane Surface(11) = {9,10};
//Attractor Point{101,200} = {0.04,0.04,2};
-Attractor Line{1,2,3,4} = {1,.05,3};
+//Attractor Line{1,2,3,4} = {1,.05,3};
Mesh.Algorithm = 2;
diff --git a/benchmarks/2d/wing-splines.geo b/benchmarks/2d/wing-splines.geo
index 0114a15d643df8a818abb07a2be96c3ed08ed7fd..02492d17efc6730aa446a76043b5014bf2c46a2c 100644
--- a/benchmarks/2d/wing-splines.geo
+++ b/benchmarks/2d/wing-splines.geo
@@ -1,6 +1,6 @@
scale = 1 ;
-lc_wing = 0.01 * scale ;
+lc_wing = 0.04 * scale ;
lc_box = 10 * scale ;
Point(3895) = {1.177410e-02*scale,-2.768003e-03*scale,0,lc_wing};