diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 04d53d6d081b0c48fa45d444edca5e979dbdb6de..f0e3b18e0ec5d337eb563598155b4911d8107d34 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.98 2007-04-20 07:11:26 geuzaine Exp $
+// $Id: CommandLine.cpp,v 1.99 2007-04-26 09:47:37 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -377,6 +377,10 @@ void Get_Options(int argc, char *argv[])
exit(1);
}
}
+ else if(!strcmp(argv[i] + 1, "c1")) {
+ i++;
+ opt_mesh_c1(0, GMSH_SET, 1);
+ }
else if(!strcmp(argv[i] + 1, "optimize_hom")) {
i++;
opt_mesh_smooth_internal_edges(0, GMSH_SET, 1);
diff --git a/Common/Context.h b/Common/Context.h
index 9af74778a8486c524f3e712191a8e3f93719b2fb..c64a09a77fc7d2ca444053de2eaa7c204fd6cfdd 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, smooth_internal_edges;
+ int order, second_order_linear, second_order_incomplete, smooth_internal_edges, c1_continuity;
int min_circ_points;
int constrained_bgmesh, lc_from_curvature;
double normals, tangents, explode;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index baa5a1281561924ef677fe7574af530594675583..91c2121a83ff8f81424e446f102b4f2eb8b5c5b7 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -910,6 +910,8 @@ StringXNumber MeshOptions_Number[] = {
{ F|O, "ElementOrder" , opt_mesh_order , 1. , // "Order" is a reserved token in the parser
"Element order (1=linear elements, N (<6) = elements of higher order)" },
+ { F|O, "C1Continuity" , opt_mesh_c1 , 0. , // "Order" is a reserved token in the parser
+ "Impose C1 continuity to high order meshes, only valid in 2D plane and ElemenOrder = 2 and 3 (todo) (Default : C0 continuity)" },
{ 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 ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 66a81897aa9b311915fffadbd3a2d3e8d501efad..d1094d8214331ecb786192e4c03f81841aee7fe6 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.340 2007-04-21 19:39:59 geuzaine Exp $
+// $Id: Options.cpp,v 1.341 2007-04-26 09:47:38 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -4815,6 +4815,17 @@ double opt_mesh_order(OPT_ARGS_NUM)
return CTX.mesh.order;
}
+double opt_mesh_c1(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET)
+ CTX.mesh.c1_continuity = (int)val;
+#if defined(HAVE_FLTK)
+ if(WID && (action & GMSH_GUI))
+ WID->mesh_butt[21]->value(CTX.mesh.c1_continuity);
+#endif
+ return CTX.mesh.c1_continuity;
+}
+
double opt_mesh_smooth_internal_edges(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 21d6bce11b2497cf38abbcc99ce30aa53f22ad4c..5d06923db737d7553b95e9ef1f366a14f48d38f0 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -469,6 +469,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_c1(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);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 5ee1cbd7b1853c684f9e5b815bfe6fa9c15ed40f..c033d81d88b3e7ba32abfa0b4e4d2f4ab9cc0db4 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.523 2007-04-21 19:39:59 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.524 2007-04-26 09:47:38 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -1090,6 +1090,7 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
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_c1(0, GMSH_SET, WID->mesh_butt[21]->value());
opt_mesh_constrained_bgmesh(0, GMSH_SET, WID->mesh_butt[5]->value());
opt_mesh_points(0, GMSH_SET, WID->mesh_butt[6]->value());
opt_mesh_lines(0, GMSH_SET, WID->mesh_butt[7]->value());
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 463a2ee8b0de95e54aa871d9232a548c66958d50..035c5cda5700e39c2d0069e29e45fdf096082077 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.603 2007-03-29 16:32:23 geuzaine Exp $
+// $Id: GUI.cpp,v 1.604 2007-04-26 09:47:38 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -1619,7 +1619,7 @@ void GUI::reset_external_view_list()
void GUI::create_option_window()
{
int width = 40 * fontsize;
- int height = 13 * BH + 5 * WB;
+ int height = 14 * BH + 5 * WB;
int L = 105 + WB;
if(opt_window) {
@@ -2288,10 +2288,14 @@ void GUI::create_option_window()
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 + 11 * 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 high order elements (8-node quads, etc.)");
mesh_butt[4]->type(FL_TOGGLE_BUTTON);
mesh_butt[4]->callback(mesh_options_ok_cb);
+ mesh_butt[21] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 12 * BH, BW, BH, "Impose C1 continuity (only for degrees 2 and 3 and 2D plane)");
+ mesh_butt[21]->type(FL_TOGGLE_BUTTON);
+ mesh_butt[21]->callback(mesh_options_ok_cb);
+
o->end();
}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index b1a90cda3838ebbdbdb2aabe1ce1af761dd706f4..3324147d221ff5bd6e49ea1588a00ecf30a3278d 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.12 2007-04-12 08:47:25 remacle Exp $
+// $Id: HighOrder.cpp,v 1.13 2007-04-26 09:47:38 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -25,14 +25,107 @@
#include "Message.h"
#include "OS.h"
-// consider 2 points with their tangent
-// x(t) = x(0) * H1(t) + dx(0) * H2(t) + x(1) * H3(t) + dx(1) * H4(t)
+void computeMidPoints(const double *p1,
+ const double *p2,
+ const double *t1,
+ const double *t2,
+ double *p3){
+ // computes the mid point p3 of the second order parametric polynomial curve according to extreme points p1,p2 and tangents t1,t2
+ // The tangent directions given in t1 and t2 are the outgoing normals !!! to the curve p1-p2
+ double x0 = p1[0];
+ double y0 = p1[1];
+ double x1 = p2[0];
+ double y1 = p2[1];
+ double t0x = -t1[0];
+ double t0y = -t1[1];
+ double t1x = t2[0];
+ double t1y = t2[1];
+
+ //computation of the rotation parameters
+ double delta_x = x1-x0;
+ double delta_y = y1-y0;
+ double L = sqrt(delta_x*delta_x+delta_y*delta_y);
+ double theta = atan2(delta_y,delta_x);
+ // double theta = atan(delta_y/delta_x);
+ // if (delta_x<0)
+ // theta += M_PI;
+ double costheta = cos(theta);
+ double sintheta = sin(theta);
+ double theta0 = atan2(t0y,t0x)-theta;
+ double theta1 = atan2(t1y,t1x)-theta;
+ double tan0 = tan(theta0);
+ double tan1 = tan(theta1);
+ delta_x = L;
+ delta_y = 0;
+
+
+ // printf("theta0 = %g theta1 = %g\n",theta0,theta1);
+
+ // computes the 6 coefficients: solving the linear system
+ // x(t) = a1 t^2 * b1 t + c1
+ // y(t) = a2 t^2 * b2 t + c2
+ // for t in [0,1]
+ // with :
+ // x(t), y(t) on pts 0 and 1 (4 conditions)
+ // and dy/dx = dy/dt*dt/dx on those 2 pts (2 conditions)
+
+ if (tan0 == tan1)
+ {
+ p3[0] = 0.5 * (p1[0]+p2[0]);
+ p3[1] = 0.5 * (p1[1]+p2[1]);
+ }
+ else
+ {
+ double b1 = 2*(delta_y-delta_x*tan1)/(tan0-tan1);
+ double a1 = delta_x-b1;
+ double b2 = tan0 * b1;
+ double a2 = delta_y - b2;
+ double c1 = 0;// in general : x0...
+ double c2 = 0;// in general : y0... but axis have changed...
+
+ // back to the initial axis
+ double ax = costheta*a1 - sintheta*a2;
+ double bx = costheta*b1 - sintheta*b2;
+ double cx = x0;//costheta*c1 - sintheta*c2 + x0;
+ double ay = sintheta*a1 + costheta*a2;
+ double by = sintheta*b1 + costheta*b2;
+ double cy = y0;//sintheta*c1 + costheta*c2 + y0;
+
+ // computes the mid point (x,y)
+ double t=0.5;
+ double tsquare=t*t;
+ p3[0] = ax*tsquare + bx*t + cx;
+ p3[1] = ay*tsquare + by*t + cy;
+ }
+ return;
+}
void Hermite2D_C1(SPoint3 &p1, SPoint3 &p2, SPoint3 &t1, SPoint3 &t2,
SPoint3 &one_third, SPoint3 &two_third)
{
- // double L = sqrt((p1.x()-p2.x())*(p1.x()-p2.x()) + (p1.y()-p2.y()) * (p1.y()-p2.y()));
- // SVector3 p1p2 (p2,p1);
+ double L = sqrt((p1.x()-p2.x())*(p1.x()-p2.x()) + (p1.y()-p2.y()) * (p1.y()-p2.y()));
+ SVector3 p1p2 (p2,p1);
+ double theta = atan2(p2.x()-p1.x(),p2.y()-p1.y());
+ double theta1 = atan2(t1.x(),t1.y());
+ double theta2 = atan2(t2.x(),t2.y());
+ double rot[2][2] =
+ {{cos(theta),-sin(theta)},
+ {sin(theta),cos(theta)}};
+
+ double ts[2] = {1./3.,2./3.};
+ SPoint3 *pt[2] = {&one_third,&two_third};
+ for (int i=0 ; i < 2; i++){
+ const double t = ts[i];
+ double H2 = L*t*(t-1)*(t-1);
+ double H4 = -L*(1-t)*t*t;
+ double x3 = t * L;
+ double y3 = H2 * theta1 + H4 * theta2;
+ // pt[i]->x() = p1.x() + x3 * rot[0][0] + y3 * rot[0][1];
+ // pt[i]->y() = p1.y() + x3 * rot[1][0] + y3 * rot[1][1];
+ }
+
+
+
// SVector3 tt1 (t1.x(),t1.y(),0);
// SVector3 tt2 (t2.x(),t2.y(),0);
// const double cost1 = p1p2 * tt1;
@@ -109,7 +202,7 @@ 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)
{
- bool hermite = false;
+ bool c1 = CTX.mesh.c1_continuity;
for(int i = 0; i < ele->getNumEdges(); i++){
MEdge edge = ele->getEdge(i);
@@ -124,33 +217,23 @@ void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
MVertex *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);
double u0 = 0., u1 = 0.;
bool reparamOK = true;
- if(hermite || (!linear && ge->geomType() != GEntity::DiscreteCurve)){
+ if(c1 || (!linear && ge->geomType() != GEntity::DiscreteCurve)){
reparamOK &= reparamOnEdge(v0, ge, u0);
reparamOK &= reparamOnEdge(v1, ge, u1);
}
- if(nPts == 2 && hermite){
+ if(nPts == 1 && c1){
SVector3 tv1 = ge->firstDer(u0);
SVector3 tv2 = ge->firstDer(u1);
- tv1.normalize();
- tv2.normalize();
- SPoint3 t1(tv1.x(), tv1.y(), 0), t2(tv2.x(), tv2.y(), 0);
- SPoint3 vv0(v0->x(), v0->y(), 0), vv1(v1->x(), v1->y(), 0);
- SPoint3 one_third, two_third;
- Hermite2D_C1(vv0, vv1, t1, t2, one_third, two_third);
- MVertex *v1 = new MVertex(one_third.x(), one_third.y(), 0);
- MVertex *v2 = new MVertex(two_third.x(), two_third.y(), 0);
- ge->mesh_vertices.push_back(v1);
- ve.push_back(v1);
- ge->mesh_vertices.push_back(v2);
- ve.push_back(v2);
- if(edge.getVertex(0) == edge.getMinVertex()){
- edgeVertices[p].push_back(v1);
- edgeVertices[p].push_back(v2);
- }
- else{
- edgeVertices[p].push_back(v2);
- edgeVertices[p].push_back(v1);
- }
+ double t1[2] = {-tv1.x(),-tv1.y()};
+ double t2[2] = { tv2.x(), tv2.y()};
+ double p1[2] = { v0->x(), v0->y()};
+ double p2[2] = { v1->x(), v1->y()};
+ double p3[2];
+ computeMidPoints(p1,p2,t1,t2,p3);
+ MVertex *hv = new MVertex(p3[0],p3[1],0);
+ ge->mesh_vertices.push_back(hv);
+ ve.push_back(hv);
+ edgeVertices[p].push_back(hv);
}
else{
std::vector<MVertex*> temp;
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index a52ab04d5c2273a4a5e82c1762528652b8fa3d58..eb4899aff09179597ae49ef26b5e45f0c0a3c101 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.36 2007-04-22 19:41:02 geuzaine Exp $
+// $Id: meshGEdge.cpp,v 1.37 2007-04-26 09:47:38 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -232,7 +232,7 @@ void meshGEdge::operator() (GEdge *ge)
N = ge->meshAttributes.nbPointsTransfinite;
}
else{
- a = Integration(ge, t_begin, t_end, F_Lc, Points, 1.e-7);
+ a = Integration(ge, t_begin, t_end, F_Lc, Points, 1.e-8);
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
}
const double b = a / (double)(N - 1);