diff --git a/Box/Main.cpp b/Box/Main.cpp
index fd007fff8e96cc67661d1f34e94c006e59a881c4..54b488599be256f362ad8ee47d6e33fba4156b13 100644
--- a/Box/Main.cpp
+++ b/Box/Main.cpp
@@ -1,4 +1,4 @@
-// $Id: Main.cpp,v 1.57 2006-02-26 01:24:44 geuzaine Exp $
+// $Id: Main.cpp,v 1.58 2006-03-08 17:04:36 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -71,7 +71,7 @@ void Info(int level, char *arg0)
 
 // Main routine for the batch (black box) version
 
-int main(int argc, char *argv[])
+int GMSHBOX(int argc, char *argv[])
 {
   ParUtil::Instance()->init(argc, argv);
 
@@ -124,11 +124,11 @@ int main(int argc, char *argv[])
       Print_Histogram(THEM->Histogram[0]);
     }
     ParUtil::Instance()->Barrier(__LINE__, __FILE__);
-    ParUtil::Instance()->Exit();
+    //    ParUtil::Instance()->Exit();
     return 1;
   }
   ParUtil::Instance()->Barrier(__LINE__, __FILE__);
-  ParUtil::Instance()->Exit();
+  //  ParUtil::Instance()->Exit();
   return 1;
 }
 
@@ -250,9 +250,21 @@ double GetValue(char *text, double defaultval)
   else
     return atof(str);
 }
-
 bool GetBinaryAnswer(const char *question, const char *yes, const char *no, 
 		     bool defaultval)
 {
-  return defaultval;
+  if(CTX.nopopup || CTX.batch )
+    return defaultval;
+
+  char answ[256];
+
+  while(1)
+    {
+      printf("%s (%s/%s)",question,yes,no);
+      
+      scanf("%s ",answ);
+      if (!strcmp(answ,yes))return true;
+      if (!strcmp(answ,no))return false;
+    }
 }
+
diff --git a/BoxMain/BoxMain.cpp b/BoxMain/BoxMain.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..af6616f2b62cb5a4db02f85e5b4aff3add676823
--- /dev/null
+++ b/BoxMain/BoxMain.cpp
@@ -0,0 +1,5 @@
+#include "Box.h"
+int main(int argc, char *argv[])
+{
+  return GMSHBOX(argc, argv);
+}
diff --git a/BoxMain/Makefile b/BoxMain/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..640b59dbc86f458fea92fcca0b8667a2e24dedcf
--- /dev/null
+++ b/BoxMain/Makefile
@@ -0,0 +1,62 @@
+# $Id: Makefile,v 1.1 2006-03-08 17:04:36 remacle Exp $
+#
+# Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+# USA.
+# 
+# Please report all bugs and problems to <gmsh@geuz.org>.
+
+include ../variables
+
+LIB     = ../lib/libGmshBoxMain.a
+INCLUDE = -I../Box
+CFLAGS  = ${OPTIM} ${FLAGS} ${INCLUDE} 
+
+SRC = BoxMain.cpp
+
+OBJ = ${SRC:.cpp=.o}
+
+.SUFFIXES: .o .cpp
+
+${LIB}: ${OBJ} 
+	${AR} ${LIB} ${OBJ} 
+	${RANLIB} ${LIB}
+
+.cpp.o:
+	${CXX} ${CFLAGS} -c $<
+
+clean:
+	rm -f *.o
+
+depend:
+	(sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \
+	${CXX} -MM ${CFLAGS} ${SRC} \
+	) >Makefile.new
+	cp Makefile Makefile.bak
+	cp Makefile.new Makefile
+	rm -f Makefile.new
+
+# DO NOT DELETE THIS LINE
+Main.o: Main.cpp ../Parallel/ParUtil.h ../Plugin/PluginManager.h \
+ ../Plugin/Plugin.h ../Common/Options.h ../Common/Message.h ../Common/Views.h \
+ ../Common/ColorTable.h ../DataStr/List.h ../Common/VertexArray.h \
+ ../Common/SmoothNormals.h ../Common/GmshMatrix.h ../Common/AdaptiveViews.h \
+ ../Common/Gmsh.h ../DataStr/Malloc.h ../DataStr/Tree.h ../DataStr/avl.h \
+ ../DataStr/Tools.h ../Common/GmshVersion.h ../Numeric/Numeric.h ../Geo/Geo.h \
+ ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h \
+ ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/Metric.h \
+ ../Mesh/Mesh.h ../Mesh/Matrix.h ../Parser/Parser.h ../Common/Context.h \
+ ../Parser/OpenFile.h ../Common/CommandLine.h
diff --git a/Geo/CAD.cpp b/Geo/CAD.cpp
index cefdc21274156c4a39773a94a97246c02d147cf8..3e8474cdff96413b800724ddda536ff0437c0935 100644
--- a/Geo/CAD.cpp
+++ b/Geo/CAD.cpp
@@ -1,4 +1,4 @@
-// $Id: CAD.cpp,v 1.95 2006-02-02 13:53:57 geuzaine Exp $
+// $Id: CAD.cpp,v 1.96 2006-03-08 17:02:50 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -2201,13 +2201,15 @@ bool ProjectPointOnCurve(Curve * c, Vertex * v, Vertex * RES, Vertex * DER)
   *RES = InterpolateCurve(CURVE, xmin, 0);
   *DER = InterpolateCurve(CURVE, xmin, 1);
   if(xmin > c->uend) {
+    xmin = c->uend;
     *RES = InterpolateCurve(CURVE, c->uend, 0);
     *DER = InterpolateCurve(CURVE, c->uend, 1);
   }
   else if(xmin < c->ubeg) {
+    xmin = c->ubeg;
     *RES = InterpolateCurve(CURVE, c->ubeg, 0);
     *DER = InterpolateCurve(CURVE, c->ubeg, 1);
-  }
+  }  
   return true;
 }
 
@@ -2308,6 +2310,8 @@ bool ProjectPointOnSurface(Surface * s, Vertex & p)
   p.Pos.X = vv.Pos.X;
   p.Pos.Y = vv.Pos.Y;
   p.Pos.Z = vv.Pos.Z;
+  p.us[0] = x[1];
+  p.us[1] = x[2];
   if(!check) {
     return false;
   }
diff --git a/Geo/Makefile b/Geo/Makefile
index 7173ab1b3d0ab8b06ac07a9bf89d91bfe5f25c9c..159a17abb022eb96fe9520a65c028f17e9fec4d9 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.72 2006-01-28 21:13:35 geuzaine Exp $
+# $Id: Makefile,v 1.73 2006-03-08 17:02:50 remacle Exp $
 #
 # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 #
@@ -19,6 +19,7 @@
 # 
 # Please report all bugs and problems to <gmsh@geuz.org>.
 
+
 include ../variables
 
 LIB     = ../lib/libGmshGeo.a
@@ -27,11 +28,21 @@ INCLUDE = -I../Common -I../DataStr -I../Geo -I../Mesh -I../Numeric\
           -I../contrib/NR
 CFLAGS  = ${OPTIM} ${FLAGS} ${INCLUDE}
 
+ifdef DEVROOT
+ INCLUDE := $(INCLUDE) -I$(DEVROOT)/MeshAdapt/model/include
+endif
+
+
 SRC = CAD.cpp \
       MinMax.cpp \
       ExtrudeParams.cpp \
       Geo.cpp \
       GeoUtils.cpp \
+      gmshModel.cpp\
+      gmshEdge.cpp\
+      gmshFace.cpp\
+      SVector3.cpp\
+      SBoundingBox3d.cpp\
       ExtractContour.cpp \
       Print_Geo.cpp
 
@@ -58,7 +69,6 @@ depend:
 	rm -f Makefile.new
 
 # DO NOT DELETE THIS LINE
-# 1 "/Users/geuzaine/.gmsh/Geo//"
 CAD.o: CAD.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../DataStr/List.h ../DataStr/Tree.h ../Numeric/Numeric.h Geo.h \
@@ -71,12 +81,10 @@ CAD.o: CAD.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../Mesh/Interpolation.h ../Mesh/Vertex.h ../Mesh/Mesh.h \
   ../Mesh/Create.h ../Mesh/Vertex.h ../Mesh/Mesh.h CAD.h ExtrudeParams.h \
   ../Common/Visibility.h ../Common/Context.h
-# 1 "/Users/geuzaine/.gmsh/Geo//"
 MinMax.o: MinMax.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
   ../Numeric/Numeric.h ../Mesh/Vertex.h ../Common/Context.h
-# 1 "/Users/geuzaine/.gmsh/Geo//"
 ExtrudeParams.o: ExtrudeParams.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
@@ -87,7 +95,6 @@ ExtrudeParams.o: ExtrudeParams.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../Common/VertexArray.h ../Common/SmoothNormals.h ../Numeric/Numeric.h \
   ../Mesh/Metric.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h \
   ../Mesh/Matrix.h ExtrudeParams.h
-# 1 "/Users/geuzaine/.gmsh/Geo//"
 Geo.o: Geo.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../DataStr/List.h ../DataStr/Tree.h ../Numeric/Numeric.h Geo.h CAD.h \
@@ -98,7 +105,6 @@ Geo.o: Geo.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Vertex.h \
   ../Mesh/Simplex.h ../Mesh/Mesh.h ../Mesh/Matrix.h ExtrudeParams.h \
   ../Parser/Parser.h ../Common/Context.h
-# 1 "/Users/geuzaine/.gmsh/Geo//"
 GeoUtils.o: GeoUtils.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
@@ -109,7 +115,31 @@ GeoUtils.o: GeoUtils.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../Common/VertexArray.h ../Common/SmoothNormals.h ../Numeric/Numeric.h \
   ../Mesh/Metric.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h \
   ../Mesh/Matrix.h ExtrudeParams.h
