From bcd0305a69e2ec20237c323627cef5e383e3cff5 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Sun, 16 Dec 2001 05:16:38 +0000
Subject: [PATCH] *** empty log message ***

---
 Common/Options.cpp             |   4 +-
 Mesh/3D_Mesh.cpp               | 136 +++++++++++++++++++++++----------
 Mesh/Makefile                  |   4 +-
 Mesh/Matrix.h                  | 125 ++++++++++++------------------
 Mesh/Metric.cpp                |  99 ++++++++++++++++++++----
 Mesh/Utils.cpp                 |   6 +-
 Parallel/Makefile              |   5 +-
 Parser/Gmsh.tab.cpp            |   8 +-
 Parser/Gmsh.yy.cpp             |  13 +++-
 benchmarks/2d/Square-Attr1.geo |   5 +-
 benchmarks/2d/conge.geo        |   2 +-
 benchmarks/2d/naca12_2d.geo    |   5 +-
 benchmarks/3d/Cube-02.geo      |   9 ++-
 examples/cube.geo              |  83 ++++++++++++++++++++
 14 files changed, 345 insertions(+), 159 deletions(-)
 create mode 100644 examples/cube.geo

diff --git a/Common/Options.cpp b/Common/Options.cpp
index b680e62628..5b6f4bad60 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.66 2001-12-03 08:41:43 geuzaine Exp $
+// $Id: Options.cpp,v 1.67 2001-12-16 05:16:37 remacle Exp $
 
 #include "Gmsh.h"
 #include "GmshUI.h"
@@ -1546,7 +1546,7 @@ double opt_mesh_constrained_bgmesh(OPT_ARGS_NUM){
 }
 double opt_mesh_degree(OPT_ARGS_NUM){
   if(action & GMSH_SET)
-    CTX.mesh.degree = 1; //(int)val;  INTERDIT POUR LE MOMENT !!!
+    CTX.mesh.degree = (int)val; // INTERDIT POUR LE MOMENT !!!
 #ifdef _FLTK
   if(WID && (action & GMSH_GUI))
     WID->mesh_butt[3]->value(CTX.mesh.degree==2);
diff --git a/Mesh/3D_Mesh.cpp b/Mesh/3D_Mesh.cpp
index a41360c6dd..00441b6653 100644
--- a/Mesh/3D_Mesh.cpp
+++ b/Mesh/3D_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: 3D_Mesh.cpp,v 1.33 2001-11-16 19:35:26 remacle Exp $
+// $Id: 3D_Mesh.cpp,v 1.34 2001-12-16 05:16:37 remacle Exp $
 
 /*
  
@@ -18,13 +18,13 @@
 
 #include <time.h>
 #include <vector>
+#include <algorithm>
 #include "Gmsh.h"
 #include "Numeric.h"
 #include "Mesh.h"
 #include "3D_Mesh.h"
 #include "Create.h"
 #include "Context.h"
-#include "mGeomSearch.h"
 
 extern Mesh       *THEM, *LOCAL;
 extern Context_T   CTX;
@@ -42,6 +42,39 @@ static int ZONEELIMINEE, Methode = 0;
 Simplex  MyNewBoundary;
 int      Alerte_Point_Scabreux;
 
+inline void cgsmpl (Simplex *s, double &x, double &y, double &z)
+{
+  x = 0.25 * ( s->V[0]->Pos.X +
+	       s->V[1]->Pos.X +
+	       s->V[2]->Pos.X +
+	       s->V[3]->Pos.X);  
+  y = 0.25 * ( s->V[0]->Pos.Y +
+	       s->V[1]->Pos.Y +
+	       s->V[2]->Pos.Y +
+	       s->V[3]->Pos.Y);  
+  z = 0.25 * ( s->V[0]->Pos.Z +
+	       s->V[1]->Pos.Z +
+	       s->V[2]->Pos.Z +
+	       s->V[3]->Pos.Z);    
+}
+
+struct closestSimplex 
+{
+  double X,Y,Z;
+  closestSimplex (double x, double y, double z)
+    : X(x),Y(y),Z(z){}
+  bool operator () (Simplex *a, Simplex *b)
+  {
+    double x1,y1,z1,x2,y2,z2;
+    cgsmpl (a,x1,y1,z1);
+    cgsmpl (b,x2,y2,z2);
+    double d1 = (x1-X)*(x1-X) + (y1-Y)*(y1-Y) + (z1-Z)*(z1-Z);
+    double d2 = (x2-X)*(x2-X) + (y2-Y)*(y2-Y) + (z2-Z)*(z2-Z);
+    return d1<d2;
+  }
+};
+
+
 
 void DebugSimplexe (Simplex * s){
   int i;
@@ -199,7 +232,33 @@ struct SimplexInBox
   }
 };
 
-AOMD::mGeomSearch < Simplex, SimplexInBox, SimplexInteriorCheck > *fastSearch = 0;
+/*
+  Recursive search;
+ */
+
+Simplex *SearchPointByNeighbor (Vertex *v, Simplex *s, Tree_T *visited, int depth)
+{
+  if(depth > 150)return 0;
+  if(Tree_Search(visited,&s))return 0;
+  Tree_Add(visited,&s);
+  if(Pt_In_Circum (s,v)) return s;
+  Simplex *S[4];
+  int k=0;
+  for( int i = 0 ; i < 4 ; i++ )
+    {
+      if (s->S[i] != &MyNewBoundary)
+	{
+	  S[k++] = s->S[i];
+	}
+    }
+  std::sort(S,S+k,closestSimplex(v->Pos.X,v->Pos.Y,v->Pos.Z));
+  for( int i = 0 ; i < k ; i++ )
+    {
+      Simplex *q = SearchPointByNeighbor (v,S[i],visited, depth+1);
+      if (q) return q;
+    }
+  return 0;
+}
 
 void Action_First_Simplexes (void *a, void *b){
   Simplex **q;
@@ -574,7 +633,6 @@ bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){
       Simplex *theNewS;
       List_Read (Simplexes_New, i, &theNewS);
       Tree_Add (m->Simplexes, &theNewS);
-      if(fastSearch)fastSearch->AddObject(theNewS);
     }
     
     /* Suppression des simplexes elimines */
@@ -583,7 +641,6 @@ bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){
       List_Read (Simplexes_Destroyed, i, &s);
       if (!Tree_Suppress (m->Simplexes, &s))
         Msg(GERROR, "Impossible to delete simplex");
-      if(fastSearch)fastSearch->RemoveObject(s);
       Free_Simplex (&s,0);
     }
     
