diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 516b538e41759911df608ac4cb4755c74af4c99e..e9c242a829cbf4aa0cd9a05f3ea1e59a254c6fd6 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -3,9 +3,12 @@
 #include "Interpolation.h"
 #include "CAD.h"
 #include "Geo.h"
+#include "Mesh.h"
+#include "Create.h"
 #include "Context.h"
 
 extern Context_T CTX;
+extern Mesh *THEM;
 
 gmshEdge::gmshEdge(GModel *model, Curve *edge, GVertex *v1, GVertex *v2)
   : GEdge(model, edge->Num, v1, v2), c(edge)
@@ -13,8 +16,11 @@ gmshEdge::gmshEdge(GModel *model, Curve *edge, GVertex *v1, GVertex *v2)
 }
 
 gmshEdge::gmshEdge(GModel *model, int num)
-  : GEdge(model, num, 0, 0), c(0)
+  : GEdge(model, num, 0, 0)
 {
+  c = Create_Curve(num, MSH_SEGM_DISCRETE, 0, NULL, NULL, -1, -1, 0., 1.);
+  Tree_Add(THEM->Curves, &c);
+  CreateReversedCurve(THEM, c);
 }
 
 gmshEdge::~gmshEdge()
@@ -23,14 +29,11 @@ gmshEdge::~gmshEdge()
 
 Range<double> gmshEdge::parBounds(int i) const
 { 
-  if(!c) return(Range<double>(0., 1.));
   return(Range<double>(c->ubeg, c->uend));
 }
 
 SBoundingBox3d gmshEdge::bounds() const
 {
-  if(!c) return SBoundingBox3d(SPoint3(0., 0., 0.));
-
   double xmin = 0., ymin = 0., zmin = 0.;
   double xmax = 0., ymax = 0., zmax = 0.;
   for (int i = 0; i < 20; i++){
@@ -59,14 +62,12 @@ SBoundingBox3d gmshEdge::bounds() const
 
 GPoint gmshEdge::point(double par) const
 {
-  if(!c) return GPoint(0., 0., 0., this, 0.);
   Vertex a = InterpolateCurve(c, par, 0);
   return GPoint(a.Pos.X,a.Pos.Y,a.Pos.Z,this,par);
 }
 
 GPoint gmshEdge::closestPoint(const SPoint3 & qp)
 {
-  if(!c) return GPoint(0., 0., 0., this, 0.);
   Vertex v;
   Vertex a;
   Vertex der;
@@ -92,7 +93,6 @@ SVector3 gmshEdge::firstDer(double par) const
 
 double gmshEdge::parFromPoint(const SPoint3 &pt) const
 {
-  if(!c) return 0;
   Vertex v;
   Vertex a;
   Vertex der;
@@ -121,8 +121,6 @@ bool gmshEdge::periodic(int dim) const
 
 GEntity::GeomType gmshEdge::geomType() const
 {
-  if(!c) return DiscreteCurve;
-
   switch (c->Typ){
   case MSH_SEGM_LINE : return Line;
   case MSH_SEGM_PARAMETRIC : return ParametricCurve;
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 407fe058fcea0fd55e02cb7e1c5f0388357bf1b2..ef7f9b911d1603183cc081af4a35c2c66d2003f1 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -4,8 +4,12 @@
 #include "Interpolation.h"
 #include "CAD.h"
 #include "Geo.h"
+#include "Mesh.h"
+#include "Create.h"
 #include "Utils.h"
 
+extern Mesh *THEM;
+
 gmshFace::gmshFace(GModel *m, Surface *face)
   : GFace(m, face->Num), s(face)
 {
@@ -22,8 +26,10 @@ gmshFace::gmshFace(GModel *m, Surface *face)
 }
 
 gmshFace::gmshFace(GModel *m, int num)
-  : GFace(m, num), s(0)
+  : GFace(m, num)
 {
+  s = Create_Surface(num, MSH_SURF_DISCRETE);
+  Tree_Add(THEM->Surfaces, &s);
 }
 
 Range<double> gmshFace::parBounds(int i) const
@@ -51,8 +57,6 @@ SBoundingBox3d gmshFace::bounds() const
 
 SVector3 gmshFace::normal(const SPoint2 &param) const
 {
-  if(!s) return SVector3(0., 0., 0.);
-
   Vertex vu = InterpolateSurface(s, param[0], param[1], 1, 1);
   Vertex vv = InterpolateSurface(s, param[0], param[1], 1, 2);
   Vertex n = vu % vv;
@@ -62,9 +66,6 @@ SVector3 gmshFace::normal(const SPoint2 &param) const
 
 Pair<SVector3,SVector3> gmshFace::firstDer(const SPoint2 &param) const
 {
-  if(!s) return Pair<SVector3,SVector3>(SVector3(0., 0., 0.),
-					SVector3(0., 0., 0.));
-
   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),
@@ -94,8 +95,6 @@ void computePlaneDatas (const GFace *gf, double VX[3],double VY[3],double &x, do
 
 GPoint gmshFace::point(double par1,double par2) const
 {
-  if(!s) return GPoint(0., 0., 0., this);
-
   double pp[2]={par1,par2};
   if(s->Typ == MSH_SURF_PLAN){
     double x,y,z,VX[3],VY[3];
@@ -112,8 +111,6 @@ GPoint gmshFace::point(double par1,double par2) const
 
 GPoint gmshFace::closestPoint(const SPoint3 & qp)
 {
-  if(!s) return GPoint(0., 0., 0., this);
-
   Vertex v;
   v.Pos.X = qp.x();
   v.Pos.Y = qp.y();
@@ -137,8 +134,6 @@ int gmshFace::containsParam(const SPoint2 &pt) const
 
 SPoint2 gmshFace::parFromPoint(const SPoint3 &qp) const
 {
-  if(!s) return SPoint2(0., 0.);
-
   double u,v;
   if (s->Typ == MSH_SURF_PLAN){
     double x,y,z,VX[3],VY[3];
@@ -171,8 +166,6 @@ bool gmshFace::degenerate(int dim) const
 
 GEntity::GeomType gmshFace::geomType() const
 {
-  if(!s) return DiscreteSurface;
-
   switch(s->Typ){
   case MSH_SURF_NURBS: return Nurb;
   case MSH_SURF_PLAN: return Plane;
@@ -187,8 +180,9 @@ int gmshFace::geomDirection() const
 }
 
 void * gmshFace::getNativePtr() const
-{ return s; }
-
+{ 
+  return s; 
+}
 
 int gmshFace::containsPoint(const SPoint3 &pt) const
 { 
diff --git a/Geo/gmshModel.cpp b/Geo/gmshModel.cpp
index 2e1537f3c93a863edc5050a40737986b93ce31bf..7c58c5da088c38a5465cb16d3e264d88c7188cfa 100644
--- a/Geo/gmshModel.cpp
+++ b/Geo/gmshModel.cpp
@@ -37,7 +37,7 @@ void gmshModel::convertFromUglyOldDataStructuresgmshModel()
     for(int i = 0; i < List_Nbr(curves); i++){
       Curve *c;
       List_Read(curves, i, &c);
-      if(c->Num >= 0){
+      if(c->Num >= 0 && c->beg && c->end){
 	if(points.find(c->beg) == points.end()){
 	  points.insert(c->beg);
 	  gmshVertex *v = new gmshVertex(this, c->beg);
diff --git a/Geo/gmshRegion.cpp b/Geo/gmshRegion.cpp
index 361ca567653d6ebf4e2d00826fbd31bc5bc8d3be..489e423cc2a9bdd663141da1b81faaafbc6ef8c9 100644
--- a/Geo/gmshRegion.cpp
+++ b/Geo/gmshRegion.cpp
@@ -5,6 +5,11 @@
 #include "Interpolation.h"
 #include "CAD.h"
 #include "Geo.h"
+#include "Mesh.h"
+#include "Create.h"
+#include "Utils.h"
+
+extern Mesh *THEM;
 
 gmshRegion::gmshRegion(GModel *m, ::Volume * volume)
   : GRegion(m, volume->Num), v(volume)
@@ -22,6 +27,8 @@ gmshRegion::gmshRegion(GModel *m, ::Volume * volume)
 }
 
 gmshRegion::gmshRegion(GModel *m, int num)
-  : GRegion(m, num), v(0)
+  : GRegion(m, num)
 {
+  v = Create_Volume(num, MSH_VOLUME_DISCRETE);
+  Tree_Add(THEM->Volumes, &v);
 }