-# 1 "/Users/geuzaine/.gmsh/Geo//"
+gmshModel.o: gmshModel.cpp gmshModel.h gmshDefs.h ../Mesh/Mesh.h \
+  ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../Mesh/Vertex.h \
+  ../Mesh/Element.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Vertex.h \
+  ../Mesh/Element.h ../Mesh/Face.h ../Mesh/Vertex.h ../Mesh/Element.h \
+  ../Mesh/Edge.h ../Mesh/Vertex.h ../Mesh/Simplex.h \
+  ../Geo/ExtrudeParams.h ../Common/VertexArray.h \
+  ../Common/SmoothNormals.h ../Numeric/Numeric.h ../Mesh/Metric.h \
+  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h ../Mesh/Matrix.h \
+  Range.h Pair.h SPoint2.h SPoint3.h SVector3.h SBoundingBox3d.h \
+  gmshVertex.h gmshEdge.h gmshFace.h ../Parser/OpenFile.h \
+  ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
+  ../Common/Message.h
+gmshEdge.o: gmshEdge.cpp gmshModel.h gmshDefs.h ../Mesh/Mesh.h \
+  ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../Mesh/Vertex.h \
+  ../Mesh/Element.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Vertex.h \
+  ../Mesh/Element.h ../Mesh/Face.h ../Mesh/Vertex.h ../Mesh/Element.h \
+  ../Mesh/Edge.h ../Mesh/Vertex.h ../Mesh/Simplex.h \
+  ../Geo/ExtrudeParams.h ../Common/VertexArray.h \
+  ../Common/SmoothNormals.h ../Numeric/Numeric.h ../Mesh/Metric.h \
+  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h ../Mesh/Matrix.h \
+  Range.h Pair.h SPoint2.h SPoint3.h SVector3.h SBoundingBox3d.h \
+  gmshVertex.h gmshEdge.h gmshFace.h ../Mesh/Interpolation.h \
+  ../Mesh/Vertex.h ../Mesh/Mesh.h CAD.h ExtrudeParams.h Geo.h
+SVector3.o: SVector3.cpp SVector3.h SPoint3.h
+SBoundingBox3d.o: SBoundingBox3d.cpp SBoundingBox3d.h SPoint3.h
 ExtractContour.o: ExtractContour.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
@@ -120,7 +150,6 @@ ExtractContour.o: ExtractContour.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../Common/VertexArray.h ../Common/SmoothNormals.h ../Numeric/Numeric.h \
   ../Mesh/Metric.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h \
   ../Mesh/Matrix.h CAD.h ExtrudeParams.h
-# 1 "/Users/geuzaine/.gmsh/Geo//"
 Print_Geo.o: Print_Geo.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
