diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index 892c4f65a8faa7155822abcf1bc02dd4c967a61c..1f3ce8f2301004b1fb3e9f867b8bcbb765a17a88 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -138,7 +138,7 @@ int GmshOpenProject(const std::string &fileName)
   return 1;
 }
 
-int GmshWriteFile(std::string fileName)
+int GmshWriteFile(const std::string &fileName)
 {
   CreateOutputFile(fileName, FORMAT_AUTO);
   return 1;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index ff5a4dc6a7efd96d353f58adf5eadee83c46d708..da6213b2fc07272e4e8609de8109c6732f96ffbe 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -126,12 +126,12 @@ void MElement::scaledJacRange(double &jmin, double &jmax)
 #endif
 }
 
-void MElement::getNode(int num, double &u, double &v, double &w)
+void MElement::getNode(int num, double &u, double &v, double &w) const
 {
   // only for MElements that don't have a lookup table for this
   // (currently only 1st order elements have)
   double uvw[3];
-  MVertex* ver = getVertex(num);
+  const MVertex* ver = getVertex(num);
   double xyz[3] = {ver->x(), ver->y(), ver->z()};
   xyz2uvw(xyz, uvw);
   u = uvw[0];
@@ -169,7 +169,7 @@ void MElement::getThirdDerivativeShapeFunctions(double u, double v, double w,
   else Msg::Error("Function space not implemented for this type of element");
 }
 
-SPoint3 MElement::barycenter_infty ()
+SPoint3 MElement::barycenter_infty () const
 {
   double xmin =  getVertex(0)->x();
   double xmax = xmin;
@@ -179,7 +179,7 @@ SPoint3 MElement::barycenter_infty ()
   double zmax = zmin;
   int n = getNumVertices();
   for(int i = 0; i < n; i++) {
-    MVertex *v = getVertex(i);
+    const MVertex *v = getVertex(i);
     xmin = std::min(xmin,v->x());
     xmax = std::max(xmax,v->x());
     ymin = std::min(ymin,v->y());
@@ -190,12 +190,12 @@ SPoint3 MElement::barycenter_infty ()
   return SPoint3(0.5*(xmin+xmax),0.5*(ymin+ymax),0.5*(zmin+zmax));
 }
 
-SPoint3 MElement::barycenter()
+SPoint3 MElement::barycenter() const
 {
   SPoint3 p(0., 0., 0.);
   int n = getNumVertices();
   for(int i = 0; i < n; i++) {
-    MVertex *v = getVertex(i);
+    const MVertex *v = getVertex(i);
     p[0] += v->x();
     p[1] += v->y();
     p[2] += v->z();
@@ -206,7 +206,7 @@ SPoint3 MElement::barycenter()
   return p;
 }
 
-SPoint3 MElement::barycenterUVW()
+SPoint3 MElement::barycenterUVW() const
 {
   SPoint3 p(0., 0., 0.);
   int n = getNumVertices();
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 8258a1f5ff250669f1b9a95307e076bcf2d8c9c2..c5123c5caae07beb922097646f8c89c253e4c901 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -196,10 +196,10 @@ class MElement
   virtual double getOuterRadius(){ return 0.; }
 
   // compute the barycenter
-  virtual SPoint3 barycenter();
-  virtual SPoint3 barycenterUVW();
+  virtual SPoint3 barycenter() const;
+  virtual SPoint3 barycenterUVW() const;
   // compute the barycenter in infinity norm
-  virtual SPoint3 barycenter_infty();
+  virtual SPoint3 barycenter_infty() const;
 
   // revert the orientation of the element
   virtual void revert(){}
@@ -225,7 +225,7 @@ class MElement
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const { return 0; }
 
   // return parametric coordinates (u,v,w) of a vertex
-  virtual void getNode(int num, double &u, double &v, double &w);
+  virtual void getNode(int num, double &u, double &v, double &w) const;
 
   // return the interpolating nodal shape functions evaluated at point
   // (u,v,w) in parametric coordinates (if order == -1, use the
@@ -286,7 +286,7 @@ class MElement
 
   // test if a point, given in parametric coordinates, belongs to the
   // element
-  virtual bool isInside(double u, double v, double w) = 0;
+  virtual bool isInside(double u, double v, double w) const = 0;
 
   // interpolate the given nodal data (resp. its gradient, curl and
   // divergence) at point (u,v,w) in parametric coordinates
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 3ff0d570b852ed5d3dc17069529cb8be19427a2e..34a24a2673a87a5450d90f7140cd1c9ea22fd094 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -90,7 +90,7 @@ void MPolyhedron::_init()
 
 }
 
-bool MPolyhedron::isInside(double u, double v, double w)
+bool MPolyhedron::isInside(double u, double v, double w) const
 {
   if(!_orig) return false;
   double uvw[3] = {u, v, w};
@@ -255,7 +255,7 @@ void MPolygon::_initVertices()
   }
 }
 
-bool MPolygon::isInside(double u, double v, double w)
+bool MPolygon::isInside(double u, double v, double w) const
 {
   if(!getParent()) return false;
   double uvw[3] = {u, v, w};
@@ -318,13 +318,13 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 
 //----------------------------------- MLineChild ------------------------------
 
-bool MLineChild::isInside(double u, double v, double w)
+bool MLineChild::isInside(double u, double v, double w) const
 {
   if(!_orig) return false;
   double uvw[3] = {u, v, w};
   double v_uvw[2][3];
   for(int i = 0; i < 2; i++) {
-    MVertex *vi = getVertex(i);
+    const MVertex *vi = getVertex(i);
     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
     _orig->xyz2uvw(v_xyz, v_uvw[i]);
   }
@@ -372,13 +372,13 @@ void MLineChild::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 
 //----------------------------------- MTriangleBorder ------------------------------
 
-bool MTriangleBorder::isInside(double u, double v, double w)
+bool MTriangleBorder::isInside(double u, double v, double w) const
 {
   if(!getParent()) return false;
   double uvw[3] = {u, v, w};
   double v_uvw[3][3];
   for(int i = 0; i < 3; i++) {
-    MVertex *vi = getVertex(i);
+    const MVertex *vi = getVertex(i);
     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
     getParent()->xyz2uvw(v_xyz, v_uvw[i]);
   }
@@ -430,13 +430,13 @@ void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 
 //-------------------------------------- MLineBorder ------------------------------
 
-bool MLineBorder::isInside(double u, double v, double w)
+bool MLineBorder::isInside(double u, double v, double w) const
 {
   if(!getParent()) return false;
   double uvw[3] = {u, v, w};
   double v_uvw[2][3];
   for(int i = 0; i < 2; i++) {
-    MVertex *vi = getVertex(i);
+    const MVertex *vi = getVertex(i);
     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
     getParent()->xyz2uvw(v_xyz, v_uvw[i]);
   }
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index d5040073d85ba9b585f7c8e2e804da506bbbe24f..b8be8f6138b3c2a77bf1ef11d2e4361c2277807a 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -160,7 +160,7 @@ class MPolyhedron : public MElement {
 
   // the parametric coordinates of the polyhedron are
   // the coordinates in the local parent element.
-  virtual bool isInside(double u, double v, double w);
+  virtual bool isInside(double u, double v, double w) const;
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
   virtual MElement *getParent() const { return _orig; }
   virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
@@ -310,7 +310,7 @@ class MPolygon : public MElement {
 
   // the parametric coordinates of the polygon are
   // the coordinates in the local parent element.
-  virtual bool isInside(double u, double v, double w);
+  virtual bool isInside(double u, double v, double w) const;
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
   virtual int getNumVerticesForMSH() {return _parts.size() * 3;}
   virtual void getVerticesIdForMSH(std::vector<int> &verts)
@@ -365,7 +365,7 @@ class MLineChild : public MLine {
   }
   // the parametric coordinates of the LineChildren are
   // the coordinates in the local parent element.
-  virtual bool isInside(double u, double v, double w);
+  virtual bool isInside(double u, double v, double w) const;
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
   virtual MElement *getParent() const { return _orig; }
   virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
@@ -401,7 +401,7 @@ class MTriangleBorder : public MTriangle {
     return NULL;
   }
   virtual int getTypeForMSH() const { return MSH_TRI_B; }
-  virtual bool isInside(double u, double v, double w);
+  virtual bool isInside(double u, double v, double w) const;
   // the integration points of the MTriangleBorder are in the parent element space
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
 };
@@ -460,7 +460,7 @@ class MLineBorder : public MLine {
     return NULL;
   }
   virtual int getTypeForMSH() const { return MSH_LIN_B; }
-  virtual bool isInside(double u, double v, double w);
+  virtual bool isInside(double u, double v, double w) const;
   // the integration points of the MLineBorder are in the parent element space
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
 };
diff --git a/Geo/MFace.h b/Geo/MFace.h
index 94569447d5d650241112d6026271f783cf49544c..0f4f095ae872c647e0a588ff1fc5b77cf4bdd4ea 100644
--- a/Geo/MFace.h
+++ b/Geo/MFace.h
@@ -60,7 +60,7 @@ class MFace {
     SPoint3 p(0., 0., 0.);
     int n = getNumVertices();
     for(int i = 0; i < n; i++) {
-      MVertex *v = getVertex(i);
+      const MVertex *v = getVertex(i);
       p[0] += v->x();
       p[1] += v->y();
       p[2] += v->z();
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index 37ea4605a073438fe8e073e8883632b9d1f0c089..ef433ca04e2467b159ba63c872416d14086063fa 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -128,7 +128,7 @@ class MHexahedron : public MElement {
     tmp = _v[4]; _v[4] = _v[6]; _v[6] = tmp;
   }
   virtual int getVolumeSign();
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     switch(num) {
     case 0 : u = -1.; v = -1.; w = -1.; break;
@@ -142,11 +142,11 @@ class MHexahedron : public MElement {
     default: u =  0.; v =  0.; w =  0.; break;
     }
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(0., 0., 0.);
   }
-  virtual bool isInside(double u, double v, double w)
+  virtual bool isInside(double u, double v, double w) const
   {
     double tol = _isInsideTolerance;
     if(u < -(1. + tol) || v < -(1. + tol) || w < -(1. + tol) ||
@@ -326,7 +326,7 @@ class MHexahedron20 : public MHexahedron {
     _vs[8] = old[10]; _vs[10] = old[8];
     _vs[9] = old[11]; _vs[11] = old[9];
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
@@ -471,7 +471,7 @@ class MHexahedron27 : public MHexahedron {
     _vs[13] = old[15]; _vs[15] = old[13];
     _vs[14] = old[16]; _vs[16] = old[14];
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
@@ -606,7 +606,7 @@ class MHexahedronN : public MHexahedron {
   {
     Msg::Error("Revert not implemented yet for MHexahedronN");
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 90354c9d9d4155af8256e5403045fef138ec4f13..f4c55d834cb8d94ca9e43f0025496d25219b0082 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -76,14 +76,14 @@ class MLine : public MElement {
   }
   virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
-  virtual bool isInside(double u, double v, double w)
+  virtual bool isInside(double u, double v, double w) const
   {
     double tol = _isInsideTolerance;
     if(u < -(1. + tol) || u > (1. + tol) || fabs(v) > tol || fabs(w) > tol)
       return false;
     return true;
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     v = w = 0.;
     switch(num) {
@@ -92,7 +92,7 @@ class MLine : public MElement {
     default: u =  0.; break;
     }
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(0, 0, 0);
   }
@@ -152,7 +152,7 @@ class MLine3 : public MLine {
   virtual int getTypeForVTK() const { return 21; }
   virtual const char *getStringForPOS() const { return "SL2"; }
   virtual const char *getStringForINP() const { return "C1D3"; }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 2 ? MLine::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
@@ -215,7 +215,7 @@ class MLineN : public MLine {
     if(_vs.size() == 9) return MSH_LIN_11;
     return 0;
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 2 ? MLine::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 1bfbdcdd9e938b6710235df5f96643c0e0aa5ad3..66ae8582c12e137e95b91cb6c11905dc9f277d8e 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -45,11 +45,11 @@ class MPoint : public MElement {
   virtual int getTypeForMSH() const { return MSH_PNT; }
   virtual int getTypeForVTK() const { return 1; }
   virtual const char *getStringForPOS() const { return "SP"; }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     u = v = w = 0.;
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(0., 0., 0.);
   }
@@ -69,7 +69,7 @@ class MPoint : public MElement {
   {
     return JacobianBasis::find(MSH_PNT);
   }
-  virtual bool isInside(double u, double v, double w)
+  virtual bool isInside(double u, double v, double w) const
   {
     double tol = _isInsideTolerance;
     if(fabs(u) > tol || fabs(v) > tol || fabs(w) > tol)
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 8a01a87900ef829a5140bb28705eb7c52e3a8988..fbdfcf78cc8be7ab4f04192766e0ef9d68df1716 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -131,7 +131,7 @@ class MPrism : public MElement {
   virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual int getVolumeSign();
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     switch(num) {
     case 0 : u = 0.; v = 0.; w = -1.; break;
@@ -143,11 +143,11 @@ class MPrism : public MElement {
     default: u = 0.; v = 0.; w =  0.; break;
     }
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(1/3., 1/3., 0.);
   }
-  virtual bool isInside(double u, double v, double w)
+  virtual bool isInside(double u, double v, double w) const
   {
     double tol = _isInsideTolerance;
     if(w > (1. + tol) || w < -(1. + tol) || u < (-tol)
@@ -307,7 +307,7 @@ class MPrism15 : public MPrism {
     tmp = _vs[2]; _vs[2] = _vs[4]; _vs[4] = tmp;
     tmp = _vs[7]; _vs[7] = _vs[8]; _vs[8] = tmp;
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 6 ? MPrism::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
@@ -434,7 +434,7 @@ class MPrism18 : public MPrism {
     // quad face vertices
     tmp = _vs[10]; _vs[10] = _vs[11]; _vs[11] = tmp;
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 6 ? MPrism::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 227a532c0b5eb851ccb88c36c987ea2272a72d88..540f73fb2d65f69d5a147922e574276b2d69ddc4 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -126,7 +126,7 @@ class MPyramid : public MElement {
     MVertex *tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
   }
   virtual int getVolumeSign();
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     switch(num) {
     case 0 : u = -1.; v = -1.; w = 0.; break;
@@ -137,11 +137,11 @@ class MPyramid : public MElement {
     default: u =  0.; v =  0.; w = 0.; break;
     }
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(0., 0., .2);
   }
-  virtual bool isInside(double u, double v, double w)
+  virtual bool isInside(double u, double v, double w) const
   {
     double tol = _isInsideTolerance;
     if(u < (w - (1. + tol)) || u > ((1. + tol) - w) || v < (w - (1. + tol)) ||
@@ -325,7 +325,7 @@ class MPyramidN : public MPyramid {
   virtual int getNumEdgesRep();
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
   virtual int getNumFacesRep();
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 5 ? MPyramid::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 2a51c1a14f7a3cea2a90817f173780d2835f5fa8..1422d1ab8456b39e3e7b165539aa6632e20385b7 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -119,7 +119,7 @@ class MQuadrangle : public MElement {
   virtual const char *getStringForINP() const { return "C2D4"; }
   virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     w = 0.;
     switch(num) {
@@ -130,7 +130,7 @@ class MQuadrangle : public MElement {
     default: u =  0.; v =  0.; break;
     }
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(0., 0., 0.);
   }
@@ -139,7 +139,7 @@ class MQuadrangle : public MElement {
   {
     MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
   }
-  virtual bool isInside(double u, double v, double w)
+  virtual bool isInside(double u, double v, double w) const
   {
     double tol = _isInsideTolerance;
     if(u < -(1. + tol) || v < -(1. + tol) || u > (1. + tol) || v > (1. + tol) ||
@@ -244,11 +244,11 @@ class MQuadrangle8 : public MQuadrangle {
     tmp = _vs[0]; _vs[0] = _vs[3]; _vs[3] = tmp;
     tmp = _vs[1]; _vs[1] = _vs[2]; _vs[2] = tmp;
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(0., 0., 0.);
   }
@@ -325,11 +325,11 @@ class MQuadrangle9 : public MQuadrangle {
     tmp = _vs[0]; _vs[0] = _vs[3]; _vs[3] = tmp;
     tmp = _vs[1]; _vs[1] = _vs[2]; _vs[2] = tmp;
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(0., 0., 0.);
   }
@@ -426,11 +426,11 @@ class MQuadrangleN : public MQuadrangle {
     inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
     _vs = inv;
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(0., 0., 0.);
   }
diff --git a/Geo/MSubElement.cpp b/Geo/MSubElement.cpp
index 5ad85d65154ead05bb0d81eb7b22f17cc16590d7..5c26f71ada6c83a8feb2a55fe8a2c1d693d2071a 100644
--- a/Geo/MSubElement.cpp
+++ b/Geo/MSubElement.cpp
@@ -99,7 +99,7 @@ void MSubTetrahedron::xyz2uvw(double xyz[3], double uvw[3]) const
   if(_orig) _orig->xyz2uvw(xyz,uvw);
 }
 
-void MSubTetrahedron::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
+void MSubTetrahedron::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w) const
 {
   if(!_orig) return;
   SPoint3 p;
@@ -110,7 +110,7 @@ void MSubTetrahedron::movePointFromParentSpaceToElementSpace(double &u, double &
   u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
 }
 
-void MSubTetrahedron::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
+void MSubTetrahedron::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w) const
 {
   if(!_orig) return;
   SPoint3 p;
@@ -121,7 +121,7 @@ void MSubTetrahedron::movePointFromElementSpaceToParentSpace(double &u, double &
   u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
 }
 
-bool MSubTetrahedron::isInside(double u, double v, double w)
+bool MSubTetrahedron::isInside(double u, double v, double w) const
 {
   if(!_orig) return false;
 
@@ -338,7 +338,7 @@ void MSubTriangle::xyz2uvw(double xyz[3], double uvw[3]) const
   if(_orig) _orig->xyz2uvw(xyz,uvw);
 }
 
-void MSubTriangle::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
+void MSubTriangle::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w) const
 {
   if(!_orig) return;
   SPoint3 p;
@@ -349,7 +349,7 @@ void MSubTriangle::movePointFromParentSpaceToElementSpace(double &u, double &v,
   u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
 }
 
-void MSubTriangle::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
+void MSubTriangle::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w) const
 {
   if(!_orig) return;
   SPoint3 p;
@@ -360,7 +360,7 @@ void MSubTriangle::movePointFromElementSpaceToParentSpace(double &u, double &v,
   u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
 }
 
-bool MSubTriangle::isInside(double u, double v, double w)
+bool MSubTriangle::isInside(double u, double v, double w) const
 {
   if(!_orig) return false;
 
@@ -566,7 +566,7 @@ void MSubLine::xyz2uvw(double xyz[3], double uvw[3]) const
   if(_orig) _orig->xyz2uvw(xyz,uvw);
 }
 
-void MSubLine::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
+void MSubLine::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w) const
 {
   if(!_orig) return;
   SPoint3 p;
@@ -577,7 +577,7 @@ void MSubLine::movePointFromParentSpaceToElementSpace(double &u, double &v, doub
   u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
 }
 
-void MSubLine::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
+void MSubLine::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w) const
 {
   if(!_orig) return;
   SPoint3 p;
@@ -588,7 +588,7 @@ void MSubLine::movePointFromElementSpaceToParentSpace(double &u, double &v, doub
   u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
 }
 
-bool MSubLine::isInside(double u, double v, double w)
+bool MSubLine::isInside(double u, double v, double w) const
 {
   if(!_orig) return false;
 
@@ -758,7 +758,7 @@ void MSubPoint::xyz2uvw(double xyz[3], double uvw[3]) const
   if(_orig) _orig->xyz2uvw(xyz,uvw);
 }
 
-void MSubPoint::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
+void MSubPoint::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w) const
 {
   if(!_orig) return;
   SPoint3 p;
@@ -769,7 +769,7 @@ void MSubPoint::movePointFromParentSpaceToElementSpace(double &u, double &v, dou
   u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
 }
 
-void MSubPoint::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
+void MSubPoint::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w) const
 {
   if(!_orig) return;
   SPoint3 p;
@@ -780,7 +780,7 @@ void MSubPoint::movePointFromElementSpaceToParentSpace(double &u, double &v, dou
   u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
 }
 
-bool MSubPoint::isInside(double u, double v, double w)
+bool MSubPoint::isInside(double u, double v, double w) const
 {
   if(!_orig) return false;
 
diff --git a/Geo/MSubElement.h b/Geo/MSubElement.h
index ee379b483ecd6bbaf352952e986b99b486d26d81..34b05983991a16291c78429afb3f52c84533a406 100644
--- a/Geo/MSubElement.h
+++ b/Geo/MSubElement.h
@@ -57,9 +57,9 @@ class MSubTetrahedron : public MTetrahedron
   virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
-  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
-  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
-  virtual bool isInside(double u, double v, double w);
+  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w) const;
+  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w) const;
+  virtual bool isInside(double u, double v, double w) const;
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
 
   virtual MElement *getParent() const { return _orig; }
@@ -115,9 +115,9 @@ class MSubTriangle : public MTriangle
   virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
-  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
-  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
-  virtual bool isInside(double u, double v, double w);
+  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w) const;
+  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w) const;
+  virtual bool isInside(double u, double v, double w) const;
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
 
   virtual MElement *getParent() const { return _orig; }
@@ -173,9 +173,9 @@ class MSubLine : public MLine
   virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
-  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
-  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
-  virtual bool isInside(double u, double v, double w);
+  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w) const;
+  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w) const;
+  virtual bool isInside(double u, double v, double w) const;
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
 
   virtual MElement *getParent() const { return _orig; }
@@ -231,9 +231,9 @@ class MSubPoint : public MPoint
   virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
-  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
-  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
-  virtual bool isInside(double u, double v, double w);
+  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w) const;
+  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w) const;
+  virtual bool isInside(double u, double v, double w) const;
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
 
   virtual MElement *getParent() const { return _orig; }
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 220e4aaae381d8657c693b694200e1d3833c6b2a..264a6b0bd5ca8c41a25bc1630e42daad6c59c027 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -131,7 +131,7 @@ class MTetrahedron : public MElement {
   void xyz2uvw(double xyz[3], double uvw[3]) const;
   virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     switch(num) {
     case 0 : u = 0.; v = 0.; w = 0.; break;
@@ -141,11 +141,11 @@ class MTetrahedron : public MElement {
     default: u = 0.; v = 0.; w = 0.; break;
     }
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(.25, .25, .25);
   }
-  virtual bool isInside(double u, double v, double w)
+  virtual bool isInside(double u, double v, double w) const
   {
     double tol = _isInsideTolerance;
     if(u < (-tol) || v < (-tol) || w < (-tol) || u > ((1. + tol) - v - w))
@@ -269,7 +269,7 @@ class MTetrahedron10 : public MTetrahedron {
     tmp = _vs[1]; _vs[1] = _vs[2]; _vs[2] = tmp;
     tmp = _vs[5]; _vs[5] = _vs[3]; _vs[3] = tmp;
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 4 ? MTetrahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
@@ -399,7 +399,7 @@ class MTetrahedronN : public MTetrahedron {
   virtual int getNumEdgesRep();
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
   virtual int getNumFacesRep();
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 4 ? MTetrahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 73b351f5ada0e2ba363841b5861a1101233dfffd..ea6e66bc647563134cd09c761348a1230712957f 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -126,7 +126,7 @@ class MTriangle : public MElement {
   }
   virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     w = 0.;
     switch(num) {
@@ -136,11 +136,11 @@ class MTriangle : public MElement {
     default: u = 0.; v = 0.; break;
     }
   }
-  virtual SPoint3 barycenterUVW()
+  virtual SPoint3 barycenterUVW() const
   {
     return SPoint3(1/3., 1/3., 0.);
   }
-  virtual bool isInside(double u, double v, double w)
+  virtual bool isInside(double u, double v, double w) const
   {
     double tol = _isInsideTolerance;
     if(u < (-tol) || v < (-tol) || u > ((1. + tol) - v) || fabs(w) > tol)
@@ -233,7 +233,7 @@ class MTriangle6 : public MTriangle {
     tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
     tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp;
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
@@ -342,7 +342,7 @@ class MTriangleN : public MTriangle {
     inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
     _vs = inv;
   }
-  virtual void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w) const
   {
     num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }