diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d5a9cfd7f681e12efd35336ca642942dc84e0121
--- /dev/null
+++ b/Geo/GEdge.cpp
@@ -0,0 +1,10 @@
+#include "GEdge.h"
+#include <algorithm>
+void GEdge::addFace ( GFace *e )
+{ 
+  l_faces.push_back (e);  
+}
+void GEdge::delFace ( GFace *e )
+{ 
+  l_faces.erase(std::find(l_faces.begin(),l_faces.end(),e));  
+}
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 370450f4ed24d6c0603bc4703c8e50741c935563..96f20d56d4a7b88b4e556ac5072fe4d3313a67fd 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -3,6 +3,9 @@
 
 #include "GEntity.h"
 #include "GVertex.h"
+#include "SVector3.h"
+#include "SPoint3.h"
+#include "SPoint2.h"
 
 /** A model edge. */
 
@@ -14,22 +17,17 @@ public:
 	GVertex *_v1)
     : GEntity (model,tag),v0(_v0),v1(_v1) 
     {
-      v0->add_edge (this);
-      v1->add_edge (this);
+      v0->addEdge (this);
+      v1->addEdge (this);
     }
   virtual ~GEdge() 
     {
-      v0->del_edge (this);
-      v1->del_edge (this);
+      v0->delEdge (this);
+      v1->delEdge (this);
     }
 
   virtual int dim() const {return 1;}
 
-  virtual std::list<GRegion*> regions() const;
-  virtual std::list<GFace*> faces() const;
-  virtual std::list<GEdge*> edges() const;
-  virtual std::list<GVertex*> vertices() const;
-
   virtual bool periodic(int dim=0) const = 0;
   virtual bool continuous(int dim=0) const = 0;
 
@@ -43,10 +41,10 @@ public:
   virtual double parFromPoint(const SPoint3 &) const = 0;
 
   /** Get the point for the given parameter location. */
-  virtual GEPoint point(double p) const = 0;  
+  virtual GPoint point(double p) const = 0;  
 
   /** Get the closest point on the edge to the given point. */
-  virtual GEPoint closestPoint(const SPoint3 & queryPoint);
+  virtual GPoint closestPoint(const SPoint3 & queryPoint);
 
   /** True if the edge contains the given parameter. */
   virtual int containsParam(double pt) const = 0;
@@ -60,15 +58,14 @@ public:
   /** reparmaterize the point onto the given face. */
   virtual SPoint2 reparamOnFace(GFace *face, double epar,int dir) const = 0;			  
 
-  void addFace ( GFace *f ){ l_faces.push_back (f);  }
-  void delFace ( GFace *f ){ l_faces.erase(std::find(l_faces.begin(),l_faces.end(),f));  }
+  void addFace ( GFace *f );
+  void delFace ( GFace *f );
 
 
 protected:
 
-  GVertex v0,v1;
+  GVertex *v0,*v1;
   std::list<GFace *> l_faces;