diff --git a/Geo/Pair.h b/Geo/Pair.h
new file mode 100644
index 0000000000000000000000000000000000000000..d5edd5b83b4e693ccb2bbf018db3691ddc62ba52
--- /dev/null
+++ b/Geo/Pair.h
@@ -0,0 +1,24 @@
+
+#ifndef H_Pair
+#define H_Pair
+
+/** A pair of values, the types of which can be different */
+template <class L, class R>
+class Pair{
+private:
+	L	Left;
+	R	Right;
+public:
+	Pair() {}
+	Pair( const L& left, const R& right) : Left(left), Right(right) {}
+	L	left() const { return Left; }
+	void left(const L& left) { Left = left; }
+	R	right() const { return Right; }
+	void right(const R& right) { Right = right; }
+
+	L	first() const { return Left; }
+	R	second() const { return Right; }
+
+};
+
+#endif
diff --git a/Geo/Range.h b/Geo/Range.h
new file mode 100644
index 0000000000000000000000000000000000000000..227e3091f62a4a0dd6c6b1ab5b1748df432f1d87
--- /dev/null
+++ b/Geo/Range.h
@@ -0,0 +1,36 @@
+#ifndef H_Range
+#define H_Range
+
+/** represents a range of values of the template type */
+template <class T>
+class Range{
+private:
+  T	Low;
+  T	High;
+public:
+  Range() {}
+  Range( const T& low, const T& high) : Low(low), High(high) {}
+  T	low() const { return Low; }
+  void low(const T& low) { Low = low; }
+  T	high() const { return High; }
+  void high(const T& high) { High = high; }
+
+  int contains(const T& value) const;
+  int contains(const Range<T> & range) const;
+
+  int operator == (const Range<T> &range) const;
+};
+
+template<class T>
+int Range<T>::contains(const T& value) const
+{ return ( (value >= Low) && (value <= High) ); }
+
+template<class T>
+int Range<T>::contains(const Range<T>& range) const
+{ return ( (range.low() >= Low) && (range.high() <= High) ); }
+
+template<class T>
+int Range<T>::operator == (const Range<T>& range) const
+{ return ( (range.low() == Low) && (range.high() == High) ); }
+
+#endif
diff --git a/Geo/SBoundingBox3d.cpp b/Geo/SBoundingBox3d.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9b24a972d5e691b040c5fed3a11ce5d3a68be50a
--- /dev/null
+++ b/Geo/SBoundingBox3d.cpp
@@ -0,0 +1,84 @@
+#include "SBoundingBox3d.h"
+#include <float.h>
+
+SBoundingBox3d::SBoundingBox3d()
+: MinPt(DBL_MAX,DBL_MAX,DBL_MAX), MaxPt(-DBL_MAX,-DBL_MAX,-DBL_MAX)
+{}
+
+SBoundingBox3d::SBoundingBox3d(const SPoint3 & pt)
+: MinPt(pt), MaxPt(pt)
+{}
+void SBoundingBox3d::operator+=(const SPoint3 &pt)
+// note: it is possible for pt[i] to be both > MaxPt[i] and < MinPt[i]
+// the first point always will be both
+{
+  if(pt[0] < MinPt[0])
+    MinPt[0] = pt[0];
+  if (pt[0] > MaxPt[0])
+    MaxPt[0] = pt[0];
+
+  if(pt[1] < MinPt[1])
+    MinPt[1] = pt[1];
+  if (pt[1] > MaxPt[1])
+    MaxPt[1] = pt[1];
+
+  if(pt[2] < MinPt[2])
+    MinPt[2] = pt[2];
+  if (pt[2] > MaxPt[2])
+    MaxPt[2] = pt[2];
+
+}
+
+void SBoundingBox3d::operator+=(const SBoundingBox3d &box)
+{
+  (*this)+=box.MinPt;
+  (*this)+=box.MaxPt;
+}
+
+SPoint3 SBoundingBox3d::min() const
+{ return MinPt; }
+
+SPoint3 SBoundingBox3d::max() const
+{ return MaxPt; }
+
+SPoint3 SBoundingBox3d::center() const
+{ return (MinPt+MaxPt)*.5; }
+
+void SBoundingBox3d::operator*=(double scale)
+{
+  SPoint3 center = (MinPt + MaxPt)*.5;
+  MaxPt -= center;
+  MinPt -= center;
+  MaxPt *= scale;
+  MinPt *= scale;
+  MaxPt += center;
+  MinPt += center;
+
+}
+
+void SBoundingBox3d::scale(double sx, double sy, double sz)
+{
+  SPoint3 center = (MinPt + MaxPt)*.5;
+  MaxPt -= center;
+  MinPt -= center;
+  MaxPt[0] *= sx;   MaxPt[1] *= sy; MaxPt[2] *= sz;
+  MinPt[0] *= sx;   MinPt[1] *= sy; MinPt[2] *= sz;
+  MaxPt += center;
+  MinPt += center;
+}
+
+void SBoundingBox3d::makeCube()
+{
+  SPoint3 len = MaxPt-MinPt;
+  double scales[3];
+  double max=-1.0;
+
+  for(int i = 0; i < 3; i++)
+    max = len[i] > max ? len[i] : max;
+
+  for(int j = 0; j < 3; j++)
+    scales[j] = max/len[j];
+
+  scale(scales[0],scales[1],scales[2]);
+  
+}
diff --git a/Geo/SBoundingBox3d.h b/Geo/SBoundingBox3d.h
new file mode 100644
index 0000000000000000000000000000000000000000..1ee0046a5e3830e6797844273311cb3ef5c1d702
--- /dev/null
+++ b/Geo/SBoundingBox3d.h
@@ -0,0 +1,37 @@
+#ifndef H_SBoundingBox3d
+#define H_SBoundingBox3d
+
+#include "SPoint3.h"
+
+/** A bounding box class - add points and it grows to be the
+  bounding box of the point set
+  */
+class SBoundingBox3d {
+public:
+  ///
+  SBoundingBox3d();
+  ///
+  SBoundingBox3d(const SPoint3 &);
+  ///
+  void operator+=(const SPoint3 &pt);
+  ///
+  void operator+=(const SBoundingBox3d &pt);
+  ///
+  void operator*=(double scale);
+  ///
+  void scale(double, double, double);
+  ///
+  SPoint3 min() const;
+  ///
+  SPoint3 max() const;
+  ///
+  SPoint3 center() const;
+  ///
+  void makeCube();
+private:
+  SPoint3 MinPt,MaxPt;
+
+};
+
+
+#endif
diff --git a/Geo/SPoint2.h b/Geo/SPoint2.h
new file mode 100644
index 0000000000000000000000000000000000000000..f856a510187e32eeda5412c4f30e307ea44afe37
--- /dev/null
+++ b/Geo/SPoint2.h
@@ -0,0 +1,95 @@
+
+#ifndef H_SPoint2
+#define H_SPoint2
+
+#include <math.h>
+
+/** A point in 2-space */
+class SPoint2 {
+protected:
+  double P[2];
+public:
+  SPoint2() {}
+  ///
+  SPoint2(double x, double y) 
+    {P[0] = x; P[1] = y;}
+  ///
+  SPoint2(double *p)
+    {P[0] = p[0]; P[1] = p[1];}
+  SPoint2(const SPoint2 &pt)
+    {P[0] = pt.P[0]; P[1] = pt.P[1]; }
+  virtual ~SPoint2() {}
+  ///
+  void setPosition(double xx, double yy);
+  ///
+  void getPosition(double *xx, double *yy) const;
+  ///
+  void position(double *) const;
+  ///
+  inline double x(void) const;
+  ///
+  inline double y(void) const;
+  
+  ///
+  double &operator[](int);
+  ///
+  double operator[](int) const;
+  SPoint2 &operator=(const SPoint2 &p);
+
+  ///
+  void operator+=(const SPoint2 &p);
+  ///
+  void operator-=(const SPoint2 &p);
+  ///
+  void operator*=(double mult);
+  ///
+  SPoint2 operator*(double mult);
+  
+  operator double *() { return P; }
+};
+
+inline SPoint2 operator + (const SPoint2 &a, const SPoint2 &b)
+{ return SPoint2(a.x()+b.x(),a.y()+b.y()); }
+
+inline SPoint2 operator - (const SPoint2 &a, const SPoint2 &b)
+{ return SPoint2(a.x()-b.x(),a.y()-b.y()); }
+
+inline void SPoint2::setPosition(double xx, double yy)
+{ P[0] = xx;  P[1] = yy; }
+
+inline void SPoint2::getPosition(double *xx, double *yy) const
+{ *xx = P[0];  *yy = P[1]; }
+
+inline void SPoint2::position(double *p) const
+{ p[0] = P[0]; p[1] = P[1]; }
+
+inline double SPoint2::x(void) const
+{ return P[0]; }
+
+inline double SPoint2::y(void) const
+{ return P[1]; }
+
+inline SPoint2 & SPoint2::operator=(const SPoint2 &p)
+{ P[0] = p.P[0]; P[1]=p.P[1]; return *this; }
+
+inline double &SPoint2::operator[](int i)
+{ return P[i]; }
+
+inline double SPoint2::operator[](int i) const
+{ return P[i]; }
+
+inline void SPoint2::operator+=(const SPoint2 &p)
+{ P[0] += p.P[0]; P[1] += p.P[1];}
+
+inline void SPoint2::operator-=(const SPoint2 &p)
+{ P[0] -= p.P[0]; P[1] -= p.P[1];}
+
+inline void SPoint2::operator*=(double mult)
+{ P[0] *= mult; P[1] *= mult; }
+
+inline SPoint2 SPoint2::operator*(double mult)
+{ return SPoint2(P[0]*mult, P[1]*mult); }
+
+
+#endif
+
diff --git a/Geo/SPoint3.h b/Geo/SPoint3.h
new file mode 100644
index 0000000000000000000000000000000000000000..33384bd7d748e234a2a47e19415be835f8982fd7
--- /dev/null
+++ b/Geo/SPoint3.h
@@ -0,0 +1,104 @@
+
+#ifndef H_SPoint3
+#define H_SPoint3
+
+#include <iostream>
+#include <math.h>
+
+/** A point in 3-space */
+class SPoint3 {
+protected:
+  double P[3];
+public:
+  SPoint3() {}
+  ///
+  SPoint3(double x, double y, double z) 
+    {P[0] = x; P[1] = y; P[2] = z;}
+  ///
+  SPoint3(const double *p)
+    {P[0] = p[0]; P[1] = p[1]; P[2] = p[2];}
+  
+  SPoint3(const SPoint3 &pt)
+    {P[0] = pt.P[0]; P[1] = pt.P[1]; P[2] = pt.P[2]; }
+
+  virtual ~SPoint3() {}
+  ///
+  void setPosition(double xx, double yy, double zz);
+  ///
+  void getPosition(double *xx, double *yy, double *zz) const;
+  ///
+  void position(double *) const;
+  
+  ///
+  inline double x(void) const;
+  ///
+  inline double y(void) const;
+  ///
+  inline double z(void) const;
+  
+  ///
+  double &operator[](int);
+  ///
+  double operator[](int) const;
+  ///
+  SPoint3 &operator=(const SPoint3 &p);
+  ///
+  void operator+=(const SPoint3 &p);
+  ///
+  void operator-=(const SPoint3 &p);
+  ///
+  void operator*=(double mult);
+  ///
+  SPoint3 operator*(double mult);
+  ///
+  operator double *() { return P; }
+
+};
+
+inline SPoint3 operator + (const SPoint3 &a, const SPoint3 &b)
+{ return SPoint3(a.x()+b.x(),a.y()+b.y(),a.z()+b.z()); }
+
+inline SPoint3 operator - (const SPoint3 &a, const SPoint3 &b)
+{ return SPoint3(a.x()-b.x(),a.y()-b.y(),a.z()-b.z()); }
+
+inline void SPoint3::setPosition(double xx, double yy, double zz)
+{ P[0] = xx;  P[1] = yy;  P[2] = zz; }
+
+inline void SPoint3::getPosition(double *xx, double *yy, double *zz) const
+{ *xx = P[0];  *yy = P[1];  *zz = P[2]; }
+
+inline void SPoint3::position(double *p) const
+{ p[0] = P[0]; p[1] = P[1]; p[2] = P[2]; }
+
+inline double SPoint3::x(void) const
+{ return P[0]; }
+
+inline double SPoint3::y(void) const
+{ return P[1]; }
+
+inline double SPoint3::z(void) const
+{ return P[2]; }
+
+inline SPoint3 & SPoint3::operator=(const SPoint3 &p)
+{ P[0] = p.P[0]; P[1]=p.P[1]; P[2]=p.P[2]; return *this; }
+
+inline void SPoint3::operator+=(const SPoint3 &p)
+{ P[0] += p.P[0]; P[1] += p.P[1]; P[2] += p.P[2];}
+
+inline void SPoint3::operator-=(const SPoint3 &p)
+{ P[0] -= p.P[0]; P[1] -= p.P[1]; P[2] -= p.P[2];}
+
+inline void SPoint3::operator*=(double mult)
+{ P[0] *= mult; P[1] *= mult; P[2] *= mult; }
+
+inline SPoint3 SPoint3::operator*(double mult)
+{ return SPoint3(P[0]*mult, P[1]*mult, P[2] *= mult); }
+
+inline double &SPoint3::operator[](int i)
+{ return P[i]; }
+
+inline double SPoint3::operator[](int i) const
+{ return P[i]; }
+
+#endif
+
diff --git a/Geo/SVector3.cpp b/Geo/SVector3.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9b50eaf29396ff389a484da23c4861124183ae73
--- /dev/null
+++ b/Geo/SVector3.cpp
@@ -0,0 +1,99 @@
+#include "SVector3.h"
+#include "math.h"
+
+
+SVector3 cross(const SVector3 &a, const SVector3 &b)
+{
+  return SVector3(a.y()*b.z()-b.y()*a.z(),
+		 -(a.x()*b.z()-b.x()*a.z()),
+		 a.x()*b.y()-b.x()*a.y());
+}
+
+double dot(const SVector3 &a, const SVector3 &b)
+{
+  return( a.x()*b.x()+a.y()*b.y()+a.z()*b.z() );
+}
+
+double angle(const SVector3 &a, const SVector3 &b, const SVector3 &n)
+{
+  double ddot = dot(a,b)/(norm(a)*norm(b));
+  if(ddot>1.0)
+    ddot = 1.0;
+  if(ddot<-1.0)
+    ddot = -1.0;
+  double dang = acos(ddot);
+  // check if we just found the angle or 2*PI-angle
+  SVector3 check = cross(a,b);
+  // ***** check this below ******
+  if( norm(check) > 0.0000001*(norm(a)+norm(b))){  // check if a,b are colinear
+    double dir = dot(check,n);
+    dang = dir < 0 ? 2*M_PI-dang : dang;
+  } else {
+    if( ddot > 0)
+      dang = 2*M_PI;
+  }
+  return dang;
+
+}
+
+double norm(const SVector3 &v)
+{ return sqrt(dot(v,v)); }
+
+double SVector3::normalize()
+{ 
+  double n = norm(*this);
+  P[0] /= n; P[1]/= n; P[2] /= n;
+  return n;
+}
+
+SVector3 operator*(double m,const SVector3 &v)
+{
+  return SVector3(v[0]*m,v[1]*m,v[2]*m);
+}
+
+SVector3 operator*(const SVector3 &v, double m)
+{
+  return SVector3(v[0]*m,v[1]*m,v[2]*m);
+}
+
+SVector3 operator*(const SVector3 &v1, const SVector3 &v2)
+{
+  return SVector3(v1[0]*v2[0],v1[1]*v2[1],v1[2]*v2[2]);
+}
+
+SVector3 operator+(const SVector3 &a,const SVector3 &b)
+{
+  return SVector3(a[0]+b[0],a[1]+b[1],a[2]+b[2]);
+}
+
+SVector3 operator-(const SVector3 &a,const SVector3 &b)
+{
+  return SVector3(a[0]-b[0],a[1]-b[1],a[2]-b[2]);
+}
+
+SVector3 & SVector3::operator += (const SVector3 &a)
+{ 
+  for(int i = 0; i < 3; i++)
+    P[i] += a[i];
+  return *this;
+}
+
+SVector3 & SVector3::operator -= (const SVector3 &a)
+{ 
+  for(int i = 0; i < 3; i++)
+    P[i] -= a[i];
+  return *this;
+}
+
+SVector3 & SVector3::operator *= (const SVector3 &a)
+{ 
+  for(int i = 0; i < 3; i++)
+    P[i] *= a[i];
+  return *this;
+}
+
+SVector3 & SVector3::operator = (double v)
+{
+  P[0] = v; P[1] = v; P[2] = v;
+  return *this;
+}
diff --git a/Geo/SVector3.h b/Geo/SVector3.h
new file mode 100644
index 0000000000000000000000000000000000000000..531ef299d9c5f7b39335d7f353915c46892f3a28
--- /dev/null
+++ b/Geo/SVector3.h
@@ -0,0 +1,86 @@
+#ifndef H_SVector3
+#define H_SVector3
+
+#include "SPoint3.h"
+
+/* concrete class for vector of size 3 */
+class SVector3 {
+public:
+  SVector3() {}
+  /// Construct from 2 SPoints, vector from p1 to p2
+  SVector3(const SPoint3 &p1, const SPoint3 &p2)
+    : P(p2-p1) {}
+  /// Construct from a single SPoint, vector from origin to p1
+  SVector3(const SPoint3 &p1)
+    : P(p1) {}
+  SVector3(double x, double y, double z)
+    : P(x,y,z) {}
+  SVector3(double v)
+    : P(v,v,v) {}
+  SVector3(const double *array)
+    : P(array) {}
+
+  ///
+  inline double x(void) const { return P.x(); }
+  ///
+  inline double y(void) const { return P.y(); }
+  ///
+  inline double z(void) const { return P.z(); }
+
+  double normalize();
+
+  // why both [] and (), why not
+  double &operator[](int);
+  double operator[](int) const;
+  double &operator()(int);
+  double operator()(int) const;
+
+  ///
+  SVector3 & operator += (const SVector3 &);
+  ///
+  SVector3 & operator -= (const SVector3 &);
+  ///
+  SVector3 & operator *= (const SVector3 &);
+  ///
+  SVector3 & operator = (double);
+
+  operator double *() { return P; }
+  
+protected:
+  SPoint3 P;
+
+};
+
+///
+double norm(const SVector3 &v);
+///
+double dot(const SVector3 &a, const SVector3 &b);
+///
+SVector3 cross(const SVector3 &a, const SVector3 &b);
+///
+double angle(const SVector3 &a, const SVector3 &b, const SVector3 &n);
+///
+SVector3 operator*(double m,const SVector3 &v);
+///
+SVector3 operator*(const SVector3 &v, double m);
+///
+SVector3 operator*(const SVector3 &v1, const SVector3 &v2);
+///
+SVector3 operator+(const SVector3 &a,const SVector3 &b);
+///
+SVector3 operator-(const SVector3 &a,const SVector3 &b);
+
+inline double &SVector3::operator[](int i)
+{ return P[i]; }
+
+inline double SVector3::operator[](int i) const
+{ return P[i]; }
+
+inline double &SVector3::operator()(int i)
+{ return P[i]; }
+
+inline double SVector3::operator()(int i) const
+{ return P[i]; }
+
+
+#endif
diff --git a/Geo/gmshDefs.h b/Geo/gmshDefs.h
new file mode 100644
index 0000000000000000000000000000000000000000..51ae88548cacff67352512b85ebac6b26ed37ee0
--- /dev/null
+++ b/Geo/gmshDefs.h
@@ -0,0 +1,73 @@
+#ifndef _H_GMSH_DEFS_
+#define _H_GMSH_DEFS_
+#ifdef _HAVE_SGMODEL_
+#include "SGModel.h"
+#include "GVPoint.h"
+#include "GEPoint.h"
+#include "GFPoint.h"
+#include "Pair.h"
+#else
+#include "SPoint2.h"
+class GRegion;
+class GFace;
+class GEdge;
+class GVertex;
+class SGModel {
+ public:
+  SGModel ( const char * ){}
+  void add(GRegion *r){}
+  void add(GFace *f){}
+  void add(GEdge *e){}
+  void add(GVertex *v){}
+};
+
+class GeoRep {
+};
+
+struct GVPoint 
+{
+  double X,Y,Z;
+  GVPoint (double _x, double _y, double _z, const GVertex *v):X(_x),Y(_y),Z(_z){};
+};
+struct GEPoint 
+{
+  double X,Y,Z;
+  GEPoint (double _x, double _y, double _z, const GEdge *e, double par):X(_x),Y(_y),Z(_z){};
+};
+struct GFPoint 
+{
+  double X,Y,Z;
+  GFPoint (double _x, double _y, double _z,
+	   const GFace *e, const SPoint2 & par):X(_x),Y(_y),Z(_z){};
+};
+
+
+class Logical{
+public:
+  enum Value{
+    False = 0,
+    True,
+    Unknown
+    };
+};
+
+class GeomType{
+public:
+  enum Value{
+    Unknown,
+    Point,
+    Line,
+    Circle,
+    Ellipse,
+    ParametricCurve,
+    Plane,
+    Nurb,
+    Cylinder,
+    Sphere,
+    Cone,
+    Torus,
+    ParametricSurface
+    };
+};
+#endif
+#endif
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1af5aed73aff5407a03b8d19402ef962142ef9d6
--- /dev/null
+++ b/Geo/gmshEdge.cpp
@@ -0,0 +1,173 @@
+#include "gmshModel.h"
+#include "gmshEdge.h"
+#include "Interpolation.h"
+#include "CAD.h"
+#include "Geo.h"
+
+gmshEdge::gmshEdge(SGModel *model,Curve *edge,GVertex *v1,GVertex *v2)
+  : GEdge ( model, edge->Num, v1, v2 ), c(edge)
+{}
+
+gmshEdge::~gmshEdge()
+{}
+
+Range<double> gmshEdge::parBounds(int i) const
+{ 
+  return( Range<double>(  c->ubeg ,c->uend ));
+}
+
+SBoundingBox3d gmshEdge::bounds() const
+{
+  double xmin,ymin,zmin;
+  double xmax,ymax,zmax;
+  for (int i=0;i<20;i++)
+    {
+      double u = c->ubeg + (i/19.) * (c->uend - c->ubeg);
+      Vertex a = InterpolateCurve(c, u, 0);
+      if (!i)
+	{
+	  xmin = xmax = a.Pos.X;
+	  ymin = ymax = a.Pos.Y;
+	  zmin = zmax = a.Pos.Z;
+	}
+      else
+	{
+	  if(a.Pos.X < xmin) xmin = a.Pos.X;
+	  if(a.Pos.Y < ymin) ymin = a.Pos.Z;
+	  if(a.Pos.Z < zmin) zmin = a.Pos.Y;
+	  if(a.Pos.X > xmax) xmax = a.Pos.X;
+	  if(a.Pos.Y > ymax) ymax = a.Pos.Z;
+	  if(a.Pos.Z > zmax) zmax = a.Pos.Y;
+	}
+    }
+  SPoint3 bmin(xmin,ymin,zmin);
+  SPoint3 bmax(xmax,ymax,zmax);
+  SBoundingBox3d bbox(bmin);
+  bbox+=bmax;
+  return bbox;
+}
+
+GEPoint gmshEdge::point(double par) const
+{
+  Vertex a = InterpolateCurve(c, par, 0);
+  return GEPoint(a.Pos.X,a.Pos.Y,a.Pos.Z,this,par);
+}
+
+GEPoint gmshEdge::closestPoint(const SPoint3 & qp)
+{
+  Vertex v;
+  Vertex a;
+  Vertex der;
+  v.Pos.X = qp.x();
+  v.Pos.Y = qp.y();
+  v.Pos.Z = qp.z();
+  ProjectPointOnCurve (c,&v,&a,&der);
+  return GEPoint(a.Pos.X,a.Pos.Y,a.Pos.Z,this,a.u);
+}
+
+int gmshEdge::containsParam(double pt) const
+{
+  Range<double> rg = parBounds(0);
+  return (pt >= rg.low() && pt <= rg.high());
+}
+
+SVector3 gmshEdge::firstDer(double par) const
+{
+  Vertex a = InterpolateCurve(c, par, 1);
+  return SVector3(a.Pos.X,a.Pos.Y,a.Pos.Z);
+}
+
+double gmshEdge::parFromPoint(const SPoint3 &pt) const
+{
+  Vertex v;
+  Vertex a;
+  Vertex der;
+  v.Pos.X = pt.x();
+  v.Pos.Y = pt.y();
+  v.Pos.Z = pt.z();
+  ProjectPointOnCurve (c,&v,&a,&der);
+  return a.u;
+}
+
+Logical::Value gmshEdge::continuous(int) const
+{ 
+  return Logical::True;
+}
+
+Logical::Value gmshEdge::degenerate(int) const
+{ 
+  return Logical::False;
+}
+
+Logical::Value gmshEdge::periodic(int dim) const
+{
+  return Logical::False;
+}
+
+int gmshEdge::isSeam(GFace *face) const
+{
+  printf("gmshEdge::isSeam() is called.\n");
+  return 0;
+}
+
+double gmshEdge::period() const
+{
+  return 0;
+}
+
+GeomType::Value gmshEdge::geomType() const
+{
+  switch (c->Typ)
+    {
+    case MSH_SEGM_LINE : return GeomType::Line;
+    case MSH_SEGM_PARAMETRIC : return GeomType::ParametricCurve;
+    case MSH_SEGM_CIRC :  
+    case MSH_SEGM_CIRC_INV : return GeomType::Circle;
+    case MSH_SEGM_ELLI:
+    case MSH_SEGM_ELLI_INV: return GeomType::Ellipse;
+    case MSH_SEGM_BSPLN:
+    case MSH_SEGM_BEZIER: 
+    case MSH_SEGM_NURBS:
+    case MSH_SEGM_SPLN: return GeomType::Nurb; 
+    default : return GeomType::Unknown;
+    }
+}
+
+int gmshEdge::geomDirection() const
+{
+  return 1;
+}
+
+double gmshEdge::tolerance() const
+{ 
+  return 1.e-14; 
+}
+
+
+void gmshEdge::nthDerivative(double param, int n, double *array) const
+{
+  throw;
+}
+
+GVertex * gmshEdge::split(double)
+{
+  throw;
+  return 0;
+}
+
+void * gmshEdge::getNativePtr() const
+{ 
+  return c;
+}
+
+// 200306
+int gmshEdge::containsPoint(const SPoint3 &pt) const
+{ 
+  throw;
+}
+
+// 200306 
+void gmshEdge::fixPeriodicPar(double &par)
+{
+  throw;
+}
diff --git a/Geo/gmshEdge.h b/Geo/gmshEdge.h
new file mode 100644
index 0000000000000000000000000000000000000000..4fad122b002c0a7b0389b3379a50c6a5bf77154c
--- /dev/null
+++ b/Geo/gmshEdge.h
@@ -0,0 +1,59 @@
+#ifndef _H_GMSH_EDGE_
+#define _H_GMSH_EDGE_
+
+#include "gmshDefs.h"
+#include "gmshModel.h"
+#include "gmshVertex.h"
+#include "Mesh.h"
+#include "Range.h"
+
+#ifdef _HAVE_SGMODEL_
+#include "GEdge.h"
+#else
+class GEdge {
+ public:
+  GEdge (SGModel *model, int tag, GVertex *v0, GVertex *v1)
+    {};
+};
+#endif
+
+class gmshEdge : public GEdge{
+public:
+  gmshEdge(SGModel *model,Curve *edge,GVertex *v1,GVertex *v2);
+  virtual ~gmshEdge();
+  int isSeam(GFace *face) const;
+  double period() const;
+  Range<double> parBounds(int i) const;
+  virtual Logical::Value periodic(int dim=0) const;
+  virtual GeomType::Value geomType() const;
+  virtual Logical::Value degenerate(int) const;
+  virtual Logical::Value continuous(int dim) const;
+  // Geometric Ops
+  SBoundingBox3d bounds() const;
+  virtual GEPoint point(double p) const;
+  virtual GEPoint closestPoint(const SPoint3 & queryPoint);
+  virtual int containsPoint(const SPoint3 &pt) const;  
+  virtual int containsParam(double pt) const;
+  virtual SVector3 firstDer(double par) const;
+  virtual void nthDerivative(double param, int n, double *array) const;
+  virtual SPoint2 reparamOnFace(GFace * face, double epar, int dir) const{
+    throw;
+  }
+  int geomDirection() const;
+  void fixPeriodicPar(double &par);
+  virtual double tolerance() const;
+  void * getNativePtr() const; 
+  virtual GVertex * split(double par);
+  Curve *c; 
+  virtual GeoRep * geometry() {return 0;}
+#ifdef _HAVE_SGMODEL_
+  virtual SSList<GEPoint> intersect(int fAxis, double fPar, GFace *f)
+  {
+    throw;
+  }
+#endif
+  virtual double parFromPoint(const SPoint3 &pt) const;
+protected:
+};
+#endif
+
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a0ebc47efcb4dd133e591e99742c056ad0995731
--- /dev/null
+++ b/Geo/gmshFace.cpp
@@ -0,0 +1,165 @@
+#include "gmshModel.h"
+#include "gmshEdge.h"
+#include "gmshFace.h"
+#include "Interpolation.h"
+#include "CAD.h"
+#include "Geo.h"
+
+Range<double> gmshFace::parBounds(int i) const
+{ 
+/*  SPAinterval range;
+  if(i ==0)
+    range = acisSurface()->equation().param_range_u();
+  if(i==1)
+    range = acisSurface()->equation().param_range_v();
+
+  if(range.finite() )
+    return( Range<double>(range.start_pt(),range.end_pt()) );
+  else{
+    //printf("*calcParBounds()* is called\n");
+    SBoundingBox2d b = calcParBounds();
+    SPoint2 min = b.min();
+    SPoint2 max = b.max();
+    return Range<double>(min[i],max[i]);
+  }
+*/
+  return Range<double>(0, 1);
+}
+
+int gmshFace::paramDegeneracies(int dir, double *par)
+{
+  return 0;
+}
+
+SBoundingBox3d gmshFace::bounds() const
+{
+  throw;
+}
+
+SVector3 gmshFace::normal(const SPoint2 &param) const
+{
+  Vertex vu = InterpolateSurface( s, param[0], param[1],1,1);
+  Vertex vv = InterpolateSurface( s, param[0], param[1],1,2);
+  Vertex n = vu % vv;
+  n.norme();
+  return SVector3(n.Pos.X,n.Pos.Y,n.Pos.Z);
+}
+
+Pair<SVector3,SVector3> gmshFace::firstDer(const SPoint2 &param) const
+{
+  Vertex vu = InterpolateSurface( s, param[0], param[1],1,1);
+  Vertex vv = InterpolateSurface( s, param[0], param[1],1,2);
+  return Pair<SVector3,SVector3>( SVector3(vu.Pos.X,vu.Pos.Y,vu.Pos.Z),
+				  SVector3(vv.Pos.X,vv.Pos.Y,vv.Pos.Z));
+}
+
+double * gmshFace::nthDerivative(const SPoint2 &param, int n, double *array) const
+{
+  throw;
+}
+
+GFPoint gmshFace::point(const SPoint2 &pt) const
+{   
+    return point(pt.x(),pt.y()); 
+}
+
+GFPoint gmshFace::point(double par1,double par2) const
+{
+  Vertex v = InterpolateSurface( s, par1, par2,0,0);
+  return GFPoint(v.Pos.X,v.Pos.Y,v.Pos.Z,this,SPoint2(par1,par2));
+}
+
+GFPoint gmshFace::closestPoint(const SPoint3 & qp)
+{
+  Vertex v;
+
+  v.Pos.X = qp.x();
+  v.Pos.Y = qp.y();
+  v.Pos.Z = qp.z();
+
+  if ( s->Typ != MSH_SURF_PLAN )
+    ProjectPointOnSurface(s, v);
+  return GFPoint(v.Pos.X,v.Pos.Y,v.Pos.Z,
+		 this,SPoint2(v.us[0],v.us[1]));
+}
+
+int gmshFace::containsParam(const SPoint2 &pt) const
+{
+  Range<double> uu = parBounds(0);
+  Range<double> vv = parBounds(1);
+  if ((pt.x()>=uu.low() && pt.x()<=uu.high()) && (pt.y()>=vv.low() && pt.y()<=vv.high()))
+     return 1;
+  else 
+     return 0;
+}
+
+
+SPoint2 gmshFace::parFromPoint(const SPoint3 &qp) const
+{
+  Vertex v;
+  v.Pos.X = qp.x();
+  v.Pos.Y = qp.y();
+  v.Pos.Z = qp.z();
+  ProjectPointOnSurface(s, v);
+  SPoint2 pt2(v.us[0],v.us[1]); 
+  return pt2;
+}
+/*
+
+GeoRep * gmshFace::geometry()
+{ return new gmshGeoRep(this,2); }
+*/
+Logical::Value gmshFace::continuous(int dim) const
+{ 
+  return Logical::True;
+}
+
+Logical::Value gmshFace::periodic(int dim) const
+{ 
+  return Logical::False;
+}
+
+Logical::Value gmshFace::degenerate(int dim) const
+{ 
+  return Logical::False;
+}
+
+GeomType::Value gmshFace::geomType() const
+{
+  int type;
+  type = s->Typ;
+  //if(type == CONE_TYPE)
+  //return GeomType::Cone;
+  if(type == MSH_SURF_NURBS)
+    return GeomType::Nurb;
+  if(type == MSH_SURF_PLAN)
+    return GeomType::Plane;
+  return GeomType::Unknown;
+
+}
+
+int gmshFace::geomDirection() const
+{
+  return 1;
+}
+
+double gmshFace::tolerance() const
+{ return 1.e-14; }
+
+
+
+void * gmshFace::getNativePtr() const
+{ return s; }
+
+
+// added 200306
+int gmshFace::containsPoint(const SPoint3 &pt) const
+{ 
+  throw;
+}
+
+// added 200306
+double gmshFace::period(int dir) const
+{
+  return 0.;
+}
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
new file mode 100644
index 0000000000000000000000000000000000000000..887f1fe9602e16aae9258ebecfac9dd02366c34d
--- /dev/null
+++ b/Geo/gmshFace.h
@@ -0,0 +1,64 @@
+#ifndef _H_GMSH_FACE_
+#define _H_GMSH_FACE_
+
+#include "gmshDefs.h"
+#include "gmshModel.h"
+#include "gmshVertex.h"
+#include "Mesh.h"
+#include "Range.h"
+
+#ifdef _HAVE_SGMODEL_
+#include "GFace.h"
+static  SSList<GEdge*> e;
+static  SSList<int> d;
+#else
+class GFace {
+ public:
+  GFace () 
+    {};
+};
+#endif
+
+class gmshFace : public GFace 
+{
+public:
+#ifdef _HAVE_SGMODEL_
+  gmshFace(SGModel *m,Surface * face):GFace (m,face->Num,e,d), s(face){}
+#else
+  gmshFace(SGModel *m,Surface * face): s(face){}
+#endif
+   virtual ~gmshFace(){}
+   Range<double> parBounds(int i) const; 
+   virtual int paramDegeneracies(int dir, double *par); 
+
+   virtual SBoundingBox3d bounds() const; 
+   virtual GFPoint point(double par1, double par2) const; 
+   virtual GFPoint point(const SPoint2 &pt) const; 
+   virtual GFPoint closestPoint(const SPoint3 & queryPoint); 
+
+   virtual int containsPoint(const SPoint3 &pt) const;  
+   virtual int containsParam(const SPoint2 &pt) const; 
+
+   virtual SVector3 normal(const SPoint2 &param) const; 
+   virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const; 
+   virtual double * nthDerivative(const SPoint2 &param, int n,  
+ 				 double *array) const; 
+
+   virtual GeomType::Value geomType() const; 
+   virtual int geomDirection() const; 
+
+   virtual Logical::Value continuous(int dim) const; 
+   virtual Logical::Value periodic(int dim) const; 
+   virtual Logical::Value degenerate(int dim) const; 
+   virtual double period(int dir) const;	// 200306 
+   virtual double tolerance() const; 
+   virtual GeoRep * geometry() {return 0;}
+   void * getNativePtr() const; 
+   Surface *s;
+   virtual Logical::Value surfPeriodic(int dim) const
+   {throw;}
+   virtual SPoint2 parFromPoint(const SPoint3 &) const;
+protected:
+};
+
+#endif
diff --git a/Geo/gmshModel.cpp b/Geo/gmshModel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..452c385e6750c97e33ac8d5593218cf72330f828
--- /dev/null
+++ b/Geo/gmshModel.cpp
@@ -0,0 +1,118 @@
+#include "gmshModel.h"
+#include "OpenFile.h"
+#include "Tools.h"
+#include "Message.h"
+
+extern Mesh *THEM;
+GFace * gmshModel::faceByTag(int n) const
+{
+  std::list<GFace*>:: const_iterator it = faces.begin();
+  std::list<GFace*>:: const_iterator end = faces.end();
+  while (it != end)
+    {
+      gmshFace *ff = (gmshFace*) (*it);
+      if ( ff->s->Num == n)return *it;
+      ++it;
+    }
+  return 0;
+}
+
+GEdge * gmshModel::edgeByTag(int n) const
+{
+  std::list<GEdge*>:: const_iterator it = edges.begin();
+  std::list<GEdge*>:: const_iterator end = edges.end();
+  while (it != end)
+    {
+      gmshEdge *ee = (gmshEdge*) (*it);
+      if ( ee->c->Num == n)return *it;
+      ++it;
+    }
+  return 0;
+}
+GVertex * gmshModel::vertexByTag(int n) const
+{
+  std::list<GVertex*>:: const_iterator it = vertices.begin();
+  std::list<GVertex*>:: const_iterator end = vertices.end();
+  while (it != end)
+    {
+      gmshVertex *vv = (gmshVertex*) (*it);
+      if ( vv->v->Num == n)return *it;
+      ++it;
+    }
+  return 0;
+}
+gmshModel::gmshModel(char *geofile)
+  : SGModel ( "toto" )
+{
+  if (geofile)
+    {
+      OpenProblem ( geofile );
+    }
+  
+  std::set<Vertex*> points;
+
+  if(Tree_Nbr(THEM->Curves)) {
+    List_T *curves = Tree2List(THEM->Curves);
+    for(int i = 0; i < List_Nbr(curves); i++){
+      Curve *c;
+      List_Read(curves, i, &c);
+      if (c->Num >=0)
+	{
+	  if (points.find(c->beg) == points.end())
+	    {
+	      points.insert(c->beg);
+	      gmshVertex *v = new gmshVertex ( this, c->beg );
+	      vertices.push_back(v); 
+	      add(v);
+	    }
+	  if (points.find(c->end) == points.end())
+	    {
+	      points.insert(c->end);
+	      gmshVertex *v = new gmshVertex ( this , c->end );
+	      vertices.push_back(v); 
+	      add(v);
+	    }
+	  gmshEdge *e = new gmshEdge ( this, c ,
+				       vertexByTag(c->beg->Num),
+				       vertexByTag(c->end->Num) );
+	  edges.push_back(e);     	  
+	  add(e);
+	}
+    }
+    List_Delete(curves);
+  }
+  if(Tree_Nbr(THEM->Surfaces)) {
+    List_T *surfaces = Tree2List(THEM->Surfaces);
+    for(int i = 0; i < List_Nbr(surfaces); i++){
+      Surface *s;
+      List_Read(surfaces, i, &s);
+      gmshFace *f = new gmshFace ( this, s );
+      faces.push_back(f); 
+      add ( f);
+    }
+    List_Delete(surfaces);
+  } 
+  if(Tree_Nbr(THEM->Volumes)) {
+    List_T *volumes = Tree2List(THEM->Volumes);
+    for(int i = 0; i < List_Nbr(volumes); i++){
+      Volume *v;
+      List_Read(volumes, i, &v);
+      gmshRegion *r = new gmshRegion ( this, v );
+      regions.push_back(r); 
+      add ( r);
+    }
+    List_Delete(volumes);
+  }
+
+  //Msg (INFO,"gmshModel Created\n");
+  //Msg (INFO,"%d Vertices\n",vertices.size());
+  //Msg (INFO,"%d Edges   \n",edges.size());
+  //Msg (INFO,"%d Faces\n",faces.size());
+  //Msg (INFO,"%d Regions\n" ,regions.size());
+
+}
+
+SGModel *createGmshModel (char *f )
+{
+  return new gmshModel (f);
+}
diff --git a/Geo/gmshModel.h b/Geo/gmshModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..37efe5d9d77ea20e6734f84ab76387e8e3321a95
--- /dev/null
+++ b/Geo/gmshModel.h
@@ -0,0 +1,56 @@
+#ifndef _H_GMSH_MODEL_
+#define _H_GMSH_MODEL_
+
+#include <list>
+#include "gmshDefs.h"
+#include "Mesh.h"
+#include "Range.h"
+#include "Pair.h"
+#include "SPoint2.h"
+#include "SPoint3.h"
+#include "SVector3.h"
+#include "SBoundingBox3d.h"
+#include "gmshVertex.h"
+#include "gmshFace.h"
+#include "gmshEdge.h"
+#include "gmshRegion.h"
+
+
+class gmshModel : public SGModel {
+public:
+  gmshModel(char *geofile = 0); 
+  virtual ~gmshModel() {};
+  virtual void setGeomTolerance(double) {};
+  virtual double tolerance() const {return 1.e-14;}
+
+  //  virtual SBoundingBox3d bounds() const;
+  /** Find the region with the given tag. */
+  // virtual GRegion * regionByTag(int n) const;  
+  /** Find the face with the given tag. */
+  virtual GFace * faceByTag(int n) const;
+  /** Find the edge with the given tag. */
+  virtual GEdge * edgeByTag(int n) const;
+  /** Find the vertex with the given tag. */
+  virtual GVertex * vertexByTag(int n) const;
+  /** Same, but with ids */
+  //  virtual GRegion * regionByID(int n) const;
+  virtual GFace * faceByID(int n) const {return faceByTag(n);}
+  virtual GEdge * edgeByID(int n) const {return edgeByTag(n);}
+  virtual GVertex * vertexByID(int n) const {return vertexByTag(n);}
+  virtual GVertex *createVertex(std::istream &in){throw;}
+  virtual GEdge *createEdge(std::istream &in){throw;}
+  virtual GFace *createFace(std::istream &in){throw;}
+  virtual GRegion *createRegion(std::istream &in){throw;}
+#ifdef _HAVE_SGMODEL_
+  virtual SString modeler() const {throw;}
+  virtual SBoundingBox3d bounds() const {throw;}
+#endif
+protected: 
+private:
+  std::list<GFace*> faces;
+  std::list<GEdge*> edges;
+  std::list<GVertex*> vertices;
+  std::list<GRegion*> regions;
+};
+
+#endif
diff --git a/Geo/gmshRegion.h b/Geo/gmshRegion.h
new file mode 100644
index 0000000000000000000000000000000000000000..71e34f11f214b4d0467f752e7043f4dc863eafc7
--- /dev/null
+++ b/Geo/gmshRegion.h
@@ -0,0 +1,34 @@
+#ifndef _H_GMSH_REGION_
+#define _H_GMSH_REGION_
+
+#include "Mesh.h"
+#include "gmshModel.h"
+#ifdef _HAVE_SGMODEL_
+#include "GRegion.h"
+#else
+class GRegion {
+ public:
+  GRegion(SGModel *m, int tag){}
+};
+#endif
+
+class gmshRegion : public GRegion {
+public:
+  gmshRegion(SGModel *m, Volume *_v)
+    : GRegion(m, _v->Num), v(_v)
+    {
+      
+    }
+  virtual ~gmshRegion() {}
+
+  virtual GeoRep * geometry(){return 0;}
+
+  virtual double tolerance() const {return 1.e-14;}
+
+  void * getNativePtr() {return v;}
+  Volume *v;
+protected:
+
+};
+#endif
+
diff --git a/Geo/gmshVertex.h b/Geo/gmshVertex.h
new file mode 100644
index 0000000000000000000000000000000000000000..1ce3e566626a69950b825a73781a366c3d839406
--- /dev/null
+++ b/Geo/gmshVertex.h
@@ -0,0 +1,34 @@
+#ifndef _H_GMSH_VERTEX_
+#define _H_GMSH_VERTEX_
+
+#include "Mesh.h"
+#include "gmshModel.h"
+
+#ifdef _HAVE_SGMODEL_
+#include "GVertex.h"
+#else
+class GVertex {
+ public:
+  GVertex(SGModel *m, int tag){}
+};
+#endif
+
+class gmshVertex : public GVertex {
+public:
+  gmshVertex(SGModel *m, Vertex *_v) : GVertex(m, _v->Num), v(_v){}
+  virtual ~gmshVertex() {}
+
+  virtual GVPoint point() const
+  {
+    return GVPoint ( v->Pos.X,v->Pos.Y,v->Pos.Z, this);
+  }
+  virtual GeoRep * geometry() {return 0;}
+
+  virtual double tolerance() const {return 1.e-14;}
+
+  void * getNativePtr() const {return v;}
+  Vertex *v;
+ protected:
+};
+
+#endif
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 7a5b9cdcce28820dac9227f9f7e4f425f6e3ee78..362c0e40b74c92096c16495627248736b44f94c0 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.49 2006-02-25 07:22:11 geuzaine Exp $
+// $Id: BDS.cpp,v 1.50 2006-03-08 17:04:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -982,36 +982,31 @@ void recur_color_plane_surf(const double eps,
 void BDS_Mesh::color_plane_surf(double eps, int NB_T)
 {
   int current_status = 100000;
-  {
-    std::set < BDS_Triangle * >all;
-    while(1) {
-      std::list < BDS_Triangle * >plane;
-      std::list < BDS_Triangle * >::iterator it = triangles.begin();
-      std::list < BDS_Triangle * >::iterator ite = triangles.end();
-      BDS_Triangle *start = 0;
-      while(it != ite) {
-        if(all.find(*it) == all.end()) {
-          start = (BDS_Triangle *) (*it);
-        }
-        if(start)
-          break;
-        ++it;
-      }
-      if(!start)
-        break;
-      recur_color_plane_surf(eps, start, start->N(), all, plane);
-      if((int)plane.size() > NB_T) {
-        // printf("plane surface found %d triangles\n",plane.size());
-        std::list < BDS_Triangle * >::iterator xit = plane.begin();
-        std::list < BDS_Triangle * >::iterator xite = plane.end();
-        while(xit != xite) {
-          (*xit)->status = current_status;
-          ++xit;
-        }
-        current_status++;
-      }
+  std::set < BDS_Triangle * >all;
+  std::list < BDS_Triangle * >::iterator it = triangles.begin();
+  std::list < BDS_Triangle * >::iterator ite = triangles.end();
+  while(it != ite) 
+    {
+      if(all.find(*it) == all.end()) 
+	{
+	  std::list < BDS_Triangle * >plane;
+	  BDS_Triangle *start = (BDS_Triangle *) (*it);
+	  recur_color_plane_surf(eps, start, start->N(), all, plane);
+	  if((int)plane.size() > NB_T) 
+	    {
+	      // printf("plane surface found %d triangles\n",plane.size());
+	      std::list < BDS_Triangle * >::iterator xit = plane.begin();
+	      std::list < BDS_Triangle * >::iterator xite = plane.end();
+	      while(xit != xite) 
+		{
+		  (*xit)->status = current_status;
+		  ++xit;
+		}
+	      current_status++;
+	    }
+	}
+      ++it;
     }
-  }
 }
 
 bool BDS_Mesh::extractVolumes()
@@ -1069,6 +1064,9 @@ void BDS_Mesh::classify(double angle, int NB_T)
     }
     geom.clear();
   }