@@ -600,8 +657,6 @@ bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){
 
 void Convex_Hull_Mesh (List_T * Points, Mesh * m){
   int i, j, N, n;
-  int Nbr_OK = 0, Nbr_NOTOK = 0;
-  double xxx[3],uuu[3];
 
   N = List_Nbr (Points);
   n = IMAX (N / 20, 1);
@@ -609,54 +664,51 @@ void Convex_Hull_Mesh (List_T * Points, Mesh * m){
   clock_t t1 = clock();
 
   Box_6_Tetraedron (Points, m);
-  //  List_Sort (Points, comparePosition);
-
-  SimplexInteriorCheck ick;
-  SimplexInBox         ibx (LC3D/30.);
+  List_Sort (Points, comparePosition);
 
-    fastSearch = new AOMD::mGeomSearch < Simplex, SimplexInBox, SimplexInteriorCheck >
-    (  m->Grid.min.X,m->Grid.max.X,
-       m->Grid.min.Y,m->Grid.max.Y,
-       m->Grid.min.Z,m->Grid.max.Z,
-       ibx,ick, 30,30,30);
-  
-  int nbEasy = 0, nbFast = 0;
+  int Nb1 = 0, Nb2 = 0, Nb3 = 0;
 
   for (i = 0; i < N; i++){
     THES = NULL;
     List_Read (Points, i, &THEV);
-    xxx[0]= THEV->Pos.X; 
-    xxx[1]= THEV->Pos.Y; 
-    xxx[2]= THEV->Pos.Z; 
-    if(fastSearch) THES = fastSearch->find(xxx,uuu); 
+    if (Simplexes_New)
+      for (j = 0; j < List_Nbr (Simplexes_New); j++)
+	{
+	  Simplex *simp;
+	  List_Read (Simplexes_New, j, &simp);
+	  if(Pt_In_Circum(simp,THEV))
+	    {
+	      THES = simp;
+		  break;
+	    }
+	}    
+    
+    if(THES) Nb1 ++;
+
     if(!THES)
-      {   
-	if (Simplexes_New)
-	  for (j = 0; j < List_Nbr (Simplexes_New); j++){
-	    Action_First_Simplexes (List_Pointer (Simplexes_New, j), NULL);	
-	    if(THES)
-	      {
-		nbEasy++;
-		break;
-	      }
-	  }
-      }    
-    else
       {
-	nbFast++;
+	if (Simplexes_New &&  List_Nbr (Simplexes_New))
+	  {
+	    Tree_T *visited = Tree_Create (sizeof (Simplex *), compareSimplex);
+	    Simplex *simp;
+	    List_Read (Simplexes_New, 0, &simp);
+	    THES = SearchPointByNeighbor (THEV, simp, visited,0);
+	    Tree_Delete(visited);
+	  }
+	if(THES) Nb2 ++;
       }
+
+
     if (!THES){
       Tree_Action (m->Simplexes, Action_First_Simplexes);
-      Nbr_OK++;
-    }
-    else{
-      Nbr_NOTOK++;
+      if(THES) Nb3 ++;
     }
+
     if (i % n == n - 1){
       volume = 0.0;
       Tree_Action (m->Simplexes, VSIM);
       Msg(STATUS3, "Nod=%d/%d Elm=%d", i+1,N,Tree_Nbr(m->Simplexes)); 
-      Msg(STATUS1, "Vol=%g (%d %d %d)",volume,nbFast,nbEasy,i-nbFast-nbEasy); 
+      Msg(STATUS1, "Vol=%g (%d %d %d)",volume,Nb1,Nb2,Nb3); 
     }
     if (!THES){
       Msg(WARNING, "Vertex (%g,%g,%g) in no simplex", THEV->Pos.X,THEV->Pos.Y,THEV->Pos.Z); 
@@ -691,8 +743,8 @@ void Convex_Hull_Mesh (List_T * Points, Mesh * m){
   }
   clock_t t2 = clock();
 
-  if(fastSearch)delete fastSearch;
-  printf("nbEasy = %d nbFast = %d N = %d t = %lf\n",nbEasy,nbFast,N,(double)(t2-t1)/CLOCKS_PER_SEC);
+  Msg(STATUS3,"Nb1 = %d Nb2 = %d Nb3 = %d N = %d t = %lf\n",Nb1,Nb2,Nb3
+      ,N,(double)(t2-t1)/CLOCKS_PER_SEC);
 }
 
 void suppress_vertex (void *data, void *dum){
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 463cb5c71d..43942c591b 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.38 2001-12-04 10:47:39 geuzaine Exp $
+# $Id: Makefile,v 1.39 2001-12-16 05:16:37 remacle Exp $
 #
 # Makefile for "libMesh.a"
 #
@@ -38,7 +38,7 @@ SRC = 1D_Mesh.cpp \
         2D_Parametric.cpp \
         2D_Mesh_Aniso.cpp \
         2D_Mesh_Shewchuk.cpp \
-      3D_Mesh_Old.cpp \
+        3D_Mesh.cpp \
         3D_SMesh.cpp \
         3D_BGMesh.cpp \
         3D_Extrude.cpp \
diff --git a/Mesh/Matrix.h b/Mesh/Matrix.h
index 36ca7e3539..aad2f67475 100644
--- a/Mesh/Matrix.h
+++ b/Mesh/Matrix.h
@@ -90,10 +90,6 @@ class Matrix3x3
 	  cout << "| " << mat[2][0] << " " << mat[2][1] << " " << mat[2][2] << "|"<<endl;
 	  */
 	}
-      Matrix3x3 ()
-      {
-	zero = 0.0;
-      }
       Matrix3x3 (Vector3 &v1, Vector3 &v2, Vector3 &v3)
       {
 	setColumn(0,v1);
@@ -112,30 +108,17 @@ class Matrix3x3
 	  (eVec[i])[2] = v3;
 	  eVal[i] = l;
 	}
-      void setMetric ()
+      void setMetric ();
+
+      Matrix3x3 (double init = 0.0)
 	{
-	  Matrix3x3 rot (eVec[0],eVec[1],eVec[2]);
-	  Matrix3x3 rotT (eVec[0],eVec[1],eVec[2]);
-	  Matrix3x3 Id (0.0);
-	  rotT.transpose();
-	  Id.Identity(1.0);
-	  for(int i=0;i<3;i++)Id.mat[i][i] = eVal[i];
-	  Matrix3x3 m = rotT * (Id * rot);
+	  zero = init;
 	  for(int i=0;i<_TAILLE_;i++)
 	    for(int j=0;j<_TAILLE_;j++)
-	      mat[i][j] = m.mat[i][j];
+	      mat[i][j] = zero;
 	}
-
-      Matrix3x3 (const double& init)
-      {
-            zero = init;
-      		for(int i=0;i<_TAILLE_;i++)
-      			for(int j=0;j<_TAILLE_;j++)
-                	mat[i][j] = zero;
-      }
       Matrix3x3 (const double& init, double z[3][3])
-      {
-	
+      {	
 	Identity(init);
 	zero = init;
 	for(int i=0;i<_TAILLE_;i++)
@@ -144,41 +127,44 @@ class Matrix3x3
       }
       Matrix3x3 operator + (const Matrix3x3 &autre)
       {
-      		Matrix3x3 m(0.);
-      		for(int i=0;i<_TAILLE_;i++)
-      			for(int j=0;j<_TAILLE_;j++)
-                	m.mat[i][j] = mat[i][j] + autre.mat[i][j];
-      		return m;
+	Matrix3x3 m(0.);
+	for(int i=0;i<_TAILLE_;i++)
+	  for(int j=0;j<_TAILLE_;j++)
+	    m.mat[i][j] = mat[i][j] + autre.mat[i][j];
+	return m;
       }
       Matrix3x3 operator - (const Matrix3x3 &autre)
-      {
-            Matrix3x3 m(0.);
-      		for(int i=0;i<_TAILLE_;i++)
-      			for(int j=0;j<_TAILLE_;j++)
-                	m.mat[i][j] = mat[i][j] - autre.mat[i][j];
-      		return m;
-      }
+	{
+	  Matrix3x3 m(0.);
+	  for(int i=0;i<_TAILLE_;i++)
+	    for(int j=0;j<_TAILLE_;j++)
+	      m.mat[i][j] = mat[i][j] - autre.mat[i][j];
+	  return m;
+	}
       double* operator [] (int i)
-      {
-             return mat[i];
-      }
+	{
+	  return mat[i];
+	}
       void Identity(double id)
-      {
-	for(int i=0;i<_TAILLE_;i++)
-	  mat[i][i] = id;
-	setEigen(0,1,0,0,1);
-	setEigen(1,0,1,0,1);
-	setEigen(2,0,0,1,1);
-      }
-
+	{
+	  for(int i=0;i<_TAILLE_;i++)
+	    for(int j=0;j<_TAILLE_;j++)
+	      {
+		if(i == j)mat[i][j] = id;
+		else mat[i][j] = 0.0;
+	      }
+	  setEigen(0,1,0,0,1);
+	  setEigen(1,0,1,0,1);
+	  setEigen(2,0,0,1,1);
+	}
+      
       void copy(double m[3][3])
       {
-
-      		for(int i=0;i<_TAILLE_;i++)
-      			for(int j=0;j<_TAILLE_;j++)
-                    m[i][j] = mat[i][j];
+	for(int i=0;i<_TAILLE_;i++)
+	  for(int j=0;j<_TAILLE_;j++)
+	    m[i][j] = mat[i][j];
       }
-
+      
       Vector3 getColumn (int i)
 	{
 	  return Vector3 (0.0,mat[0][i],mat[1][i],mat[2][i]);
@@ -200,19 +186,7 @@ class Matrix3x3
 	  mat[i][j] = k;
 	}
 
-      Matrix3x3 operator * (const Matrix3x3 &autre)
-      {
-
-                Matrix3x3 m(0.);
-      		for(int i=0;i<_TAILLE_;i++)
-                for(int j=0;j<_TAILLE_;j++)
-                {
-                    m.mat[i][j] = zero;
-                    for(int k=0;k<_TAILLE_;k++)
-                    	m.mat[i][j] += mat[i][k] * autre.mat[k][j];
-                }
-      		return m;
-      }
+      Matrix3x3 operator * (const Matrix3x3 &autre);
       /*
       bool invert ()
       {
@@ -237,19 +211,18 @@ class Matrix3x3
 
       void transpose()
       {
-            double temp;
-            for(int i=0;i<_TAILLE_;i++)
-                for(int j=0;j<i;j++)
-                {
-                    if(i!=j)
-                    {
-                    	temp = mat[i][j];
-                        mat[i][j] = mat[j][i];
-                        mat[j][i] = temp;
-                    }
-                }
+	double temp;
+	for(int i=0;i<_TAILLE_;i++)
+	  for(int j=0;j<i;j++)
+	    {
+	      if(i!=j)
+		{
+		  temp = mat[i][j];
+		  mat[i][j] = mat[j][i];
+		  mat[j][i] = temp;
+		}
+	    }
       }
 };
 
-#undef _TAILLE_
 #endif
diff --git a/Mesh/Metric.cpp b/Mesh/Metric.cpp
index 54f16c279d..8cd744aa9a 100644
--- a/Mesh/Metric.cpp
+++ b/Mesh/Metric.cpp
@@ -1,5 +1,5 @@
-// $Id: Metric.cpp,v 1.7 2001-09-04 16:25:05 geuzaine Exp $
-
+// $Id: Metric.cpp,v 1.8 2001-12-16 05:16:37 remacle Exp $
+#include <time.h>
 #include "Gmsh.h"
 #include "Numeric.h"
 #include "Geo.h"
@@ -8,6 +8,61 @@
 #include "Matrix.h"
 #include "Interpolation.h"
 
+Matrix3x3 Matrix3x3::operator * (const Matrix3x3 &autre)
+{  
+  Matrix3x3 m(0.);
+  for(int i=0;i<3;i++)
+    {
+      for(int j=0;j<3;j++)
+	{
+	  for(int k=0;k<3;k++)
+	    m.mat[i][j] += mat[i][k] * autre.mat[k][j];
+	  if (fabs(m.mat[i][j]) < 1.e-15)m.mat[i][j] = 0.0;
+	}
+    }
+  return m;
+}
+
+void Matrix3x3::setMetric()
+{
+  Matrix3x3 rot (eVec[0],eVec[1],eVec[2]);
+  Matrix3x3 rotT (eVec[0],eVec[1],eVec[2]);
+  Matrix3x3 Id (0.0);
+  rotT.transpose();
+  Id.Identity(1.0);
+  for(int i=0;i<3;i++)Id.mat[i][i] = eVal[i];
+  Matrix3x3 kk = (Id * rot);
+  Matrix3x3 m = rotT * kk;
+  /*
+  Msg(INFO,"rot: %g %g %g %g %g %g %g %g %g \n  ", 
+      rot[0][0], rot[0][1], rot[0][2],
+      rot[1][0], rot[1][1], rot[1][2],
+      rot[2][0], rot[2][1], rot[2][2]);	      	      
+  Msg(INFO,"rotT: %g %g %g %g %g %g %g %g %g \n  ", 
+      rotT[0][0], rotT[0][1], rotT[0][2],
+      rotT[1][0], rotT[1][1], rotT[1][2],
+      rotT[2][0], rotT[2][1], rotT[2][2]);	      	      
+  Msg(INFO,"id: %g %g %g %g %g %g %g %g %g \n  ", 
+      Id[0][0], Id[0][1], Id[0][2],
+      Id[1][0], Id[1][1], Id[1][2],
+      Id[2][0], Id[2][1], Id[2][2]);	      	      
+
+  Msg(INFO,"kk: %g %g %g %g %g %g %g %g %g \n  ", 
+      kk[0][0], kk[0][1], kk[0][2],
+      kk[1][0], kk[1][1], kk[1][2],
+      kk[2][0], kk[2][1], kk[2][2]);	      	      
+
+  Msg(INFO,"m: %g %g %g %g %g %g %g %g %g \n  ", 
+      m[0][0], m[0][1], m[0][2],
+      m[1][0], m[1][1], m[1][2],
+      m[2][0], m[2][1], m[2][2]);	      	      
+  */
+  for(int i=0;i<3;i++)
+    for(int j=0;j<3;j++)
+      mat[i][j] = m.mat[i][j];
+}
+
+
 GMSHMetric::GMSHMetric ()
 {
   Identity ();
@@ -99,6 +154,7 @@ Local_Metric_Of_Attractors ( double X, double Y, double Z,
 	  
 	}
       double d1 = d * a->Radius;
+      if(fabs(d1) < a->Radius * 1.e-6)d1 = 0.0;
       u = exp (-(d1 * d1));
       x1 = (1. - u) + u * a->lc1;
       x2 = (1. - u) + u * a->lc2;
@@ -137,6 +193,19 @@ Local_Metric_Of_Attractors ( double X, double Y, double Z,
 	      myMetric.setEigen (1,d2.Pos.X,d2.Pos.Y,d2.Pos.Z,1./(x2*x2));
 	      myMetric.setEigen (2,d3.Pos.X,d3.Pos.Y,d3.Pos.Z,1./(x1*x1));
 	      myMetric.setMetric();
+	      /*
+	      Msg(INFO,"%d %12.5E heigens : %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E \n  ", a->c->Num,d,
+		  metr.Pos.X, metr.Pos.Y, metr.Pos.Z,
+		  d2.Pos.X,d2.Pos.Y,d2.Pos.Z,
+		  d3.Pos.X,d3.Pos.Y,d3.Pos.Z);
+
+	      Msg(INFO,"heigenv : %12.5E %12.5E \n",x1,x2);
+
+	      Msg(INFO,"curve Metric: %g %g %g %g %g %g %g %g %g \n  ", 
+		  myMetric[0][0], myMetric[0][1], myMetric[0][2],
+		  myMetric[1][0], myMetric[1][1], myMetric[1][2],
+		  myMetric[2][0], myMetric[2][1], myMetric[2][2]);	      	      
+	      */
 	    }
 	  Matrix3x3 *M[2];
 	  M[0] = &myMetric;
@@ -145,7 +214,7 @@ Local_Metric_Of_Attractors ( double X, double Y, double Z,
 	  // myOldMetric = myMetric;
 	}
     }
-  myOldMetric.copy(m);
+  myOldMetric.copy(m);  
   return 1.0;
 }
 
@@ -285,7 +354,7 @@ setSimplexQuality (Simplex * s, Surface * surf)
       s->Quality = DMAX (DMAX (q1, q2), q3) / (RacineDeTrois);
     }
   else
-    {
+    {      
       s->Center_Ellipsum_2D (m);
       s->Quality = 3. * s->Radius / (s->V[0]->lc + s->V[1]->lc + s->V[2]->lc);
     }
@@ -381,12 +450,13 @@ getLc (double u, Curve * c)
   Vertex du = InterpolateCurve (c, u, 1);
   Local_Metric_Of_Attractors (v.Pos.X, v.Pos.Y, v.Pos.Z, NULL);
   l = LengthVector (&du);
-  /*
-  printf("GetLC : u = %g l=%g  lc=%g  return=%g\n", u, l, v.lc, l/v.lc);
-  printf("Metric: %g %g %g %g %g %g %g %g %g \n  ", 
-	 m[0][0], m[0][1], m[0][2],
-	 m[1][0], m[1][1], m[1][2],
-	 m[2][0], m[2][1], m[2][2]);
+  /*   
+    Msg(INFO,"GetLC : u = %g l=%g  lc=%g  return=%g\n", u, l, v.lc, l/v.lc);
+    Msg(INFO,"Metric: %g %g %g %g %g %g %g %g %g \n  ", 
+    m[0][0], m[0][1], m[0][2],
+    m[1][0], m[1][1], m[1][2],
+    m[2][0], m[2][1], m[2][2]);
+    Msg(INFO,"du = %g %g %g\n",du.Pos.X,du.Pos.Y,du.Pos.Z);
   */
   return l / v.lc;
 }
@@ -394,10 +464,11 @@ getLc (double u, Curve * c)
 double GMSHMetric::
 LengthVector (Vertex * v)
 {
-  Vertex mult (v->Pos.X * m[0][0] + v->Pos.Y * m[0][1] + v->Pos.Z * m[0][2],
-	       v->Pos.X * m[1][0] + v->Pos.Y * m[1][1] + v->Pos.Z * m[1][2],
-	       v->Pos.X * m[2][0] + v->Pos.Y * m[2][1] + v->Pos.Z * m[2][2]);
-  return sqrt (mult * (*v));
+  double qqq = v->Pos.X *  (v->Pos.X * m[0][0] + v->Pos.Y * m[0][1] + v->Pos.Z * m[0][2]) +
+    v->Pos.Y * (v->Pos.X * m[1][0] + v->Pos.Y * m[1][1] + v->Pos.Z * m[1][2]) +
+    v->Pos.Z * (v->Pos.X * m[2][0] + v->Pos.Y * m[2][1] + v->Pos.Z * m[2][2]);
+  //  Msg(INFO,"qqq = %g\n",qqq);
+  return sqrt (qqq);
 }
 
 double GMSHMetric::
diff --git a/Mesh/Utils.cpp b/Mesh/Utils.cpp
index b03dba3716..9799e73d40 100644
--- a/Mesh/Utils.cpp
+++ b/Mesh/Utils.cpp
@@ -1,4 +1,4 @@
-// $Id: Utils.cpp,v 1.7 2001-11-30 14:11:47 geuzaine Exp $
+// $Id: Utils.cpp,v 1.8 2001-12-16 05:16:37 remacle Exp $
 
 #include "Gmsh.h"
 #include "Numeric.h"
@@ -515,8 +515,8 @@ void RecursiveIntegration (IntPoint * from, IntPoint * to, double (*f) (double X
   val3 = trapeze (&P, to);
 
   err = fabs (val1 - val2 - val3);
-
-  if ((err < Prec) && (*depth > 1)){
+  //  Msg(INFO,"Int %22.15 E %22.15 E %22.15 E\n", val1,val2,val3);
+  if (((err < Prec) && (*depth > 1)) || (*depth > 25)){
     List_Read (pPoints, List_Nbr (pPoints) - 1, &p1);
     P.p = p1.p + val2;
     List_Add (pPoints, &P);
diff --git a/Parallel/Makefile b/Parallel/Makefile
index 52825380c3..fb031842df 100644
--- a/Parallel/Makefile
+++ b/Parallel/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.6 2001-10-30 09:52:38 geuzaine Exp $
+# $Id: Makefile,v 1.7 2001-12-16 05:16:37 remacle Exp $
 #
 # Makefile for "libParallel.a"
 #
@@ -8,14 +8,13 @@ AR       = ar ruvs
 RM       = rm
 RANLIB   = ranlib
 LIB      = ../lib/libParallel.a
+CC       = gcc
 
 OPT_FLAGS     = -g
 OS_FLAGS      = 
 VERSION_FLAGS = 
 
 RMFLAGS  = -f
-CFLAGS   = $(OPT_FLAGS) $(OS_FLAGS) $(VERSION_FLAGS) $(INCLUDE)\
-          -I/usr/local/mpich/latest/include/
 
 SRC = ParUtil.cpp
 
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index ae8fb39c40..91d3e81b6a 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -177,7 +177,7 @@
 #line 1 "Gmsh.y"
  
 
-// $Id: Gmsh.tab.cpp,v 1.133 2001-12-03 10:39:39 gyselinc Exp $
+// $Id: Gmsh.tab.cpp,v 1.134 2001-12-16 05:16:37 remacle Exp $
 
 #include <stdarg.h>
 #ifndef _NOPLUGIN
@@ -2685,7 +2685,7 @@ static const short yycheck[] = {    23,
    175,    -1,    -1,    -1,    -1,    -1,   181,    -1,   183
 };
 /* -*-C-*-  Note some compilers choke on comments on `#line' lines.  */
-#line 3 "/usr/lib/bison.simple"
+#line 3 "/usr/share/bison.simple"
 /* This file comes from bison-1.28.  */
 
 /* Skeleton output parser for bison,
@@ -2899,7 +2899,7 @@ __yy_memcpy (char *to, char *from, unsigned int count)
 #endif
 #endif
 
-#line 217 "/usr/lib/bison.simple"
+#line 217 "/usr/share/bison.simple"
 
 /* The user can define YYPARSE_PARAM as the name of an argument to be passed
    into yyparse.  The argument should have type void *.
@@ -6214,7 +6214,7 @@ case 385:
     break;}
 }
    /* the action file gets copied in in place of this dollarsign */
-#line 543 "/usr/lib/bison.simple"
+#line 543 "/usr/share/bison.simple"
 
   yyvsp -= yylen;
   yyssp -= yylen;
diff --git a/Parser/Gmsh.yy.cpp b/Parser/Gmsh.yy.cpp
index 7d4b2f1919..f21a86a637 100644
--- a/Parser/Gmsh.yy.cpp
+++ b/Parser/Gmsh.yy.cpp
@@ -2,7 +2,7 @@
 /* A lexical scanner generated by flex */
 
 /* Scanner skeleton version:
- * $Header: /cvsroot/gmsh/Parser/Gmsh.yy.cpp,v 1.133 2001-12-03 10:39:39 gyselinc Exp $
+ * $Header: /cvsroot/gmsh/Parser/Gmsh.yy.cpp,v 1.134 2001-12-16 05:16:38 remacle Exp $
  */
 
 #define FLEX_SCANNER
@@ -10,7 +10,6 @@
 #define YY_FLEX_MINOR_VERSION 5
 
 #include <stdio.h>
-#include <unistd.h>
 
 
 /* cfront 1.2 defines "c_plusplus" instead of "__cplusplus" */
@@ -24,6 +23,7 @@
 #ifdef __cplusplus
 
 #include <stdlib.h>
+#include <unistd.h>
 
 /* Use prototypes in function declarations. */
 #define YY_USE_PROTOS
@@ -1000,7 +1000,7 @@ char *yytext;
 #define INITIAL 0
 #line 2 "Gmsh.l"
 
-// $Id: Gmsh.yy.cpp,v 1.133 2001-12-03 10:39:39 gyselinc Exp $
+// $Id: Gmsh.yy.cpp,v 1.134 2001-12-16 05:16:38 remacle Exp $
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -1191,7 +1191,7 @@ YY_MALLOC_DECL
 YY_DECL
 	{
 	register yy_state_type yy_current_state;
-	register char *yy_cp = NULL, *yy_bp = NULL;
+	register char *yy_cp, *yy_bp;
 	register int yy_act;
 
 #line 63 "Gmsh.l"
@@ -2790,6 +2790,11 @@ YY_BUFFER_STATE b;
 	}
 
 
+#ifndef YY_ALWAYS_INTERACTIVE
+#ifndef YY_NEVER_INTERACTIVE
+extern int isatty YY_PROTO(( int ));
+#endif
+#endif
 
 #ifdef YY_USE_PROTOS
 void yy_init_buffer( YY_BUFFER_STATE b, FILE *file )
diff --git a/benchmarks/2d/Square-Attr1.geo b/benchmarks/2d/Square-Attr1.geo
index 6088290011..877312bdac 100644
--- a/benchmarks/2d/Square-Attr1.geo
+++ b/benchmarks/2d/Square-Attr1.geo
@@ -12,8 +12,9 @@ Line(3) = {1,4};
 Line(4) = {4,3};                    
 Line Loop(5) = {1,2,3,4};                    
 Plane Surface(6) = {5};      
-Point(11) = {0.5,0.5,-1,lc};                    
+Point(11) = {0.5,0.5,0,lc};                    
 Point(22) = {0.5,0.5,1,lc};                    
 Line(5) = {11,22};                    
-Attractor Line{5} = {1,0.1,7};                             
+//Attractor Point{1,2,3,4,11} = {.01,.01,7};                             
+Attractor Line{1,2,3,4} = {1,.1,7};                             
 Mesh.Algorithm = 2;
diff --git a/benchmarks/2d/conge.geo b/benchmarks/2d/conge.geo
index b7cb6024f3..a7557d5df6 100644
--- a/benchmarks/2d/conge.geo
+++ b/benchmarks/2d/conge.geo
@@ -72,7 +72,7 @@ Line Loop(21) = {17,-15,18,19,-20,16};
 Plane Surface(22) = {21};
 Line Loop(23) = {11,-12,13,14,1,2,-3,4,5,6,7,-8,9,10};
 Plane Surface(24) = {23,21};
-
+Recombine Surface{24,21};
 
 
 
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index 2e19f9bbb9..b89500c343 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -228,6 +228,7 @@ Line(9) = {1001,1002};
 Line Loop(10) = {-6,-5,-9,-8,-7,3,4};
 Plane Surface(11) = {10};
 
-Attractor Point{101} = {0.01,0.01,2};
-Attractor Line{1,2,3,4} = {0.05,0.05,3};
+Attractor Point{101,200} = {0.04,0.04,2};
+Attractor Line{1,2,3,4} = {.1,.001,3};
+Mesh.Algorithm = 2;
 
diff --git a/benchmarks/3d/Cube-02.geo b/benchmarks/3d/Cube-02.geo
index 9d838badf4..fc5c29197d 100644
--- a/benchmarks/3d/Cube-02.geo
+++ b/benchmarks/3d/Cube-02.geo
@@ -2,10 +2,11 @@
 cube meshed  
 non uniformly  
 *****************************/  
-Point(1) = {0.0,0.0,0.0,.04};          
-Point(2) = {1,0.0,0.0,.2};          
-Point(3) = {1,1,0.0,.1};          
-Point(4) = {0,1,0.0,.2};          
+k = .3;
+Point(1) = {0.0,0.0,0.0,.04*k};          
+Point(2) = {1,0.0,0.0,.2*k};          
+Point(3) = {1,1,0.0,.1*k};          
+Point(4) = {0,1,0.0,.2*k};          
 Line(1) = {4,3};          
 Line(2) = {3,2};          
 Line(3) = {2,1};          
diff --git a/examples/cube.geo b/examples/cube.geo
new file mode 100644
index 0000000000..decdd4a696
--- /dev/null
+++ b/examples/cube.geo
@@ -0,0 +1,83 @@
+lcar1 = .0275;
+
+length = 0.25;
+height = 1.0;
+depth = 0.25;
+
+Point(newp) = {length/2,height/2,depth,lcar1}; /* Point      1 */
+Point(newp) = {length/2,height/2,0,lcar1}; /* Point      2 */
+Point(newp) = {-length/2,height/2,depth,lcar1}; /* Point      3 */
+Point(newp) = {-length/2,-height/2,depth,lcar1}; /* Point      4 */
+Point(newp) = {length/2,-height/2,depth,lcar1}; /* Point      5 */
+Point(newp) = {length/2,-height/2,0,lcar1}; /* Point      6 */
+Point(newp) = {-length/2,height/2,0,lcar1}; /* Point      7 */
+Point(newp) = {-length/2,-height/2,0,lcar1}; /* Point      8 */
+Line(1) = {3,1};
+Line(2) = {3,7};
+Line(3) = {7,2};
+Line(4) = {2,1};
+Line(5) = {1,5};
+Line(6) = {5,4};
+Line(7) = {4,8};
+Line(8) = {8,6};
+Line(9) = {6,5};
+Line(10) = {6,2};
+Line(11) = {3,4};
+Line(12) = {8,7};
+Line Loop(13) = {-6,-5,-1,11};
+Plane Surface(14) = {13};
+Line Loop(15) = {4,5,-9,10};
+Plane Surface(16) = {15};
+Line Loop(17) = {-3,-12,8,10};
+Plane Surface(18) = {17};
+Line Loop(19) = {7,12,-2,11};
+Plane Surface(20) = {19};
+Line Loop(21) = {-4,-3,-2,1};
+Plane Surface(22) = {21};
+Line Loop(23) = {8,9,6,7};
+Plane Surface(24) = {23};
+
+Surface Loop(25) = {14,24,-18,22,16,-20};
+Complex Volume(26) = {25};
+
+n = 21;
+
+/*
+sizex = 4;
+sizey = 6;
+sizez = 4;
+*/
+sizex = n*length;
+sizey = n*height;
+sizez = n*depth;
+
+sizex = 9;
+sizey = 33;
+sizez = 9;
+
+
+Transfinite Line {2,4,9,7} = sizez;
+Transfinite Line {11,5,10,12} = sizey;
+Transfinite Line {3,1,6,8} = sizex;
+
+Transfinite Surface {14} = {5,4,3,1};
+Transfinite Surface {16} = {6,2,1,5};
+Transfinite Surface {18} = {6,2,7,8};
+Transfinite Surface {20} = {3,7,8,4};
+Transfinite Surface {22} = {3,7,2,1};
+Transfinite Surface {24} = {4,8,6,5};
+
+Recombine Surface {14,16,18,20,22,24};
+
+Transfinite Volume {26} = {3,7,2,1,4,8,6,5};
+
+Physical Surface(200) = {14,16,18,20,24};
+
+Physical Volume(100) = {26};
+/*
+Mesh.use_cut_plane = 1;
+Mesh.cut_planea = 0;
+Mesh.cut_planeb = 0;
+Mesh.cut_planec = 1;
+Mesh.cut_planed = -0.125;
+*/
\ No newline at end of file
-- 
GitLab