diff --git a/Geo/GModelIO_GEO.cpp b/Geo/GModelIO_GEO.cpp
index baabc9ff4c0a8209a388ab2b6a469139a31d8f53..c86089b72707297f450e339b8dc958479750fae1 100644
--- a/Geo/GModelIO_GEO.cpp
+++ b/Geo/GModelIO_GEO.cpp
@@ -483,6 +483,10 @@ static bool skipVertex(GVertex *gv)
 int GModel::writeGEO(const std::string &name, bool printLabels, bool onlyPhysicals)
 {
   FILE *fp = Fopen(name.c_str(), "w");
+  if(!fp){
+    Msg::Error("Could not open file '%s'", name.c_str());
+    return 0;
+  }
 
   std::map<double, std::string> meshSizeParameters;
   int cpt = 0;
@@ -533,6 +537,6 @@ int GModel::writeGEO(const std::string &name, bool printLabels, bool onlyPhysica
   if(getFields()->getBackgroundField() > 0)
     fprintf(fp, "Background Field = %i;\n", getFields()->getBackgroundField());
 
-  if(fp) fclose(fp);
+  fclose(fp);
   return 1;
 }
diff --git a/Mesh/BackgroundMesh2D.cpp b/Mesh/BackgroundMesh2D.cpp
index c023bc0362eba2b0cd628c9a28e9774cc2f17d66..f3a062ca774bf4c0de63525bf5c622296a19ede5 100644
--- a/Mesh/BackgroundMesh2D.cpp
+++ b/Mesh/BackgroundMesh2D.cpp
@@ -60,13 +60,21 @@ void normalizeAngle(double &angle)
 
 void backgroundMesh2D::create_face_mesh()
 {
-  quadsToTriangles(dynamic_cast<GFace*>(gf), 100000);
+  GFace *face = dynamic_cast<GFace*>(gf);
+  if(!face){
+    Msg::Error("Entity is not a face in background mesh");
+    return;
+  }
+
+  quadsToTriangles(face, 100000);
 
   // storing the initial mesh from GFace
   tempTR.clear();
-  GFace *face = dynamic_cast<GFace*>(gf);
+
   for(unsigned int i = 0; i < face->triangles.size(); i++)
-    tempTR.push_back(new MTriangle(face->triangles[i]->getVertex(0), face->triangles[i]->getVertex(1), face->triangles[i]->getVertex(2)));
+    tempTR.push_back(new MTriangle(face->triangles[i]->getVertex(0),
+                                   face->triangles[i]->getVertex(1),
+                                   face->triangles[i]->getVertex(2)));
 
   // avoid computing curvatures on the fly : only on the
   // BGM computes once curvatures at each node
@@ -121,19 +129,23 @@ void backgroundMesh2D::reset(bool erase_2D3D)
   }
 }
 
