*** empty log message ***

Common/CommandLine.cpp

-// $Id: CommandLine.cpp,v 1.125 2008-03-24 21:39:46 geuzaine Exp $
+// $Id: CommandLine.cpp,v 1.126 2008-03-26 09:37:49 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -82,7 +82,7 @@ void Print_Usage(const char *name)
   Msg(DIRECT, "  -format string        Set output mesh format (msh, msh1, msh2, unv, vrml, stl, mesh,");
   Msg(DIRECT, "                          bdf, p3d, cgns, med)");
   Msg(DIRECT, "  -bin                  Use binary format when available");  
-  Msg(DIRECT, "  -algo string          Select mesh algorithm (iso, netgen, tetgen)");
+  Msg(DIRECT, "  -algo string          Select mesh algorithm (de, del2d, frontal, iso, netgen, tetgen)");
   Msg(DIRECT, "  -smooth int           Set number of mesh smoothing steps");
   Msg(DIRECT, "  -optimize[_netgen]    Optimize quality of tetrahedral elements");
   Msg(DIRECT, "  -order int            Set mesh order (1, ..., 5)");
@@ -527,6 +527,8 @@ void Get_Options(int argc, char *argv[])
             CTX.mesh.algo2d = ALGO_3D_TETGEN_DELAUNAY;
           else if(!strncmp(argv[i], "netgen", 6))
             CTX.mesh.algo3d = ALGO_3D_NETGEN;
+          else if(!strncmp(argv[i], "frontal", 7))
+            CTX.mesh.algo2d = ALGO_2D_FRONTAL;
           else if(!strncmp(argv[i], "del2d", 5) || !strncmp(argv[i], "tri", 3))
             CTX.mesh.algo2d = ALGO_2D_DELAUNAY;
           else if(!strncmp(argv[i], "bds", 3))

Common/GmshDefines.h

@@ -94,6 +94,7 @@
 #define ALGO_2D_TRIANGLE            3 // unused
 #define ALGO_2D_MESHADAPT           4
 #define ALGO_2D_DELAUNAY            5
+#define ALGO_2D_FRONTAL             6
 
 // 3D mesh algorithms
 #define ALGO_3D_TETGEN_DELAUNAY    1

Mesh/meshGFace.cpp

-// $Id: meshGFace.cpp,v 1.128 2008-03-25 20:25:35 remacle Exp $
+// $Id: meshGFace.cpp,v 1.129 2008-03-26 09:37:49 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -156,7 +156,7 @@ void remeshUnrecoveredEdges(std::set<EdgeToRecover> &edgesNotRecovered,
 bool AlgoDelaunay2D(GFace *gf)
 {
   if(noseam(gf) && /*gf->getNativeType() == GEntity::GmshModel &&*/ 
-     CTX.mesh.algo2d == ALGO_2D_DELAUNAY /*&& gf->geomType() == GEntity::Plane*/)
+     (CTX.mesh.algo2d == ALGO_2D_DELAUNAY || CTX.mesh.algo2d == ALGO_2D_FRONTAL) /*&& gf->geomType() == GEntity::Plane*/)
     return true;
   return false;
 }