-
 };
 
 
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2c0306d1d5fe1da7c6284555495ab9409a0d1d4c
--- /dev/null
+++ b/Geo/GFace.cpp
@@ -0,0 +1,13 @@
+#include "GFace.h"
+#include "GEdge.h"
+
+GFace::~GFace ()
+{ 
+  std::list<GEdge*>::iterator it = l_edges.begin();
+
+  while (it != l_edges.end())
+    {
+      (*it)->delFace(this);
+      ++it;
+    }
+}
diff --git a/Geo/GFace.h b/Geo/GFace.h
new file mode 100644
index 0000000000000000000000000000000000000000..4913ca6de011c4ae2c64e6c26b8a30e0edc2011f
--- /dev/null
+++ b/Geo/GFace.h
@@ -0,0 +1,67 @@
+#ifndef H_GFace
+#define H_GFace
+
+class GRegion;
+#include "GPoint.h"
+#include "GEntity.h"
+#include "SPoint2.h"
+#include "SVector3.h"
+#include "Pair.h"
+
+/** A model face. 
+ */
+class GFace : public GEntity 
+{
+protected: 
+  std::list<GEdge *> l_edges;
+  std::list<int>     l_dirs;
+  GRegion *r1, *r2;
+public:
+  GFace(GModel *model, int tag) : GEntity (model,tag),r1(0),r2(0){}
+  virtual ~GFace();
+
+  void addRegion ( GRegion *r ){ r1?r2=r:r1=r;  }
+  void delRegion ( GRegion *r ){ if(r1==r)r1=r2;r2=0;  }
+
+  virtual int dim() const {return 2;}
+
+  /** Get the location of any parametric degeneracies on the face in the
+    given parametric direction. */
+  virtual int paramDegeneracies(int dir, double *par) = 0;
+
+  /** Return the point on the face corresponding to the given parameter. */
+  virtual GPoint point(double par1, double par2) const = 0;
+  virtual GPoint point(const SPoint2 &pt) const = 0;
+
+  /** Return the parmater location on the face given a point in space that
+    is on the face. */
+  virtual SPoint2 parFromPoint(const SPoint3 &) const = 0;
+
+  /** True if the parameter value is interior to the face. */
+  virtual int containsParam(const SPoint2 &pt) const = 0;
+
+  /** Period of the face in the given direction. */
+  virtual double period(int dir) const;
+
+  /** Return the point on the face closest to the given point. */
+  virtual GPoint closestPoint(const SPoint3 & queryPoint);
+
+  /** Return the normal to the face at the given parameter location. */
+  virtual SVector3 normal(const SPoint2 &param) const = 0;
+
+  /** Return the first derivate of the face at the parameter location. */
+  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const = 0;
+  
+  /** Return the nth derivate of the face at the parametric location. */
+  virtual double * nthDerivative(const SPoint2 &param, int n, 
+				 double *array) const;
+  /* true if the surface underlying the face is periodic and we
+     need to worry about that. */
+  virtual bool surfPeriodic(int dim) const = 0;
+protected:
+};
+
+
+#endif
+
+
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7791e61595c2af22f42278f9f2ab2d4345514f47
--- /dev/null
+++ b/Geo/GModel.cpp
@@ -0,0 +1,70 @@
+#include "GModel.h"
+
+int GModel::numRegion() const
+{
+  return regions.size();
+}
+int GModel::numFace  () const
+{
+  return faces.size();
+}
+int GModel::numEdge  () const
+{
+  return edges.size();
+}
+int GModel::numVertex() const
+{
+  return vertices.size();
+}
+
+GRegion * GModel::regionByTag(int n) const
+{
+  std::list<GRegion*>:: const_iterator it = regions.begin();
+  std::list<GRegion*>:: const_iterator end = regions.end();
+  while (it != end)
+    {
+      GRegion *rr = (GRegion*) (*it);
+      if ( rr->tag() == n)return *it;
+      ++it;
+    }
+  return 0;
+}
+
+GFace * GModel::faceByTag(int n) const
+{
+  std::list<GFace*>:: const_iterator it = faces.begin();
+  std::list<GFace*>:: const_iterator end = faces.end();
+  while (it != end)
+    {
+      GFace *ff = (GFace*) (*it);
+      if ( ff->tag() == n)return *it;
+      ++it;
+    }
+  return 0;
+}
+
+GEdge * GModel::edgeByTag(int n) const
+{
+  std::list<GEdge*>:: const_iterator it = edges.begin();
+  std::list<GEdge*>:: const_iterator end = edges.end();
+  while (it != end)
+    {
+      GEdge *ee = (GEdge*) (*it);
+      if ( ee->tag() == n)return *it;
+      ++it;
+    }
+  return 0;
+}
+GVertex * GModel::vertexByTag(int n) const
+{
+  std::list<GVertex*>:: const_iterator it = vertices.begin();
+  std::list<GVertex*>:: const_iterator end = vertices.end();
+  while (it != end)
+    {
+      GVertex *vv = (GVertex*) (*it);
+      if ( vv->tag() == n)return *it;
+      ++it;
+    }
+  return 0;
+}
+
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 4d09093d56a0caf448ba5d31525949bc129b9934..1105b575b7c7d9ceb18a5c8d667757ae997c8a98 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -1,6 +1,7 @@
 #ifndef H_GModel
 #define H_GModel
 
+#include <algorithm>
 #include "GVertex.h"
 #include "GEdge.h"
 #include "GFace.h"
@@ -13,10 +14,10 @@ class GModel
 public:
   virtual ~GModel();
 
-  typedef std::list<GRegion*>::const_iterator riter;
-  typedef std::list<GFace*>::const_iterator   fiter;
-  typedef std::list<GEdge*>::const_iterator   eiter;
-  typedef std::list<GVertex*>::const_iterator viter;
+  typedef std::list<GRegion*>::iterator riter;
+  typedef std::list<GFace*>::iterator   fiter;
+  typedef std::list<GEdge*>::iterator   eiter;
+  typedef std::list<GVertex*>::iterator viter;
 
   /** Returns the geometric tolerance for the entire model. */
   virtual double tolerance() const =0;
@@ -28,39 +29,39 @@ public:
   int numVertex() const;
 
   /** Get the nth region in this model. */
-  GRegion * region(int n) const;
-  GFace   * face  (int n) const;
-  GEdge   * edge  (int n) const;
-  GVertex * vertex(int n) const;
+  //  GRegion * region(int n) const;
+  //  GFace   * face  (int n) const;
+  //  GEdge   * edge  (int n) const;
+  //  GVertex * vertex(int n) const;
 
   /** Get an iterator initialized to the first entity in this model. */