+
+  color_plane_surf(3.1415/200, 40);
+
   {
     std::list < BDS_Edge * >::iterator it = edges.begin();
     std::list < BDS_Edge * >::iterator ite = edges.end();
@@ -1278,6 +1276,7 @@ void BDS_Mesh::classify(double angle, int NB_T)
   Msg(INFO, "Computing curvatures");
   compute_curvatures(edges);
   Msg(INFO, "Reverse engineering surfaces");
+
   reverseEngineerCAD();
   Msg(INFO, "Creating search structures");
   createSearchStructures();
@@ -1640,6 +1639,27 @@ bool BDS_Mesh::read_mesh(const char *filename)
             Max[2] = (Max[2] > z) ? Max[2] : z;
             add_point(i + 1, x, y, z);
           }
+
+	  // NORMALIZE HERE
+
+	  LC = sqrt((Min[0] - Max[0]) * (Min[0] - Max[0]) +
+		    (Min[1] - Max[1]) * (Min[1] - Max[1]) +
+		    (Min[2] - Max[2]) * (Min[2] - Max[2]));
+	  	  
+	  std::set < BDS_Point *, PointLessThan >::iterator it = points.begin();
+	  std::set < BDS_Point *, PointLessThan >::iterator ite = points.end();
+	  
+	  while(it != ite) 
+	    {
+	      (*it)->X /= LC;
+	      (*it)->Y /= LC;
+	      (*it)->Z /= LC;
+	      ++it;
+	    }
+	  Min[0]/=LC;Min[1]/=LC;Min[2]/=LC;
+	  Max[0]/=LC;Max[1]/=LC;Max[2]/=LC;
+	  LC = 1;
+
           MAXPOINTNUMBER = nbv + 1;
         }
         else if(!strcmp(name, "Triangles")) {
@@ -1671,9 +1691,8 @@ bool BDS_Mesh::read_mesh(const char *filename)
       }
     }
 