-
-void backgroundMesh2D::unset(){
+void backgroundMesh2D::unset()
+{
   for (unsigned int i = 0; i < vertices.size(); i++) delete vertices[i];
   for (unsigned int i = 0; i < getNumMeshElements(); i++) delete elements[i];
   if (octree)delete octree;
   octree=NULL;
 }
 
-
-void backgroundMesh2D::create_mesh_copy(){
+void backgroundMesh2D::create_mesh_copy()
+{
   // TODO: useful to extend it to other elements ???
   //std::set<SPoint2> myBCNodes;
   GFace *face = dynamic_cast<GFace*>(gf);
+  if(!face){
+    Msg::Error("Entity is not a face in background mesh");
+    return;
+  }
   for (unsigned int i = 0; i < face->triangles.size(); i++){
     MTriangle *e = face->triangles[i];
     MVertex *news[3];
@@ -158,8 +170,14 @@ void backgroundMesh2D::create_mesh_copy(){
 }
 
 
-GPoint backgroundMesh2D::get_GPoint_from_MVertex(const MVertex *v)const{
-  return dynamic_cast<GFace*>(gf)->point(SPoint2(v->x(),v->y()));
+GPoint backgroundMesh2D::get_GPoint_from_MVertex(const MVertex *v)const
+{
+  GFace *face = dynamic_cast<GFace*>(gf);
+  if(!face){
+    Msg::Error("Entity is not a face in background mesh");
+    return GPoint();
+  }
+  return face->point(SPoint2(v->x(),v->y()));
 }
 
 
@@ -171,7 +189,10 @@ backgroundMesh2D::backgroundMesh2D(GFace *_gf, bool erase_2D3D):BGMBase(2,_gf),s
     // now, the new mesh has been copied in local in backgroundMesh2D, deleting the mesh
     // from GFace, back to the previous one !
     GFace *face = dynamic_cast<GFace*>(gf);
-    face->triangles = tempTR;
+    if(!face)
+      Msg::Error("Entity is not a face in background mesh");
+    else
+      face->triangles = tempTR;
   }
 }
 
@@ -181,7 +202,10 @@ backgroundMesh2D::~backgroundMesh2D()
 }
 
 
-void backgroundMesh2D::propagateValues(DoubleStorageType &dirichlet, simpleFunction<double> &eval_diffusivity, bool in_parametric_plane){
+void backgroundMesh2D::propagateValues(DoubleStorageType &dirichlet,
+                                       simpleFunction<double> &eval_diffusivity,
+                                       bool in_parametric_plane)
+{
 #if defined(HAVE_SOLVER)
   linearSystem<double> *_lsys = 0;
 #if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
@@ -207,12 +231,16 @@ void backgroundMesh2D::propagateValues(DoubleStorageType &dirichlet, simpleFunct
   // Number vertices
   std::set<MVertex*> vs;
   GFace *face = dynamic_cast<GFace*>(gf);
+  if(!face){
+    Msg::Error("Entity is not a face in background mesh");
+    delete _lsys;
+    return;
+  }
   for (unsigned int k = 0; k < face->triangles.size(); k++)
     for (int j=0;j<3;j++)vs.insert(face->triangles[k]->getVertex(j));
   for (unsigned int k = 0; k < face->quadrangles.size(); k++)
     for (int j=0;j<4;j++)vs.insert(face->quadrangles[k]->getVertex(j));
 
-
   std::map<MVertex*,SPoint3> theMap;
   if ( in_parametric_plane) {
     for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
@@ -254,10 +282,15 @@ void backgroundMesh2D::propagateValues(DoubleStorageType &dirichlet, simpleFunct
 #endif
 }
 
-
-void backgroundMesh2D::computeSizeField(){
+void backgroundMesh2D::computeSizeField()
+{
   GFace *face = dynamic_cast<GFace*>(gf);
-  list<GEdge*> e;// = face->edges();
+  if(!face){
+    Msg::Error("Entity is not a face in background mesh");
+    return;
+  }
+
+  list<GEdge*> e;
   replaceMeshCompound(face, e);
   list<GEdge*>::const_iterator it = e.begin();
   DoubleStorageType sizes;
@@ -303,7 +336,8 @@ inline double myAngle (const SVector3 &a, const SVector3 &b, const SVector3 &d){
 }
 
 
-void backgroundMesh2D::updateSizes(){
+void backgroundMesh2D::updateSizes()
+{
   DoubleStorageType::iterator itv = sizeField.begin();
   for ( ; itv != sizeField.end(); ++itv){
     SPoint2 p;
@@ -319,6 +353,10 @@ void backgroundMesh2D::updateSizes(){
     }
     else{
       GFace *face = dynamic_cast<GFace*>(gf);
+      if(!face){
+        Msg::Error("Entity is not a face in background mesh");
+        return;
+      }
       reparamMeshVertexOnFace(v, face, p);
       lc = sizeFactor * BGM_MeshSize(face, p.x(), p.y(), v->x(), v->y(), v->z());
     }
@@ -363,7 +401,10 @@ frameFieldBackgroundMesh2D::frameFieldBackgroundMesh2D(GFace *_gf):backgroundMes
   // now, the new mesh has been copied in local in backgroundMesh2D, deleting the mesh
   // from GFace, back to the previous one !
   GFace *face = dynamic_cast<GFace*>(gf);
-  face->triangles = tempTR;
+  if(!face)
+    Msg::Error("Entity is not a face in background mesh");
+  else
+    face->triangles = tempTR;
 }
 
 frameFieldBackgroundMesh2D::~frameFieldBackgroundMesh2D(){}
@@ -437,6 +478,11 @@ void frameFieldBackgroundMesh2D::computeCrossField(simpleFunction<double> &eval_
 
   list<GEdge*> e;
   GFace *face = dynamic_cast<GFace*>(gf);
+  if(!face){
+    Msg::Error("Entity is not a face in background mesh");
+    return;
+  }
+
   replaceMeshCompound(face, e);
 
   list<GEdge*>::const_iterator it = e.begin();
@@ -523,6 +569,10 @@ void frameFieldBackgroundMesh2D::eval_crossfield(MVertex *vert, STensor3 &cf)
 {
   SPoint2 parampoint;
   GFace *face = dynamic_cast<GFace*>(gf);
+  if(!face){
+    Msg::Error("Entity is not a face in background mesh");
+    return;
+  }
   reparamMeshVertexOnFace(vert, face, parampoint);
   return eval_crossfield(parampoint[0], parampoint[1], cf);
 }
@@ -596,10 +646,16 @@ void frameFieldBackgroundMesh2D::exportCrossField(const std::string &filename)
 }
 
 // returns the cross field as a pair of othogonal vectors (NOT in parametric coordinates, but real 3D coordinates)
-Pair<SVector3, SVector3> frameFieldBackgroundMesh2D::compute_crossfield_directions(double u,double v, double angle_current)
+Pair<SVector3, SVector3> frameFieldBackgroundMesh2D::compute_crossfield_directions(double u, double v,
+                                                                                   double angle_current)
 {
   // get the unit normal at that point
   GFace *face = dynamic_cast<GFace*>(gf);
+  if(!face){
+    Msg::Error("Entity is not a face in background mesh");
+    return Pair<SVector3,SVector3>(SVector3(), SVector3());
+  }
+
   Pair<SVector3, SVector3> der = face->firstDer(SPoint2(u,v));
   SVector3 s1 = der.first();
   SVector3 s2 = der.second();
@@ -623,7 +679,6 @@ Pair<SVector3, SVector3> frameFieldBackgroundMesh2D::compute_crossfield_directio
 
 bool frameFieldBackgroundMesh2D::compute_RK_infos(double u,double v, double x, double y, double z, RK_form &infos)
 {
-
   // check if point is in domain
   if (!inDomain(u,v)) return false;
 
@@ -635,6 +690,11 @@ bool frameFieldBackgroundMesh2D::compute_RK_infos(double u,double v, double x, d
 
   // get the unit normal at that point
   GFace *face = dynamic_cast<GFace*>(gf);
+  if(!face){
+    Msg::Error("Entity is not a face in background mesh");
+    return false;
+  }
+
   Pair<SVector3, SVector3> der = face->firstDer(SPoint2(u,v));
   SVector3 s1 = der.first();
   SVector3 s2 = der.second();
diff --git a/Mesh/BackgroundMesh2D.h b/Mesh/BackgroundMesh2D.h
index 3b1712dbde45ea46cc593b473403c5074594195f..1bfba1b1847dcd211cb785c7734cbea69c7fd82a 100644
--- a/Mesh/BackgroundMesh2D.h
+++ b/Mesh/BackgroundMesh2D.h
@@ -24,117 +24,99 @@ using namespace std;
    Those triangles are local to the backgroundMesh2D so that
    they do not depend on the actual mesh that can be deleted.
    It extends the sizefield from the edges.
- */
+*/
 class backgroundMesh2D : public BGMBase {
-  protected:
-    virtual void propagateValues(DoubleStorageType &dirichlet, simpleFunction<double> &eval_diffusivity, bool in_parametric_plane = false);
-    virtual void computeSizeField();
-
-    virtual inline unsigned int getNumMeshElements()const{return elements.size();}
-    virtual const MElement* getElement(unsigned int i)const;
-
-    virtual GPoint get_GPoint_from_MVertex(const MVertex *) const;
-
-    // only 2D:
-    void updateSizes();
-    // copy the mesh stored in GFace in local
-    void create_mesh_copy();
-    // creates a mesh of GFace and store it in local !!!, does not store the mesh in GFace !
-    void create_face_mesh();
-
-    double sizeFactor;
-
-
-    std::vector<MTriangle*> tempTR;
-    vector<MElement*> elements;
-    vector<MVertex*> vertices;
-
-    map<MVertex*,MVertex*> _3Dto2D;
-    map<MVertex*,MVertex*> _2Dto3D;
-
-
-  public:
-    backgroundMesh2D(GFace *, bool erase_2D3D=true);
-    virtual ~backgroundMesh2D();
-
-    virtual MElementOctree* getOctree();
-
-    // TODO: only 2D
-    virtual void reset(bool erase_2D3D=true);// deletes everything and rebuild with GFace*
-    virtual void unset();// deletes everything... called by destructor.
-
-    // not used !!!! TODO !!!
-    void setSizeFactor (double s) {sizeFactor = s;}
-
-
-    virtual vector<MVertex*>::iterator beginvertices(){return vertices.begin();}
-    virtual vector<MVertex*>::iterator endvertices(){return vertices.end();}
-    virtual vector<MElement*>::iterator beginelements(){return elements.begin();}
-    virtual vector<MElement*>::iterator endelements(){return elements.end();}
-    virtual vector<MVertex*>::const_iterator beginvertices()const{return vertices.begin();}
-    virtual vector<MVertex*>::const_iterator endvertices()const{return vertices.end();}
-    virtual vector<MElement*>::const_iterator beginelements()const{return elements.begin();}
-    virtual vector<MElement*>::const_iterator endelements()const{return elements.end();}
-
+protected:
+  virtual void propagateValues(DoubleStorageType &dirichlet,
+                               simpleFunction<double> &eval_diffusivity,
+                               bool in_parametric_plane = false);
+  virtual void computeSizeField();
+  virtual inline unsigned int getNumMeshElements()const{return elements.size();}
+  virtual const MElement* getElement(unsigned int i)const;
+  virtual GPoint get_GPoint_from_MVertex(const MVertex *) const;
+
+  // only 2D:
+  void updateSizes();
+  // copy the mesh stored in GFace in local
+  void create_mesh_copy();
+  // creates a mesh of GFace and store it in local !!!, does not store the mesh in GFace !
+  void create_face_mesh();
+
+  double sizeFactor;
+  std::vector<MTriangle*> tempTR;
+  vector<MElement*> elements;
+  vector<MVertex*> vertices;
+  map<MVertex*,MVertex*> _3Dto2D;
+  map<MVertex*,MVertex*> _2Dto3D;
+
+public:
+  backgroundMesh2D(GFace *, bool erase_2D3D=true);
+  virtual ~backgroundMesh2D();
+  virtual MElementOctree* getOctree();
+
+  // TODO: only 2D
+  virtual void reset(bool erase_2D3D=true);// deletes everything and rebuild with GFace*
+  virtual void unset();// deletes everything... called by destructor.
+
+  // not used !!!! TODO !!!
+  void setSizeFactor (double s) {sizeFactor = s;}
+
+  virtual vector<MVertex*>::iterator beginvertices(){ return vertices.begin(); }
+  virtual vector<MVertex*>::iterator endvertices(){ return vertices.end(); }
+  virtual vector<MElement*>::iterator beginelements(){ return elements.begin(); }
+  virtual vector<MElement*>::iterator endelements(){ return elements.end(); }
+  virtual vector<MVertex*>::const_iterator beginvertices() const { return vertices.begin(); }
+  virtual vector<MVertex*>::const_iterator endvertices() const { return vertices.end(); }
+  virtual vector<MElement*>::const_iterator beginelements() const { return elements.begin(); }
+  virtual vector<MElement*>::const_iterator endelements() const { return elements.end(); }
 };
 
-
-class RK_form{// informations for RK at one point
-  public:
-    RK_form(){
-    }
-//    RK_form(const RK_form &other){
-//      t1 = other.t1;
-//      t2 = other.t2;
-//      h.first = other.h.first;
-//      h.second = other.h.second;
-//      paramh.first = other.paramh.first;
-//      paramh.second = other.paramh.second;
-//      paramt1 = other.paramt1;
-//      paramt2 = other.paramt2;
-//      angle = other.angle;
-//    }
-
-    SMetric3 metricField;
-    SVector3 t1, t2;// 3D cross field directions
-    SVector3 normal;// 3D cross field directions
-    pair<double,double> h;// 3D sizes
-    pair<double,double> paramh;// sizes in parametric domain
-    SPoint2 paramt1,paramt2;
-    double angle,localsize;
+class RK_form{ // informations for RK at one point
+public:
+  RK_form() : angle(0.), localsize(0.)
+  {
+  }
+  SMetric3 metricField;
+  SVector3 t1, t2;// 3D cross field directions
+  SVector3 normal;// 3D cross field directions
+  pair<double,double> h;// 3D sizes
+  pair<double,double> paramh;// sizes in parametric domain
+  SPoint2 paramt1, paramt2;
+  double angle, localsize;
 };
 
+class frameFieldBackgroundMesh2D : public backgroundMesh2D {
+private:
+  // specification for cross field
+  DoubleStorageType angles;  // an attached angle
+  DoubleStorageType smoothness;
+  void computeCrossField(simpleFunction<double> &eval_diffusivity);
+  void computeSmoothness();
 
+public:
+  frameFieldBackgroundMesh2D(GFace *_gf);
+  virtual ~frameFieldBackgroundMesh2D();
 
-class frameFieldBackgroundMesh2D : public backgroundMesh2D{
-  private:
-    // specification for cross field
-    DoubleStorageType angles;  // an attached angle
-    DoubleStorageType smoothness;
-
-    void computeCrossField(simpleFunction<double> &eval_diffusivity);
-    void computeSmoothness();
-
-  public:
-    frameFieldBackgroundMesh2D(GFace *_gf);
-    virtual ~frameFieldBackgroundMesh2D();
-
-    virtual void reset(bool erase_2D3D=true);
+  virtual void reset(bool erase_2D3D=true);
 
-    double angle(double u, double v);
-    double angle(MVertex *v);
+  double angle(double u, double v);
+  double angle(MVertex *v);
 
-    double get_smoothness(double u, double v);
-    double get_smoothness(MVertex *v);
+  double get_smoothness(double u, double v);
+  double get_smoothness(MVertex *v);
 
-    void eval_crossfield(double u, double v, STensor3 &cf);
-    void eval_crossfield(MVertex *vert, STensor3 &cf);
+  void eval_crossfield(double u, double v, STensor3 &cf);
+  void eval_crossfield(MVertex *vert, STensor3 &cf);
 
-    void exportCrossField(const string &filename);
-    Pair<SVector3, SVector3> compute_crossfield_directions(double u,double v, double angle_current);
-    bool compute_RK_infos(double u,double v, double x, double y, double z,RK_form &infos);
+  void exportCrossField(const string &filename);
+  Pair<SVector3, SVector3> compute_crossfield_directions(double u,double v,
+                                                         double angle_current);
+  bool compute_RK_infos(double u,double v, double x, double y, double z,RK_form &infos);
 
-    inline void exportSmoothness(const string &filename) const{export_scalar(filename,smoothness);};
+  inline void exportSmoothness(const string &filename) const
+  {
+    export_scalar(filename,smoothness);
+  }
 };
 
 #endif
diff --git a/Mesh/BackgroundMesh3D.cpp b/Mesh/BackgroundMesh3D.cpp
index 9ff8296e71ea566f58eeb5a20ad545b12b9e79fa..1541328a45c302af06ecc9697c77eab19e0c3a2c 100644
--- a/Mesh/BackgroundMesh3D.cpp
+++ b/Mesh/BackgroundMesh3D.cpp
@@ -61,6 +61,10 @@ void backgroundMesh3D::computeSizeField()
 
   // fills dirichlet BC
   GRegion *region = dynamic_cast<GRegion*>(gf);
+  if(!region){
+    Msg::Error("Entity is not a region in background mesh");
+    return;
+  }
   list<GFace*> faces = region->faces();
   MVertex *v;
   MElement *e;
@@ -68,6 +72,10 @@ void backgroundMesh3D::computeSizeField()
   for (list<GFace*>::iterator it=faces.begin();it!=faces.end();it++){// for all GFace
     GFace *face = *it;
     frameFieldBackgroundMesh2D *bgm2d = dynamic_cast<frameFieldBackgroundMesh2D*>(BGMManager::get(face));
+    if(!bgm2d){
+      Msg::Error("Background mesh is not a 2D frame field");
+      return;
+    }
     for (unsigned int i=0;i<face->getNumMeshElements();i++){// for all elements
       e = face->getMeshElement(i);
       if (e->getDim()!=2) continue;// of dim=2
@@ -94,6 +102,10 @@ void backgroundMesh3D::propagateValues(DoubleStorageType &dirichlet,
 {
   // same as Size_field::solve()
   GRegion *gr = dynamic_cast<GRegion*>(gf);
+  if(!gr){
+    Msg::Error("Entity is not a region in background mesh");
+    return;
+  }
 
 #if defined(HAVE_PETSC)
   linearSystem<double>* system = 0;
@@ -195,6 +207,10 @@ MElementOctree* backgroundMesh3D::getOctree()
 {
   if(!octree){
     GRegion *gr = dynamic_cast<GRegion*>(gf);
+    if(!gr){
+      Msg::Error("Entity is not a region in background mesh");
+      return 0;
+    }
     Msg::Debug("Rebuilding BackgroundMesh element octree");
     std::vector<MElement*> copy;// !!!
     for (std::vector<MTetrahedron*>::iterator it = gr->tetrahedra.begin();
@@ -578,6 +594,10 @@ void frameFieldBackgroundMesh3D::initiate_crossfield()
 
   // first, for boundaries :
   GRegion *gr = dynamic_cast<GRegion*>(gf);
+  if(!gr){
+    Msg::Error("Entity is not a region in background mesh");
+    return;
+  }
   list<GFace*> faces = gr->faces();
   // here, not using the gm2D since we are interested by the new 2D vertices,
   // not the old (now erased) ones... alternative would be to reset the 2DBGM...
@@ -585,6 +605,10 @@ void frameFieldBackgroundMesh3D::initiate_crossfield()
     GFace *face = *it;
     frameFieldBackgroundMesh2D *bgm2d =
       dynamic_cast<frameFieldBackgroundMesh2D*>(BGMManager::get(face));
+    if(!bgm2d){
+      Msg::Error("Background mesh is not a 2D frame field");
+      return;
+    }
     // initializing the vertices on the faces
     for (unsigned int i=0;i<face->getNumMeshElements();i++){// for all elements
       MElement *e = face->getMeshElement(i);
@@ -721,6 +745,10 @@ double frameFieldBackgroundMesh3D::get_vectorial_smoothness(const int idir,
 void frameFieldBackgroundMesh3D::build_vertex_to_element_table()
 {
   GRegion *gr = dynamic_cast<GRegion*>(gf);
+  if(!gr){
+    Msg::Error("Entity is not a region in background mesh");
+    return;
+  }
   MElement *e;
   MVertex *v;
 
@@ -746,12 +774,20 @@ void frameFieldBackgroundMesh3D::build_vertex_to_element_table()
 const MElement* backgroundMesh3D::getElement(unsigned int i)const
 {
   GRegion *gr = dynamic_cast<GRegion*>(gf);
+  if(!gr){
+    Msg::Error("Entity is not a region in background mesh");
+    return 0;
+  }
   return gr->getMeshElement(i);
 }
 
 unsigned int backgroundMesh3D::getNumMeshElements()const
 {
   GRegion *gr = dynamic_cast<GRegion*>(gf);
+  if(!gr){
+    Msg::Error("Entity is not a region in background mesh");
+    return 0;
+  }
   return gr->getNumMeshElements();
 }
 
diff --git a/Mesh/pointInsertionRTreeTools.h b/Mesh/pointInsertionRTreeTools.h
index 42945a89ee9496482d89547b55d46ab0ac4656a5..b9ffba8105c1d8a89635be5f0f29a1752733f60c 100644
--- a/Mesh/pointInsertionRTreeTools.h
+++ b/Mesh/pointInsertionRTreeTools.h
@@ -106,23 +106,29 @@ private:
   SVector3 *vec;
   bool ok;
 public:
-  Wrapper3D():ok(true){};
-  Wrapper3D(MVertex* _i,MVertex* _p):individual(_i), parent(_p),ok(true){};
-  ~Wrapper3D(){};
-  void set_ok(bool b){ok=b;};
-  void set_individual(MVertex *vertex){individual=vertex;};
-  void set_parent(MVertex *vertex){parent=vertex;};
-  void set_size(double *h){size=h;};
-  void set_crossfield(STensor3 *_cf){cf=_cf;};
-  void set_direction(SVector3 *_v){vec=_v;};
-  bool get_ok(){return ok;};
-  void set_bgm(frameFieldBackgroundMesh3D *bgm){bgmesh = bgm;};
-  frameFieldBackgroundMesh3D * bgm(){return bgmesh;};
-  MVertex* get_individual(){return individual;};
-  MVertex* get_parent(){return parent;};
-  STensor3* get_crossfield(){return cf;};
-  SVector3* get_direction(){return vec;};
-  double* get_size(){return size;};
+  Wrapper3D()
+    : individual(0), parent(0), size(0), cf(0), vec(0), ok(true)
+  {
+  }
+  Wrapper3D(MVertex* _i,MVertex* _p)
+    : individual(_i), parent(_p), size(0), cf(0), vec(0), ok(true)
+  {
+  }
+  ~Wrapper3D(){}
+  void set_ok(bool b){ ok = b; }
+  void set_individual(MVertex *vertex){ individual = vertex; }
+  void set_parent(MVertex *vertex){ parent = vertex; }
+  void set_size(double *h){ size = h; }
+  void set_crossfield(STensor3 *_cf){ cf = _cf; }
+  void set_direction(SVector3 *_v){ vec = _v; }
+  bool get_ok(){ return ok; }
+  void set_bgm(frameFieldBackgroundMesh3D *bgm){ bgmesh = bgm; }
+  frameFieldBackgroundMesh3D * bgm(){ return bgmesh; }
+  MVertex* get_individual(){ return individual; }
+  MVertex* get_parent(){ return parent; }
+  STensor3* get_crossfield(){ return cf; }
+  SVector3* get_direction(){ return vec; }
+  double* get_size(){ return size; }
 };
 
 extern double infinity_distance_3D(const MVertex *v1,const MVertex *v2,STensor3 &cf);
@@ -134,8 +140,8 @@ extern void fill_min_max(double x,double y,double z,double h,double *min,double
 // listOfPoints AND in RTree: larger memory footprint but less CPU time...
 class smoothness_vertex_pair{
 public:
-  smoothness_vertex_pair(){};
-  ~smoothness_vertex_pair(){};
+  smoothness_vertex_pair() : rank(0.), size(0.), v(0), dir(0), layer(0) {}
+  ~smoothness_vertex_pair(){}
   STensor3 cf;
   SVector3 direction;
   double rank, size;
@@ -180,30 +186,30 @@ public:
 
 class listOfPointsScalarSmoothness : public listOfPoints{
 public:
-  listOfPointsScalarSmoothness(){ };
+  listOfPointsScalarSmoothness(){ }
   virtual ~listOfPointsScalarSmoothness()
   {
     while (!empty())
       erase_first();
-  };
-  virtual void insert(smoothness_vertex_pair *svp){ points.insert(svp); };
-  virtual unsigned int size(){ return points.size(); };
-  virtual MVertex* get_first_vertex(){ return (*points.begin())->v; };
-  virtual STensor3 get_first_crossfield(){ return (*points.begin())->cf; };
-  virtual double get_first_size(){ return (*points.begin())->size; };
-  virtual int get_first_layer(){ return (*points.begin())->layer; };
+  }
+  virtual void insert(smoothness_vertex_pair *svp){ points.insert(svp); }
+  virtual unsigned int size(){ return points.size(); }
+  virtual MVertex* get_first_vertex(){ return (*points.begin())->v; }
+  virtual STensor3 get_first_crossfield(){ return (*points.begin())->cf; }
+  virtual double get_first_size(){ return (*points.begin())->size; }
+  virtual int get_first_layer(){ return (*points.begin())->layer; }
   virtual SVector3 get_first_direction()
   {
     Msg::Error("listOfPointsScalarSmoothness::get_first_direction NOT applicable");
     return SVector3(0.);
-  };
+  }
   virtual void erase_first()
   {
     smoothness_vertex_pair *ptr = *(points.begin());
     points.erase(points.begin());
     delete ptr;
-  };
-  virtual bool empty(){ return points.empty(); };
+  }
+  virtual bool empty(){ return points.empty(); }
 
 protected:
   std::set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
@@ -211,12 +217,12 @@ protected:
 
 class listOfPointsVectorialSmoothness : public listOfPointsScalarSmoothness{
 public:
-  listOfPointsVectorialSmoothness(){};
+  listOfPointsVectorialSmoothness(){}
   virtual ~listOfPointsVectorialSmoothness(){
     while (!empty())
       erase_first();
-  };
-  virtual SVector3 get_first_direction(){ return (*points.begin())->direction; };
+  }
+  virtual SVector3 get_first_direction(){ return (*points.begin())->direction; }
 protected:
   std::set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
 };
@@ -228,24 +234,24 @@ public:
     while (!empty())
       erase_first();
   };
-  virtual void insert(smoothness_vertex_pair *svp){ points.push(svp); };
-  virtual unsigned int size(){ return points.size(); };
-  virtual MVertex* get_first_vertex(){ return (points.front())->v; };
-  virtual STensor3 get_first_crossfield(){ return (points.front())->cf; };
-  virtual double get_first_size(){ return (points.front())->size; };
-  virtual int get_first_layer(){ return (points.front())->layer; };
+  virtual void insert(smoothness_vertex_pair *svp){ points.push(svp); }
+  virtual unsigned int size(){ return points.size(); }
+  virtual MVertex* get_first_vertex(){ return (points.front())->v; }
+  virtual STensor3 get_first_crossfield(){ return (points.front())->cf; }
+  virtual double get_first_size(){ return (points.front())->size; }
+  virtual int get_first_layer(){ return (points.front())->layer; }
   virtual SVector3 get_first_direction()
   {
     Msg::Error("listOfPointsFifo::get_first_direction NOT applicable");
     return SVector3(0.);
-  };
+  }
   virtual void erase_first()
   {
     smoothness_vertex_pair *ptr = points.front();
     points.pop();
     delete ptr;
-  };
-  virtual bool empty(){ return points.empty(); };
+  }
+  virtual bool empty(){ return points.empty(); }
 
 protected:
   std::queue<smoothness_vertex_pair*> points;
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index 0be2818e142c3d38f631aa1dc5b3381cbbc1c85e..942e91b5de35f3981f98461f77e45a9b916b0a0a 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -110,6 +110,10 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
 
   Msg::Info("Writing %s", _fileName.c_str());
   FILE *fName = Fopen(_fileName.c_str(),"w");
+  if(!fName){
+    Msg::Error("Could not open file '%s'", _fileName.c_str());
+    return;
+  }
   fprintf(fName, "View \"distance \"{\n");
 
   for (unsigned int ii=0; ii<_entities.size(); ii++) {
@@ -201,7 +205,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
   std::vector<GEntity*> _entities;
   GModel::current()->getEntities(_entities);
   if (!_entities.size() || !_entities[_entities.size()-1]->getMeshElement(0)) {
-    Msg::Error("This plugin needs a mesh !");
+    Msg::Error("This plugin needs a mesh!");
     return view;
   }
 
@@ -495,7 +499,16 @@ PView *GMSH_DistancePlugin::execute(PView *v)
     data2->setName("ortogonal field");
 
     Msg::Info("Writing  orthogonal.pos");
-    FILE * f5 = Fopen("orthogonal.pos","w");
+    FILE *f5 = Fopen("orthogonal.pos","w");
+    if(!f5){
+      Msg::Error("Could not open file 'orthogonal.pos'");
+#if defined(HAVE_SOLVER)
+      delete lsys;
+      delete dofView;
+#endif
+      return view;
+    }
+
     fprintf(f5,"View \"orthogonal\"{\n");
     for (std::vector<MElement* >::iterator it = allElems.begin();
          it != allElems.end(); it++){
diff --git a/utils/misc/package_gmsh_getdp.sh b/utils/misc/package_gmsh_getdp.sh
index 0be07843cad2f7f4c90a767563a45985fdbbdd24..710e986683d37ce230a0680ef8e55d6c54756e81 100755
--- a/utils/misc/package_gmsh_getdp.sh
+++ b/utils/misc/package_gmsh_getdp.sh
@@ -20,7 +20,7 @@ GETDP=svn
 #GMSH=2.8.5
 #GETDP=2.4.4
 
-MODELS='machines relay inductor indheat magnetometer antennas acoustic_scattering time_reversal shielding waveguides transfo_simple ddm_wave_simple bloch_periodic_waveguides magnets thermal_conduction'
+MODELS='machines relay inductor indheat magnetometer antennas acoustic_scattering time_reversal shielding waveguides transfo_simple ddm_wave_simple bloch_periodic_waveguides magnets thermal_conduction magnetostriction'
 
 # get onelab models
 mkdir /tmp/models