-  riter firstRegion() const {return regions.begin();}
-  fiter firstFace() const {return faces.begin();}
-  eiter firstEdge() const {return edges.begin();}
-  viter firstVertex() const {return vertices.begin();}
-  riter lastRegion() const {return regions.end();}
-  fiter lastFace() const {return faces.end();}
-  eiter lastEdge() const {return edges.end();}
-  viter lastVertex() const {return vertices.end();}
+  riter firstRegion() {return regions.begin();}
+  fiter firstFace() {return faces.begin();}
+  eiter firstEdge()  {return edges.begin();}
+  viter firstVertex() {return vertices.begin();}
+  riter lastRegion() {return regions.end();}
+  fiter lastFace() {return faces.end();}
+  eiter lastEdge() {return edges.end();}
+  viter lastVertex() {return vertices.end();}
 
   /** Find the region with the given tag. */
   virtual GRegion * regionByTag(int n) const;
-  virtual GFace * faceByTag(int n) const;
-  virtual GEdge * edgeByTag(int n) const;
+  virtual GFace   * faceByTag  (int n) const;
+  virtual GEdge   * edgeByTag  (int n) const;
   virtual GVertex * vertexByTag(int n) const;
 
-  virtual GRegion * regionByID(int n) const;
-  virtual GFace * faceByID(int n) const;
-  virtual GEdge * edgeByID(int n) const;
-  virtual GVertex * vertexByID(int n) const;
+  //  virtual GRegion * regionByID(int n) const;
+  //  virtual GFace   * faceByID  (int n) const;
+  //  virtual GEdge   * edgeByID  (int n) const;
+  //  virtual GVertex * vertexByID(int n) const;
 
   virtual void setGeomTolerance(double) {};
-  void setDisplayCoherence(int ); // default is coherent
+  //  void setDisplayCoherence(int ); // default is coherent
 
   void add(GRegion *r){regions.push_back(r);}
-  void add(GFace *f)  {regions.push_back(f);}
-  void add(GEdge *e)  {regions.push_back(e);}
-  void add(GVertex *v){regions.push_back(v);}
+  void add(GFace *f)  {faces.push_back(f);}
+  void add(GEdge *e)  {edges.push_back(e);}
+  void add(GVertex *v){vertices.push_back(v);}
 
   void remove(GRegion *r){regions.erase(std::find(firstRegion(),lastRegion(),r));}
   void remove(GFace *f){faces.erase(std::find(firstFace(),lastFace(),f));}
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c478abd5a6e73516647a53de3b88437bd16e572b
--- /dev/null
+++ b/Geo/GRegion.cpp
@@ -0,0 +1,13 @@
+#include "GRegion.h"
+#include "GFace.h"
+
+GRegion::~GRegion ()
+{ 
+  std::list<GFace*>::iterator it = l_faces.begin();
+
+  while (it != l_faces.end())
+    {
+      (*it)->delRegion(this);
+      ++it;
+    }
+}
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
new file mode 100644
index 0000000000000000000000000000000000000000..8e7f06ae1d48d4cef51bf398d57cd51d79b9757d
--- /dev/null
+++ b/Geo/GRegion.h
@@ -0,0 +1,17 @@
+#ifndef H_GRegion
+#define H_GRegion
+
+#include "GEntity.h"
+
+/** A model region. */
+class GRegion : public GEntity {
+  std::list<GFace *> l_faces;
+  std::list<int>     l_dirs;
+public:
+  GRegion(GModel *model, int tag) : GEntity (model,tag) {}
+  virtual ~GRegion();
+  virtual int dim() const {return 3;}
+};
+
+#endif
+
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..aa6526e6d1d9e8fce35faec979457ef92bdb8227
--- /dev/null
+++ b/Geo/GVertex.cpp
@@ -0,0 +1,10 @@
+#include "GVertex.h"
+#include <algorithm>
+void GVertex::addEdge ( GEdge *e )
+{ 
+  l_edges.push_back (e);  
+}
+void GVertex::delEdge ( GEdge *e )
+{ 
+  l_edges.erase(std::find(l_edges.begin(),l_edges.end(),e));  
+}
diff --git a/Geo/GVertex.h b/Geo/GVertex.h
new file mode 100644
index 0000000000000000000000000000000000000000..a04de813884b5a99972f89bd22ff6f144d1e3dfd
--- /dev/null
+++ b/Geo/GVertex.h
@@ -0,0 +1,21 @@
+#ifndef H_GVertex
+#define H_GVertex
+
+#include "GEntity.h"
+#include "GPoint.h"
+
+/** A model vertex. */
+class GVertex  : public GEntity 
+{
+public:
+  GVertex(GModel *m, int tag) : GEntity (m,tag) {}
+  virtual ~GVertex() {}
+  virtual GPoint point() const = 0;
+  void addEdge ( GEdge *e );
+  void delEdge ( GEdge *e );
+protected:
+  std::list<GEdge*> l_edges;
+};
+
+#endif
+