-    LC = sqrt((Min[0] - Max[0]) * (Min[0] - Max[0]) +
-              (Min[1] - Max[1]) * (Min[1] - Max[1]) +
-              (Min[2] - Max[2]) * (Min[2] - Max[2]));
+
+
   }
   else {
     throw;
@@ -2378,7 +2397,7 @@ void BDS_Mesh::compute_metric_edge_lengths(const BDS_Metric & metric)
 	BDS_GeomEntity *g = e->g;
 	if (g && g->classif_degree == 1)
 	  {
-	    double l = (*it)->length() * 4;
+	    double l = (*it)->length() * 1000;
 	    if (l < (*it)->target_length)
 	      (*it)->target_length =l;
 	  }   
@@ -2387,6 +2406,8 @@ void BDS_Mesh::compute_metric_edge_lengths(const BDS_Metric & metric)
     }
   }
 
+  
+
   {
     std::list < BDS_Edge * >::iterator it = edges.begin();
     std::list < BDS_Edge * >::iterator ite = edges.end();
@@ -2421,29 +2442,37 @@ void BDS_Mesh::compute_metric_edge_lengths(const BDS_Metric & metric)
   }
   //  printf("smoothing\n");
   const int NITER = 3;
+  const double BETA = metric.beta;
+
+  std::vector<BDS_Point *> temp(points.size());
+  std::copy( points.begin(),points.end(),temp.begin());
 
   for(int I = 0; I < NITER; ++I) {
-    const double BETA = metric.beta;
-    std::set < BDS_Point *, PointLessThan >::iterator it = points.begin();
-    std::set < BDS_Point *, PointLessThan >::iterator ite = points.end();
+
+    std::vector < BDS_Point * >::const_iterator it = temp.begin();
+    std::vector < BDS_Point * >::const_iterator ite = temp.end();
+
     while(it != ite) {
 
       std::list < BDS_Edge * >::iterator eit = (*it)->edges.begin();
+      std::list < BDS_Edge * >::iterator eite = (*it)->edges.end();
 
-      double l_min = metric._max;
+      double l_min = metric._max;      
 
-      while(eit != (*it)->edges.end()) {
-        if(l_min > (*eit)->target_length)
-          l_min = (*eit)->target_length;
+      while(eit != eite) {
+	BDS_Edge *ee = (*eit);
+        if(l_min > ee->target_length)
+          l_min = ee->target_length;
         ++eit;
       }
 
       l_min /= BETA;
 
       eit = (*it)->edges.begin();
-      while(eit != (*it)->edges.end()) {
-        if((*eit)->target_length > l_min)
-          (*eit)->target_length = l_min;
+      while(eit != eite) {
+	BDS_Edge *ee = (*eit);
+        if(ee->target_length > l_min)
+          ee->target_length = l_min;
         ++eit;
       }
       ++it;
@@ -2541,7 +2570,13 @@ int BDS_Mesh::adapt_mesh(const BDS_Metric & metric, bool smooth,
 	    double qb2 = quality_triangle ( (*it)->p2 , op[0] , op[1] );
 	    double qa = (qa1<qa2)?qa1:qa2; 
 	    double qb = (qb1<qb2)?qb1:qb2; 
-	    if (qb > qa)	      
+
+	    double dd = dist_droites_gauches((*it)->p1 , (*it)->p2,
+					     op[0],op[1]);
+	    
+	    double ll = (*it)->length();
+
+	    if ((qb > qa && dd < 0.1 * ll) || (qb > 5 * qa))	      
 	      {
 		nb_modif++;
 		swap_edge ( *it );
@@ -2554,18 +2589,22 @@ int BDS_Mesh::adapt_mesh(const BDS_Metric & metric, bool smooth,
   cleanup();  
   if (smooth ){
     Msg(INFO,"smoothing %d points\n",points.size());
-    std::set<BDS_Point*, PointLessThan>::iterator it   = points.begin();
-    std::set<BDS_Point*, PointLessThan>::iterator ite  = points.end();
-    while (it != ite)
+
+    std::vector <BDS_Point *> temp_l(points.size());
+    std::copy( points.begin(),points.end(),temp_l.begin());
+    std::vector < BDS_Point * >::iterator itx = temp_l.begin();
+    std::vector < BDS_Point * >::iterator itxe = temp_l.end();
+    while (itx != itxe)
       {
-	smooth_point(*it,geom_mesh);
-	++it;
+	smooth_point(*itx,geom_mesh);
+	++itx;
       }
+    printf("coucouc1\n");
   }
   
   Msg(INFO,"%d snaps have succeeded , %d have failed\n",SNAP_SUCCESS,SNAP_FAILURE);
   // outputScalarField (triangles,"b.pos");
-  // applyOptimizationPatterns(); // FIXME: this is buggy
+   applyOptimizationPatterns(); // FIXME: this is buggy
   cleanup();  
   return nb_modif;
 }
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 5d48497f04405c8b2201004feb4e8d667e1ca28f..32d45adb3dacdeba33ce37b09d8c928c4607277f 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -379,6 +379,17 @@ public:
 		      (n[0]->Y+n[1]->Y+n[2]->Y)/3.,
 		      (n[0]->Z+n[1]->Z+n[2]->Z)/3.);
   }
+
+  inline double longest_edge_length() const
+  {
+    double l1 = e1->length();
+    double l2 = e2->length();
+    double l3 = e3->length();
+
+    if (l1  > l2 && l1 > l3) return l1;
+    if (l2  > l1 && l2 > l3) return l2;
+    return l3;
+  }
   
   inline void _update()
   { 
@@ -608,20 +619,26 @@ public:
   std::list<BDS_Edge*>      edges; 
   std::list<BDS_Triangle*>   triangles; 
   std::list<BDS_Tet*>        tets; 
+  // Points
   BDS_Point * add_point(int num , double x, double y,double z);
-  BDS_Edge  * add_edge(int p1, int p2);
   void del_point(BDS_Point *p);
+  BDS_Point *find_point(int num);
+  // Edges
+  BDS_Edge  * add_edge(int p1, int p2);
   void del_edge(BDS_Edge *e);
+  BDS_Edge  *find_edge(int p1, int p2);
+  BDS_Edge  *find_edge(BDS_Point *p1, BDS_Point *p2, BDS_Triangle *t)const;
+  // Triangles
   BDS_Triangle *add_triangle(int p1, int p2, int p3); 
   BDS_Triangle *add_triangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3);
-  BDS_Tet *add_tet(int p1, int p2, int p3,int p4);
   void del_triangle(BDS_Triangle *t);
-  void add_geom(int degree, int tag);
-  BDS_Point *find_point(int num);
-  BDS_Edge  *find_edge(int p1, int p2);
   BDS_Triangle  *find_triangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3);
-  BDS_Edge  *find_edge(BDS_Point *p1, BDS_Point *p2, BDS_Triangle *t)const;
+  // Tets
+  BDS_Tet *add_tet(int p1, int p2, int p3,int p4);
+  // Geom entities
+  void add_geom(int degree, int tag);
   BDS_GeomEntity *get_geom(int p1, int p2);
+  // 2D operators
   bool swap_edge(BDS_Edge *);
   bool collapse_edge(BDS_Edge *, BDS_Point*, const double eps);
   void snap_point(BDS_Point* , BDS_Mesh *geom = 0);
@@ -629,6 +646,7 @@ public:
   bool smooth_point_b(BDS_Point*); 
   bool move_point(BDS_Point *p , double X, double Y, double Z);
   bool split_edge(BDS_Edge *, double coord, BDS_Mesh *geom = 0);
+  // Global operators
   void classify(double angle, int nb = -1); 
   void color_plane_surf(double eps , int nb);
   void reverseEngineerCAD() ;
diff --git a/Mesh/DiscreteSurface.cpp b/Mesh/DiscreteSurface.cpp
index 1fd66bd64b9e5a931cbdd5c77b55e1cb42827690..2e62d1ecc0bae63aef55184ae4e8f794656b73a0 100644
--- a/Mesh/DiscreteSurface.cpp
+++ b/Mesh/DiscreteSurface.cpp
@@ -1,4 +1,4 @@
-// $Id: DiscreteSurface.cpp,v 1.37 2006-02-26 16:26:09 geuzaine Exp $
+// $Id: DiscreteSurface.cpp,v 1.38 2006-03-08 17:04:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -224,9 +224,17 @@ void BDS_To_Mesh_2(Mesh * m)
       m->bds_mesh->points.begin();
     std::set < BDS_Point *, PointLessThan >::iterator ite =
       m->bds_mesh->points.end();
+
+
     while(it != ite) {
+      //      double dx = 1.e-3 * m->bds_mesh->LC * (double)rand() / (double)RAND_MAX;
+      //      double dy = 1.e-3 * m->bds_mesh->LC * (double)rand() / (double)RAND_MAX;
+      //      double dz = 1.e-3 * m->bds_mesh->LC * (double)rand() / (double)RAND_MAX;
+      double dx = 0;
+      double dy = 0;
+      double dz = 0;
       Vertex *vert =
-        Create_Vertex((*it)->iD, (*it)->X, (*it)->Y, (*it)->Z, (*it)->min_edge_length(), 0.0);
+        Create_Vertex((*it)->iD, (*it)->X+dx, (*it)->Y+dy, (*it)->Z+dz, (*it)->min_edge_length(), 0.0);
       Tree_Add(m->Vertices, &vert);
       ++it;
     }
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index fdd238bdd479e37797e4c299782817c42be33e49..0ae02dca4869dfc87a13c90eb66b45c32a23d105 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.81 2006-02-26 16:26:09 geuzaine Exp $
+// $Id: Generator.cpp,v 1.82 2006-03-08 17:04:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -250,7 +250,6 @@ bool TooManyElements(Mesh *M, int dim){
   Tree_Action(M->Points, GetSumOfAllLc);
   SumOfAllLc /= (double)Tree_Nbr(M->Points);
   if(pow(CTX.lc / SumOfAllLc, dim) < 1.e7) return false;
-
   return !GetBinaryAnswer("Your choice of characteristic lengths will likely produce\n"
 			  "a very large mesh. Do you really want to continue?\n\n"
 			  "(To disable this warning in the future, select `Enable\n"
diff --git a/Mesh/Read_Mesh.cpp b/Mesh/Read_Mesh.cpp
index ea7418abd6c8db5c80446a18989eca235d1bfae5..518e56ad63e87458962098ec807e78d78c4fd25e 100644
--- a/Mesh/Read_Mesh.cpp
+++ b/Mesh/Read_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Read_Mesh.cpp,v 1.98 2006-01-17 17:09:05 geuzaine Exp $
+// $Id: Read_Mesh.cpp,v 1.99 2006-03-08 17:04:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -19,6 +19,7 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
+#include <map>
 #include "Gmsh.h"
 #include "Geo.h"
 #include "CAD.h"
@@ -162,6 +163,34 @@ double SetLC(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4 = 0)
   return lc;
 }
 
+int getPartition ( const std::multimap<int,int> &nod2proc , int nbNod, Vertex *verts )
+{
+
+  std::map<int,int> proc2nbnod;
+  for (int i=0;i<nbNod;i++)
+    {
+      std::multimap<int,int>::const_iterator beg = nod2proc.lower_bound(verts[i].Num);      
+      std::multimap<int,int>::const_iterator end = nod2proc.upper_bound(verts[i].Num);
+      while ( beg != end )
+	{
+	  int iProc = (*beg).second;	  
+	  if (proc2nbnod.find(iProc) == proc2nbnod.end()) proc2nbnod[iProc] = 1;
+	  else proc2nbnod[iProc] = proc2nbnod[iProc]+1;
+	  ++beg;
+	}
+    }
+  {
+    std::map<int,int>::const_iterator beg = proc2nbnod.begin();      
+    std::map<int,int>::const_iterator end = proc2nbnod.end();
+    while (beg!=end)
+      {
+	if ( (*beg).second == nbNod ) return (*beg).first;
+	beg++;
+      }    
+  }
+  return 0;  
+}
+
 void Read_Mesh_MSH(Mesh * M, FILE * fp)
 {
   char String[256];
@@ -183,6 +212,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
   Tree_T *Duplicates = NULL;
   Tree_T *groups = List2Tree(M->PhysicalGroups, comparePhysicalGroup);
 
+  std::multimap<int,int> nod2proc;
   while(1) {
     do {
       if(!fgets(String, sizeof(String), fp))
@@ -194,6 +224,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
     if(feof(fp))
       break;
 
+    Msg(INFO, "%s\n", &String[1]);
     /*  F o r m a t  */
 
     if(!strncmp(&String[1], "MeshFormat", 10)) {
@@ -236,6 +267,23 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
       }
     }
 
+    /*  NODE'S PROCESSORS */
+
+    else if(!strncmp(&String[1], "PARA", 4))
+      {
+	fscanf(fp, "%d", &Nbr_Nodes);
+	Msg(INFO,"%d parallel nodes\n",Nbr_Nodes);
+	for(i_Node = 0; i_Node < Nbr_Nodes; i_Node++) {
+	  int nbProc;
+	  fscanf(fp, "%d %d",&Num, &nbProc);
+	  for (int iProc=0;iProc<nbProc;iProc++)
+	    {
+	      int iProcNum;
+	      fscanf(fp, "%d",&iProcNum);
+	      nod2proc.insert(std::pair<int, int>(Num,iProcNum ));
+	    }
+	}
+      }
     /*  NODES  */
 
     else if(!strncmp(&String[1], "NOD", 3) ||
@@ -323,8 +371,15 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
 	}
 
         for(j = 0; j < Nbr_Nodes; j++)
-          fscanf(fp, "%d", &verts[j].Num);
-	
+	  {
+	    fscanf(fp, "%d", &verts[j].Num);
+	  }
+
+	if (nod2proc.size())
+	  {
+	    Partition = getPartition ( nod2proc ,Nbr_Nodes , verts );
+	  }
+
         for(i = 0; i < Nbr_Nodes; i++) {
           vertsp[i] = &verts[i];
           if(!(vertspp = (Vertex **) Tree_PQuery(M->Vertices, &vertsp[i])))
diff --git a/Mesh/Vertex.h b/Mesh/Vertex.h
index da09302928ab7dbefa6b09e6008f0b3e6a64a0a1..3a98170caa8f2f52a304ba125dbf4718c233223f 100644
--- a/Mesh/Vertex.h
+++ b/Mesh/Vertex.h
@@ -22,9 +22,9 @@
 
 #include "List.h"
 
-typedef struct {
+struct Coord{
   double X,Y,Z;
-}Coord;
+};
 
 class Vertex {
   public :
diff --git a/configure.in b/configure.in
index e0e0533ca6cbc4cf8bd0aa5620d8df0ece6c2a98..166a89eac7cbfbf41b5816eea5428256a471d89f 100644
--- a/configure.in
+++ b/configure.in
@@ -1,4 +1,4 @@
-dnl $Id: configure.in,v 1.93 2006-02-25 15:59:07 geuzaine Exp $
+dnl $Id: configure.in,v 1.94 2006-03-08 17:03:35 remacle Exp $
 dnl
 dnl Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 dnl
@@ -249,8 +249,8 @@ if test "x$enable_gui" != "xno"; then
 
 else
 
-  GMSH_DIRS="${GMSH_DIRS} Box"
-  GMSH_LIBS="-Llib -lGmshBox -lGmshParser -lGmshMesh -lGmshGeo -lGmshCommon"
+  GMSH_DIRS="${GMSH_DIRS} Box BoxMain"
+  GMSH_LIBS="-Llib -lGmshBoxMain -lGmshBox -lGmshParser -lGmshMesh -lGmshGeo -lGmshCommon"
   GMSH_LIBS="${GMSH_LIBS} -lGmshDataStr -lGmshPlugin -lGmshNumeric -lGmshParallel"
 
 fi
@@ -332,6 +332,30 @@ if test "x$enable_contrib" != "xno"; then
     fi
   fi
 
+
+  dnl Check for MODEL
+  AC_CHECK_FILE(${DEVROOT}/MeshAdapt/model/model/model/SGModel.h, MODEL="yes", MODEL="no")
+  if test "x${MODEL}" = "xyes"; then
+    if test "x$enable_model" != "xno"; then
+       GMSH_LIBS="${GMSH_LIBS} ${DEVROOT}/MeshAdapt/model/lib/x86_linux/libmodel-O.a"
+       FLAGS="-D_HAVE_SGMODEL_ ${FLAGS}"
+       echo "********************************************************************"
+       echo "You are building a version of Gmsh that contains model, the"
+       echo "Modeler Interface."
+       echo "Please note that by doing so, you agree with model's licensing"
+       echo "requirements."
+       echo "To disable model, run configure again with the --disable-model"
+       echo "option."
+       echo "********************************************************************"
+    fi
+  else
+    if test "x$enable_model" != "xno"; then
+       echo "********************************************************************"
+       echo "If you want to use model for doing modeling interface, contact us."
+       echo "********************************************************************"
+    fi
+  fi
+
   dnl Check for METIS
   AC_CHECK_FILE(./contrib/Metis/metis.h, METIS="yes", METIS="no")
   if test "x${METIS}" = "xyes"; then