@@ -721,7 +721,7 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
   // vertices
   if(AlgoDelaunay2D(gf)){
     gmshBowyerWatson(gf);
-    //    laplaceSmoothing(gf);
+    laplaceSmoothing(gf);
   }
   else if (debug){
     char name[256];
@@ -1325,6 +1325,13 @@ void meshGFace::operator() (GFace *gf)
     else
       algo = "MeshAdapt+Delaunay";
     break;
+  case ALGO_2D_FRONTAL:
+    // FIXME: Delaunay not available in all cases at the moment
+    if(AlgoDelaunay2D(gf))
+      algo = "Frontal";
+    else
+      algo = "MeshAdapt+Delaunay";
+    break;
   case ALGO_2D_MESHADAPT_DELAUNAY:
     algo = "MeshAdapt+Delaunay";
     break;

Mesh/meshGFaceDelaunayInsertion.cpp

-// $Id: meshGFaceDelaunayInsertion.cpp,v 1.16 2008-03-25 20:25:35 remacle Exp $
+// $Id: meshGFaceDelaunayInsertion.cpp,v 1.17 2008-03-26 09:37:49 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -31,6 +31,9 @@
 #include <map>
 #include <algorithm>
 
+#include "Context.h"
+extern Context_T CTX;
+
 void circumCenterXY(double *p1, double *p2, double *p3, double *res)
 {
   double d, a1, a2, a3;
@@ -545,6 +548,10 @@ static void insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
 
 void gmshBowyerWatson(GFace *gf)
 {
+  if (CTX.mesh.algo2d == ALGO_2D_FRONTAL){
+    gmshBowyerWatsonFrontal(gf);
+    return;
+  }
   std::set<MTri3*,compareTri3Ptr> AllTris;
   std::vector<double> vSizes, vSizesBGM, Us, Vs;
 
@@ -635,7 +642,7 @@ static double length_metric ( const double p[2], const double q[2], const double
 */
 
 
-void gmshBowyerWatson2(GFace *gf){
+void gmshBowyerWatsonFrontal(GFace *gf){
 //void gmshFrontalDelaunay (GFace *gf){
   std::set<MTri3*,compareTri3Ptr> AllTris;
   std::vector<double> vSizes, vSizesBGM, Us, Vs;
@@ -692,21 +699,29 @@ void gmshBowyerWatson2(GFace *gf){
     // we try to find a point that would produce a perfect triangle while
     // connecting the 2 points of the active edge
 
-    const double p    = 0.5*length_metric (P,Q,metric);
-    const double q    = length_metric (center,midpoint,metric);
-    const double rhoM = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] + vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3);
+
+    double dir[2] = {center[0]-midpoint[0],center[1]-midpoint[1]};
+    double q = sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
+    dir[0]/=q;
+    dir[1]/=q;
+    const double RATIO = sqrt(  dir[0]*dir[0]*metric[0]+
+			      2*dir[1]*dir[0]*metric[1]+
+			        dir[1]*dir[1]*metric[2]);    
+
+    const double p    = 0.5*sqrt(DSQR(P[0]-Q[0])+DSQR(P[1]-Q[1]));//length_metric (P,Q,metric);
+    //    const double q    = length_metric (center,midpoint,metric);
+    const double rhoM = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] + vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3) / RATIO;
 
     //    printf("%g vs %g\n",2*p,rhoM);
 
     //    const double rhoM_hat = std::max(rhoM,2*p);
     const double rhoM_hat = std::min(std::max(rhoM,p),(p*p+q*q)/(2*q));
     const double d = rhoM_hat + sqrt (rhoM_hat*rhoM_hat - p*p);
-    double dir[2] = {center[0]-midpoint[0],center[1]-midpoint[1]};
-    double L = sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
+    
     double newPoint[2] = 
       {
-	midpoint[0] + d * dir[0]/L,
-	midpoint[1] + d * dir[1]/L
+	midpoint[0] + d * dir[0],
+	midpoint[1] + d * dir[1]
       };
 
     //    printf("%g %g -- %g %g -- %g %g\n",midpoint[0],midpoint[1],pa[0],pa[1],newPoint[0],newPoint[1]);

Mesh/meshGFaceDelaunayInsertion.h

@@ -106,6 +106,7 @@ void connectTriangles(std::list<MTri3*> &);
 void connectTriangles(std::vector<MTri3*> &);
 void connectTriangles(std::set<MTri3*,compareTri3Ptr> &AllTris);
 void gmshBowyerWatson(GFace *gf);
+void gmshBowyerWatsonFrontal(GFace *gf);
 
 struct edgeXface
 {

Plugin/GeoEarthImport.cpp

@@ -31,7 +31,8 @@ void GeoEarthImport::add_point(const SPoint3 &point){
   SPoint2 stereo(-point.x()/(1+point.z()),-point.y()/(1+point.z()));
   midpoint+=lastpoint;
   midpoint/=2;  
-  if( ip==first_point_in_loop || !size_field || point.distance(lastpoint)*6361e3 > (*size_field)(midpoint[0],midpoint[1],midpoint[2])){
+  //  printf("%g\n",sqrt(midpoint[0]*midpoint[0]+midpoint[1]*midpoint[1]+midpoint[2]*midpoint[2]));
+  if( ip==first_point_in_loop || !size_field || point.distance(lastpoint)*6361e3 > (*size_field)(midpoint[0]*6361e3,midpoint[1]*6361e3,midpoint[2]*6361e3)){
     loop_buff<<"Point ( IP + "<<ip++<<" ) = {"<<stereo.x()<<", "<<stereo.y()<<", "<<0<<", 1};\n";
     lastpoint=point;
   }

demos/gen_earth_100km.geo

 // you need gshhs data file, available on ftp://ftp.soest.hawaii.edu/pwessel/gshhs/gshhs_1.6.zip
-Field[1] = MathEval;
-Field[1].F= "100.e3";
-/*
+Field[1] = Box;
+Field[1].VIn = 2000;
+Field[1].VOut = 100000;
+Field[1].XMax = 6000000;
+Field[1].XMin = 3000000;
+Field[1].YMin = -2000000;
+Field[1].YMax = 1000000;
+Field[1].ZMin = 3220000;
+Field[1].ZMax = 6000000;
+//Field[2] = MathEval;
+//Field[2].F= "1.e3";
+
+
 Plugin(GSHHS).iField=1;
 Plugin(GSHHS).InFileName="gshhs_f.b";
-Plugin(GSHHS).OutFileName="earth_100km.geo";
+Plugin(GSHHS).OutFileName="earth_adapt.geo";
 Plugin(GSHHS).Run;
-*/
-Merge "earth_100km.geo";
+/*
+Merge "earth_4km.geo";
 Field[3] = Threshold;
-Field[3].LcMin = 10e3;
+Field[3].LcMin = 4e3;
 Field[3].LcMax = 500e3;
 Field[3].DistMax = 2000e3;
 Field[3].DistMin = 0e3;
@@ -17,4 +27,4 @@ Field[3].IField = 2;
 Field[4]=MathEval;
 Field[4].F = "80e3";
 Background Field = 3;
-
+*/