From 208a1cbc32f9184709c8a499f120e75215924fd7 Mon Sep 17 00:00:00 2001
From: Gauthier Becker <gauthierbecker@gmail.com>
Date: Mon, 18 Mar 2013 00:37:52 +0000
Subject: [PATCH] Declare/Define as const some functions of MElement that are
 const but were not declared/defined like this. I Hope that nothing is broken.
 (Maybe I forget to put const in a child class...). The const is needed for
 compatibility with Summit.

---
 Geo/MElement.cpp     | 27 ++++++------
 Geo/MElement.h       | 34 ++++++++-------
 Geo/MElementCut.h    | 44 ++++++++++++++------
 Geo/MHexahedron.h    |  8 +++-
 Geo/MLine.h          |  5 ++-
 Geo/MPoint.h         |  7 ++--
 Geo/MPrism.h         |  5 ++-
 Geo/MPyramid.h       |  4 +-
 Geo/MQuadrangle.h    |  6 ++-
 Geo/MSubElement.cpp  | 99 ++++++++++++++++++++++++++++----------------
 Geo/MSubElement.h    | 68 ++++++++++++++++--------------
 Geo/MTetrahedron.cpp |  2 +-
 Geo/MTetrahedron.h   |  9 ++--
 Geo/MTriangle.cpp    |  2 +-
 Geo/MTriangle.h      | 11 +++--
 15 files changed, 203 insertions(+), 128 deletions(-)

diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 6ab47463b6..ff5a4dc6a7 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -139,15 +139,14 @@ void MElement::getNode(int num, double &u, double &v, double &w)
   w = uvw[2];
 }
 
-void MElement::getShapeFunctions(double u, double v, double w, double s[], int o)
+void MElement::getShapeFunctions(double u, double v, double w, double s[], int o) const
 {
   const nodalBasis* fs = getFunctionSpace(o);
   if(fs) fs->f(u, v, w, s);
   else Msg::Error("Function space not implemented for this type of element");
 }
 
-void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3],
-                                     int o)
+void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3],int o) const
 {
   const nodalBasis* fs = getFunctionSpace(o);
   if(fs) fs->df(u, v, w, s);
@@ -155,7 +154,7 @@ void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3]
 }
 
 void MElement::getHessShapeFunctions(double u, double v, double w, double s[][3][3],
-                                     int o)
+                                     int o) const
 {
   const nodalBasis* fs = getFunctionSpace(o);
   if(fs) fs->ddf(u, v, w, s);
@@ -163,7 +162,7 @@ void MElement::getHessShapeFunctions(double u, double v, double w, double s[][3]
 }
 
 void MElement::getThirdDerivativeShapeFunctions(double u, double v, double w,
-                                                double s[][3][3][3], int o)
+                                                double s[][3][3][3], int o) const
 {
   const nodalBasis* fs = getFunctionSpace(o);
   if(fs) fs->dddf(u, v, w, s);
@@ -261,7 +260,7 @@ std::string MElement::getInfoString()
   return std::string(tmp);
 }
 
-static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
+static double _computeDeterminantAndRegularize(const MElement *ele, double jac[3][3])
 {
   double dJ = 0;
 
@@ -328,7 +327,7 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
   return dJ;
 }
 
-double MElement::getJacobian(double u, double v, double w, double jac[3][3])
+double MElement::getJacobian(double u, double v, double w, double jac[3][3]) const
 {
   jac[0][0] = jac[0][1] = jac[0][2] = 0.;
   jac[1][0] = jac[1][1] = jac[1][2] = 0.;
@@ -349,7 +348,7 @@ double MElement::getJacobian(double u, double v, double w, double jac[3][3])
   return _computeDeterminantAndRegularize(this, jac);
 }
 
-double MElement::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
+double MElement::getJacobian(const fullMatrix<double> &gsf, double jac[3][3]) const
 {
   jac[0][0] = jac[0][1] = jac[0][2] = 0.;
   jac[1][0] = jac[1][1] = jac[1][2] = 0.;
@@ -367,7 +366,7 @@ double MElement::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
 }
 
 double MElement::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
-{
+const {
   jac[0][0] = jac[0][1] = jac[0][2] = 0.;
   jac[1][0] = jac[1][1] = jac[1][2] = 0.;
   jac[2][0] = jac[2][1] = jac[2][2] = 0.;
@@ -424,7 +423,7 @@ void MElement::getNodesCoord(fullMatrix<double> &nodesXYZ)
   }
 }
 
-void MElement::pnt(double u, double v, double w, SPoint3 &p)
+void MElement::pnt(double u, double v, double w, SPoint3 &p) const
 {
   double x = 0., y = 0., z = 0.;
   double sf[1256];
@@ -438,7 +437,7 @@ void MElement::pnt(double u, double v, double w, SPoint3 &p)
   p = SPoint3(x, y, z);
 }
 
-void MElement::pnt(const std::vector<double> &sf, SPoint3 &p)
+void MElement::pnt(const std::vector<double> &sf, SPoint3 &p) const
 {
   double x = 0., y = 0., z = 0.;
   for (int j = 0; j < getNumShapeFunctions(); j++) {
@@ -464,7 +463,7 @@ void MElement::primaryPnt(double u, double v, double w, SPoint3 &p)
   p = SPoint3(x,y,z);
 }
 
-void MElement::xyz2uvw(double xyz[3], double uvw[3])
+void MElement::xyz2uvw(double xyz[3], double uvw[3]) const
 {
   // general Newton routine for the nonlinear case (more efficient
   // routines are implemented for simplices, where the basis functions
@@ -480,7 +479,7 @@ void MElement::xyz2uvw(double xyz[3], double uvw[3])
     double sf[1256];
     getShapeFunctions(uvw[0], uvw[1], uvw[2], sf);
     for (int i = 0; i < getNumShapeFunctions(); i++) {
-      MVertex *v = getShapeFunctionNode(i);
+      const MVertex *v = getShapeFunctionNode(i);
       xn += v->x() * sf[i];
       yn += v->y() * sf[i];
       zn += v->z() * sf[i];
@@ -501,7 +500,7 @@ void MElement::xyz2uvw(double xyz[3], double uvw[3])
   }
 }
 
-void MElement::xyzTouvw(fullMatrix<double> *xu)
+void MElement::xyzTouvw(fullMatrix<double> *xu) const
 {
   double _xyz[3] = {(*xu)(0,0),(*xu)(0,1),(*xu)(0,2)}, _uvw[3];
   xyz2uvw(_xyz, _uvw);
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 533c877a49..76429591a1 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -68,6 +68,7 @@ class MElement
 
   // get & set the vertices
   virtual int getNumVertices() const = 0;
+  virtual const MVertex *getVertex(int num) const= 0;
   virtual MVertex *getVertex(int num) = 0;
   void getVertices(std::vector<MVertex*> &verts)
   {
@@ -117,7 +118,7 @@ class MElement
 
   // get the edges
   virtual int getNumEdges() = 0;
-  virtual MEdge getEdge(int num) = 0;
+  virtual MEdge getEdge(int num) const= 0;
 
   // give an MEdge as input and get its local number and sign
   virtual void getEdgeInfo(const MEdge & edge, int &ithEdge, int &sign) const
@@ -230,24 +231,24 @@ class MElement
   // (u,v,w) in parametric coordinates (if order == -1, use the
   // polynomial order of the element)
   virtual void getShapeFunctions(double u, double v, double w, double s[],
-                                 int order=-1);
+                                 int order=-1) const;
 
   // return the gradient of of the nodal shape functions evaluated at
   // point (u,v,w) in parametric coordinates (if order == -1, use the
   // polynomial order of the element)
   virtual void getGradShapeFunctions(double u, double v, double w, double s[][3],
-                                     int order=-1);
+                                     int order=-1) const;
   virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3],
-                                     int order=-1);
+                                     int order=-1) const;
   virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3],
-                                     int order=-1);
+                                     int order=-1) const;
   // return the Jacobian of the element evaluated at point (u,v,w) in
   // parametric coordinates
-  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
+  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]) const;
   // To be compatible with _vgrads of functionSpace without having to put under fullMatrix form
-  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
-  virtual double getJacobian(double u, double v, double w, double jac[3][3]);
-  inline double getJacobian(double u, double v, double w, fullMatrix<double> &j){
+  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])const ;
+  virtual double getJacobian(double u, double v, double w, double jac[3][3]) const;
+  inline double getJacobian(double u, double v, double w, fullMatrix<double> &j) const{
     double JAC[3][3];
     const double detJ = getJacobian (u,v,w,JAC);
     for (int i=0;i<3;i++){
@@ -264,19 +265,20 @@ class MElement
   }
   void getSignedJacobian(fullVector<double> &jacobian, int o = -1);
   void getNodesCoord(fullMatrix<double> &nodesXYZ);
-  virtual int getNumShapeFunctions(){ return getNumVertices(); }
-  virtual int getNumPrimaryShapeFunctions(){ return getNumPrimaryVertices(); }
-  virtual MVertex *getShapeFunctionNode(int i){ return getVertex(i); }
+  virtual int getNumShapeFunctions() const{ return getNumVertices(); }
+  virtual int getNumPrimaryShapeFunctions() { return getNumPrimaryVertices(); }
+  virtual const MVertex *getShapeFunctionNode(int i) const{ return getVertex(i); }
+  virtual MVertex *getShapeFunctionNode(int i) { return getVertex(i); }
 
   // get the point in cartesian coordinates corresponding to the point
   // (u,v,w) in parametric coordinates
-  virtual void pnt(double u, double v, double w, SPoint3 &p);
-  virtual void pnt(const std::vector<double> &sf,SPoint3 &p); // To be compatible with functionSpace without changing form
+  virtual void pnt(double u, double v, double w, SPoint3 &p) const;
+  virtual void pnt(const std::vector<double> &sf,SPoint3 &p) const; // To be compatible with functionSpace without changing form
   virtual void primaryPnt(double u, double v, double w, SPoint3 &p);
 
   // invert the parametrisation
-  virtual void xyz2uvw(double xyz[3], double uvw[3]);
-  void xyzTouvw(fullMatrix<double> *xu);
+  virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
+  void xyzTouvw(fullMatrix<double> *xu) const;
 
   // move point between parent and element parametric spaces
   virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index 73e03e27bd..d5040073d8 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -66,8 +66,12 @@ class MPolyhedron : public MElement {
     return (num < (int)_vertices.size()) ?
       _vertices[num] : _innerVertices[num - _vertices.size()];
   }
+  virtual const MVertex *getVertex(int num) const
+  {
+    return (num < (int)_vertices.size()) ? _vertices[num] : _innerVertices[num - _vertices.size()];
+  }
   virtual int getNumEdges() { return _edges.size(); }
-  virtual MEdge getEdge(int num) { return _edges[num]; }
+  virtual MEdge getEdge(int num) const{ return _edges[num]; }
   virtual int getNumEdgesRep() { return _edges.size(); }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
   {
@@ -125,19 +129,19 @@ class MPolyhedron : public MElement {
   {
     return (_orig ? _orig->getJacobianFuncSpace(order) : 0);
   }
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int o) const
   {
     if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
   }
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o) const
   {
     if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
   }
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o) const
   {
     if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
   }
-  virtual int getNumShapeFunctions()
+  virtual int getNumShapeFunctions() const
   {
     return (_orig ? _orig->getNumShapeFunctions() : 0);
   }
@@ -145,10 +149,15 @@ class MPolyhedron : public MElement {
   {
     return (_orig ? _orig->getNumPrimaryShapeFunctions() : 0);
   }
+  virtual const MVertex *getShapeFunctionNode(int i) const
+  {
+    return (_orig ? _orig->getShapeFunctionNode(i) : 0);
+  }
   virtual MVertex *getShapeFunctionNode(int i)
   {
     return (_orig ? _orig->getShapeFunctionNode(i) : 0);
   }
+
   // the parametric coordinates of the polyhedron are
   // the coordinates in the local parent element.
   virtual bool isInside(double u, double v, double w);
@@ -214,8 +223,12 @@ class MPolygon : public MElement {
     return (num < (int)_vertices.size()) ?
       _vertices[num] : _innerVertices[num - _vertices.size()];
   }
+  virtual const MVertex *getVertex(int num) const
+  {
+    return (num < (int)_vertices.size()) ? _vertices[num] : _innerVertices[num - _vertices.size()];
+  }
   virtual int getNumEdges() { return _edges.size(); }
-  virtual MEdge getEdge(int num) { return _edges[num]; }
+  virtual MEdge getEdge(int num) const{ return _edges[num]; }
   virtual int getNumEdgesRep() { return getNumEdges(); }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
   {
@@ -266,19 +279,19 @@ class MPolygon : public MElement {
   {
     return (_orig ? _orig->getJacobianFuncSpace(order) : 0);
   }
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int o) const
   {
     if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
   }
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o) const
   {
     if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
   }
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o) const
   {
     if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
   }
-  virtual int getNumShapeFunctions()
+  virtual int getNumShapeFunctions() const
   {
     return (_orig ? _orig->getNumShapeFunctions() : 0);
   }
@@ -286,10 +299,15 @@ class MPolygon : public MElement {
   {
     return (_orig ? _orig->getNumPrimaryShapeFunctions() : 0);
   }
+  virtual const MVertex *getShapeFunctionNode(int i) const
+  {
+    return (_orig ? _orig->getShapeFunctionNode(i) : 0);
+  }
   virtual MVertex *getShapeFunctionNode(int i)
   {
     return (_orig ? _orig->getShapeFunctionNode(i) : 0);
   }
+
   // the parametric coordinates of the polygon are
   // the coordinates in the local parent element.
   virtual bool isInside(double u, double v, double w);
@@ -333,15 +351,15 @@ class MLineChild : public MLine {
     if(_orig) return _orig->getJacobianFuncSpace(order);
     return 0;
   }
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int o) const
   {
     if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
   }
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o) const
   {
     if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
   }
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o) const
   {
     if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
   }
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index 7abda6995d..37ea4605a0 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -57,6 +57,7 @@ class MHexahedron : public MElement {
   virtual int getDim() const { return 3; }
   virtual int getNumVertices() const { return 8; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual const MVertex *getVertex(int num) const { return _v[num]; }  
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
   virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
@@ -66,7 +67,7 @@ class MHexahedron : public MElement {
     return getVertex(map[num]);
   }
   virtual int getNumEdges(){ return 12; }
-  virtual MEdge getEdge(int num)
+  virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_hexa(num, 0)], _v[edges_hexa(num, 1)]);
   }
@@ -228,6 +229,7 @@ class MHexahedron20 : public MHexahedron {
   virtual int getPolynomialOrder() const { return 2; }
   virtual int getNumVertices() const { return 20; }
   virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
+  virtual const MVertex *getVertex(int num) const { return num < 8 ? _v[num] : _vs[num - 8]; }  
   virtual MVertex *getVertexUNV(int num)
   {
     static const int map[20] = {0, 8, 1, 11, 2, 13, 3, 9, 10, 12,
@@ -374,6 +376,8 @@ class MHexahedron27 : public MHexahedron {
   virtual int getPolynomialOrder() const { return 2; }
   virtual int getNumVertices() const { return 27; }
   virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
+  virtual const MVertex *getVertex(int num) const { return num < 8 ? _v[num] : _vs[num - 8]; }
+
   virtual MVertex *getVertexDIFF(int num)
   {
     static const int map[27] = {6, 8, 26, 24, 0, 2, 20, 18, 7, 15, 3, 17, 5, 25,
@@ -524,6 +528,8 @@ class MHexahedronN : public MHexahedron {
   virtual int getPolynomialOrder() const { return (int)_order; }
   virtual int getNumVertices() const { return 8 + _vs.size(); }
   virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
+  virtual const MVertex *getVertex(int num) const { return num < 8 ? _v[num] : _vs[num - 8]; }
+
   virtual int getNumEdgeVertices() const { return 12 * (_order - 1); }
   virtual int getNumFaceVertices() const { return 6 * (_order - 1)*(_order - 1); }
   virtual int getNumVolumeVertices() const
diff --git a/Geo/MLine.h b/Geo/MLine.h
index d548a13882..90354c9d9d 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -38,6 +38,7 @@ class MLine : public MElement {
   virtual int getDim() const { return 1; }
   virtual int getNumVertices() const { return 2; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual const MVertex *getVertex(int num) const{ return _v[num]; }  
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
   virtual double getInnerRadius(); // half-length of segment line
   virtual double getLength(); // length of segment line
@@ -47,7 +48,7 @@ class MLine : public MElement {
     ithVertex = _v[0] == vertex ? 0 : 1;
   }
   virtual int getNumEdges(){ return 1; }
-  virtual MEdge getEdge(int num){ return MEdge(_v[0], _v[1]); }
+  virtual MEdge getEdge(int num) const{ return MEdge(_v[0], _v[1]); }
   virtual int getNumEdgesRep(){ return 1; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
   {
@@ -124,6 +125,7 @@ class MLine3 : public MLine {
   virtual int getPolynomialOrder() const { return 2; }
   virtual int getNumVertices() const { return 3; }
   virtual MVertex *getVertex(int num){ return num < 2 ? _v[num] : _vs[num - 2]; }
+  virtual const MVertex *getVertex(int num) const{ return num < 2 ? _v[num] : _vs[num - 2]; }  
   virtual MVertex *getVertexUNV(int num)
   {
     static const int map[3] = {0, 2, 1};
@@ -184,6 +186,7 @@ class MLineN : public MLine {
   virtual int getPolynomialOrder() const { return _vs.size() + 1; }
   virtual int getNumVertices() const { return _vs.size() + 2; }
   virtual MVertex *getVertex(int num){ return num < 2 ? _v[num] : _vs[num - 2]; }
+  virtual const MVertex *getVertex(int num) const{ return num < 2 ? _v[num] : _vs[num - 2]; }
   virtual int getNumEdgeVertices() const { return _vs.size(); }
   virtual int getNumEdgesRep(){ return _vs.size() + 1; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 75a6dfa20d..1bfbdcdd9e 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -31,9 +31,10 @@ class MPoint : public MElement {
   virtual int getDim() const { return 0; }
   virtual int getNumVertices() const { return 1; }
   virtual MVertex *getVertex(int num){ return _v[0]; }
+  virtual const MVertex *getVertex(int num) const { return _v[0]; }
   virtual void setVertex(int num,  MVertex *v){ _v[0] = v; }
   virtual int getNumEdges(){ return 0; }
-  virtual MEdge getEdge(int num){ return MEdge(); }
+  virtual MEdge getEdge(int num) const{ return MEdge(); }
   virtual int getNumEdgesRep(){ return 0; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n){}
   virtual int getNumFaces(){ return 0; }
@@ -52,11 +53,11 @@ class MPoint : public MElement {
   {
     return SPoint3(0., 0., 0.);
   }
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int o) const
   {
     s[0] = 1.;
   }
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o) const
   {
     s[0][0] = s[0][1] = s[0][2] = 0.;
   }
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index ad0a5ee529..8a01a87900 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -64,9 +64,10 @@ class MPrism : public MElement {
   virtual int getNumVertices() const { return 6; }
   virtual double getInnerRadius();
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual const MVertex *getVertex(int num)const{ return _v[num]; }  
   virtual void setVertex(int num,   MVertex *v){ _v[num] = v; }
   virtual int getNumEdges(){ return 9; }
-  virtual MEdge getEdge(int num)
+  virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_prism(num, 0)], _v[edges_prism(num, 1)]);
   }
@@ -227,6 +228,7 @@ class MPrism15 : public MPrism {
   virtual int getPolynomialOrder() const { return 2; }
   virtual int getNumVertices() const { return 15; }
   virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
+  virtual const MVertex *getVertex(int num) const { return num < 6 ? _v[num] : _vs[num - 6]; }  
   virtual MVertex *getVertexUNV(int num)
   {
     static const int map[15] = {0, 6, 1, 9, 2, 7, 8, 10, 11, 3, 12, 4, 14, 5, 13};
@@ -356,6 +358,7 @@ class MPrism18 : public MPrism {
   virtual int getPolynomialOrder() const { return 2; }
   virtual int getNumVertices() const { return 18; }
   virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
+  virtual const MVertex *getVertex(int num) const{ return num < 6 ? _v[num] : _vs[num - 6]; }  
   virtual int getNumEdgeVertices() const { return 9; }
   virtual int getNumFaceVertices() const { return 3; }
   virtual int getNumEdgesRep(){ return 18; }
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index c5111ea20f..227a532c0b 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -67,11 +67,12 @@ class MPyramid : public MElement {
   virtual int getDim() const { return 3; }
   virtual int getNumVertices() const { return 5; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual const MVertex *getVertex(int num) const{ return _v[num]; }  
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
   virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual int getNumEdges(){ return 8; }
-  virtual MEdge getEdge(int num)
+  virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_pyramid(num, 0)], _v[edges_pyramid(num, 1)]);
   }
@@ -244,6 +245,7 @@ class MPyramidN : public MPyramid {
   virtual int getPolynomialOrder() const { return _order; }
   virtual int getNumVertices() const { return 5 + _vs.size(); }
   virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
+  virtual const MVertex *getVertex(int num) const{ return num < 5 ? _v[num] : _vs[num - 5]; }
   virtual int getNumEdgeVertices() const { return 8 * (_order - 1); }
   virtual int getNumFaceVertices() const
   {
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index bd6ee6056c..2a51c1a14f 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -56,6 +56,7 @@ class MQuadrangle : public MElement {
   virtual int getDim() const { return 2; }
   virtual int getNumVertices() const { return 4; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual const MVertex *getVertex(int num) const{ return _v[num]; }  
   virtual void setVertex(int num, MVertex *v){ _v[num] = v; }
   virtual MVertex *getVertexDIFF(int num)
   {
@@ -63,7 +64,7 @@ class MQuadrangle : public MElement {
     return getVertex(map[num]);
   }
   virtual int getNumEdges(){ return 4; }
-  virtual MEdge getEdge(int num)
+  virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_quad(num, 0)], _v[edges_quad(num, 1)]);
   }
@@ -199,6 +200,7 @@ class MQuadrangle8 : public MQuadrangle {
   virtual int getPolynomialOrder() const { return 2; }
   virtual int getNumVertices() const { return 8; }
   virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+  virtual const MVertex *getVertex(int num) const { return num < 4 ? _v[num] : _vs[num - 4]; }  
   virtual MVertex *getVertexUNV(int num)
   {
     static const int map[8] = {0, 4, 1, 5, 2, 6, 3, 7};
@@ -285,6 +287,7 @@ class MQuadrangle9 : public MQuadrangle {
   virtual int getPolynomialOrder() const { return 2; }
   virtual int getNumVertices() const { return 9; }
   virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+  virtual const MVertex *getVertex(int num) const { return num < 4 ? _v[num] : _vs[num - 4]; }
   virtual MVertex *getVertexDIFF(int num)
   {
     static const int map[9] = {0, 2, 8, 6, 1, 5, 7, 3, 4};
@@ -369,6 +372,7 @@ class MQuadrangleN : public MQuadrangle {
   virtual int getPolynomialOrder() const { return _order; }
   virtual int getNumVertices() const {return 4 + _vs.size(); }
   virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+  virtual const MVertex *getVertex(int num) const{ return num < 4 ? _v[num] : _vs[num - 4]; }
   virtual int getNumFaceVertices() const
   {
     if(_order > 1 && (int)_vs.size() + 4 == (_order + 1) * (_order + 1))
diff --git a/Geo/MSubElement.cpp b/Geo/MSubElement.cpp
index 11f82cd009..9a6c1e24b5 100644
--- a/Geo/MSubElement.cpp
+++ b/Geo/MSubElement.cpp
@@ -29,37 +29,37 @@ const JacobianBasis* MSubTetrahedron::getJacobianFuncSpace(int order) const
   return 0;
 }
 
-void MSubTetrahedron::getShapeFunctions(double u, double v, double w, double s[], int order)
+void MSubTetrahedron::getShapeFunctions(double u, double v, double w, double s[], int order) const
 {
   if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
 }
 
-void MSubTetrahedron::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
+void MSubTetrahedron::getGradShapeFunctions(double u, double v, double w, double s[][3], int order) const
 {
   if(_orig) _orig->getGradShapeFunctions(u, v, w, s, order);
 }
 
-void MSubTetrahedron::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
+void MSubTetrahedron::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order) const
 {
   if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
 }
 
-void MSubTetrahedron::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
+void MSubTetrahedron::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order) const
 {
   if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
 }
 
-double MSubTetrahedron::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
+double MSubTetrahedron::getJacobian(const fullMatrix<double> &gsf, double jac[3][3]) const
 {
   if(_orig) return _orig->getJacobian(gsf, jac);
   return 0;
 }
-double MSubTetrahedron::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
+double MSubTetrahedron::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]) const
 {
   if(_orig) return _orig->getJacobian(gsf, jac);
   return 0;
 }
-double MSubTetrahedron::getJacobian(double u, double v, double w, double jac[3][3])
+double MSubTetrahedron::getJacobian(double u, double v, double w, double jac[3][3]) const
 {
   if(_orig) return _orig->getJacobian(u, v, w, jac);
   return 0;
@@ -70,7 +70,7 @@ double MSubTetrahedron::getPrimaryJacobian(double u, double v, double w, double
   return 0;
 }
 
-int MSubTetrahedron::getNumShapeFunctions()
+int MSubTetrahedron::getNumShapeFunctions() const
 {
   if(_orig) return _orig->getNumShapeFunctions();
   return 0;
@@ -82,12 +82,18 @@ int MSubTetrahedron::getNumPrimaryShapeFunctions()
   return 0;
 }
 
-MVertex* MSubTetrahedron::getShapeFunctionNode(int i)
+const MVertex* MSubTetrahedron::getShapeFunctionNode(int i) const
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
   return 0;
 }
 
+MVertex* MSubTetrahedron::getShapeFunctionNode(int i)
+{
+	  if(_orig) return _orig->getShapeFunctionNode(i);
+	    return 0;
+}
+
 void MSubTetrahedron::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
 {
   if(!_orig) return;
@@ -209,12 +215,12 @@ const JacobianBasis* MSubTriangle::getJacobianFuncSpace(int order) const
   return 0;
 }
 
-void MSubTriangle::getShapeFunctions(double u, double v, double w, double s[], int order)
+void MSubTriangle::getShapeFunctions(double u, double v, double w, double s[], int order) const
 {
   if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
 }
 
-void MSubTriangle::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
+void MSubTriangle::getGradShapeFunctions(double u, double v, double w, double s[][3], int order) const
 {
   if(!_orig)
     return;
@@ -266,28 +272,28 @@ void MSubTriangle::getGradShapeFunctions(double u, double v, double w, double s[
   }
 }
 
-void MSubTriangle::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
+void MSubTriangle::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order) const
 {
   if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
 }
 
-void MSubTriangle::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
+void MSubTriangle::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order) const
 {
   if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
 }
 
-double MSubTriangle::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
+double MSubTriangle::getJacobian(const fullMatrix<double> &gsf, double jac[3][3]) const
 {
   if(_orig) return _orig->getJacobian(gsf, jac);
   return 0;
 }
-double MSubTriangle::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
+double MSubTriangle::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]) const
 {
   if(_orig) return _orig->getJacobian(gsf, jac);
   return 0;
 }
 double MSubTriangle::getJacobian(double u, double v, double w, double jac[3][3])
-{
+const{
   if(_orig) return _orig->getJacobian(u, v, w, jac);
   return 0;
 }
@@ -297,7 +303,7 @@ double MSubTriangle::getPrimaryJacobian(double u, double v, double w, double jac
   return 0;
 }
 
-int MSubTriangle::getNumShapeFunctions()
+int MSubTriangle::getNumShapeFunctions() const
 {
   if(_orig) return _orig->getNumShapeFunctions();
   return 0;
@@ -309,12 +315,19 @@ int MSubTriangle::getNumPrimaryShapeFunctions()
   return 0;
 }
 
-MVertex* MSubTriangle::getShapeFunctionNode(int i)
+const MVertex* MSubTriangle::getShapeFunctionNode(int i) const
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
   return 0;
 }
 
+MVertex* MSubTriangle::getShapeFunctionNode(int i)
+{
+	  if(_orig) return _orig->getShapeFunctionNode(i);
+	    return 0;
+}
+
+
 void MSubTriangle::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
 {
   if(!_orig) return;
@@ -434,12 +447,12 @@ const JacobianBasis* MSubLine::getJacobianFuncSpace(int order) const
   return 0;
 }
 
-void MSubLine::getShapeFunctions(double u, double v, double w, double s[], int order)
+void MSubLine::getShapeFunctions(double u, double v, double w, double s[], int order) const
 {
   if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
 }
 
-void MSubLine::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
+void MSubLine::getGradShapeFunctions(double u, double v, double w, double s[][3], int order) const
 {
   if(!_orig)
     return;
@@ -481,28 +494,28 @@ void MSubLine::getGradShapeFunctions(double u, double v, double w, double s[][3]
   }
 }
 
-void MSubLine::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
+void MSubLine::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order) const
 {
   if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
 }
 
-void MSubLine::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
+void MSubLine::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order) const
 {
   if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
 }
 
-double MSubLine::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
+double MSubLine::getJacobian(const fullMatrix<double> &gsf, double jac[3][3]) const
 {
   if(_orig) return _orig->getJacobian(gsf, jac);
   return 0;
 }
-double MSubLine::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
+double MSubLine::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])const
 {
   if(_orig) return _orig->getJacobian(gsf, jac);
   return 0;
 }
 double MSubLine::getJacobian(double u, double v, double w, double jac[3][3])
-{
+const{
   if(_orig) return _orig->getJacobian(u, v, w, jac);
   return 0;
 }
@@ -512,7 +525,7 @@ double MSubLine::getPrimaryJacobian(double u, double v, double w, double jac[3][
   return 0;
 }
 
-int MSubLine::getNumShapeFunctions()
+int MSubLine::getNumShapeFunctions() const
 {
   if(_orig) return _orig->getNumShapeFunctions();
   return 0;
@@ -524,12 +537,19 @@ int MSubLine::getNumPrimaryShapeFunctions()
   return 0;
 }
 
-MVertex* MSubLine::getShapeFunctionNode(int i)
+const MVertex* MSubLine::getShapeFunctionNode(int i) const
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
   return 0;
 }
 
+MVertex* MSubLine::getShapeFunctionNode(int i)
+{
+	  if(_orig) return _orig->getShapeFunctionNode(i);
+	    return 0;
+}
+
+
 void MSubLine::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
 {
   if(!_orig) return;
@@ -649,38 +669,38 @@ const JacobianBasis* MSubPoint::getJacobianFuncSpace(int order) const
   return 0;
 }
 
-void MSubPoint::getShapeFunctions(double u, double v, double w, double s[], int order)
+void MSubPoint::getShapeFunctions(double u, double v, double w, double s[], int order) const
 {
   if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
 }
 
-void MSubPoint::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
+void MSubPoint::getGradShapeFunctions(double u, double v, double w, double s[][3], int order) const
 {
   if(_orig) _orig->getGradShapeFunctions(u, v, w, s, order);
 }
 
-void MSubPoint::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
+void MSubPoint::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order) const
 {
   if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
 }
 
-void MSubPoint::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
+void MSubPoint::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order) const
 {
   if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
 }
 
-double MSubPoint::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
+double MSubPoint::getJacobian(const fullMatrix<double> &gsf, double jac[3][3]) const
 {
   if(_orig) return _orig->getJacobian(gsf, jac);
   return 0;
 }
-double MSubPoint::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
+double MSubPoint::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])const
 {
   if(_orig) return _orig->getJacobian(gsf, jac);
   return 0;
 }
 double MSubPoint::getJacobian(double u, double v, double w, double jac[3][3])
-{
+const{ 
   if(_orig) return _orig->getJacobian(u, v, w, jac);
   return 0;
 }
@@ -690,7 +710,7 @@ double MSubPoint::getPrimaryJacobian(double u, double v, double w, double jac[3]
   return 0;
 }
 
-int MSubPoint::getNumShapeFunctions()
+int MSubPoint::getNumShapeFunctions() const
 {
   if(_orig) return _orig->getNumShapeFunctions();
   return 0;
@@ -702,12 +722,19 @@ int MSubPoint::getNumPrimaryShapeFunctions()
   return 0;
 }
 
-MVertex* MSubPoint::getShapeFunctionNode(int i)
+const MVertex* MSubPoint::getShapeFunctionNode(int i) const
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
   return 0;
 }
 
+MVertex* MSubPoint::getShapeFunctionNode(int i)
+{
+	  if(_orig) return _orig->getShapeFunctionNode(i);
+	    return 0;
+}
+
+
 void MSubPoint::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
 {
   if(!_orig) return;
diff --git a/Geo/MSubElement.h b/Geo/MSubElement.h
index eaf6fef9d5..95de4fcecb 100644
--- a/Geo/MSubElement.h
+++ b/Geo/MSubElement.h
@@ -44,16 +44,17 @@ class MSubTetrahedron : public MTetrahedron
   virtual const nodalBasis* getFunctionSpace(int order=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
   // the parametric coordinates are the coordinates in the local parent element
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
-  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
-  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
-  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
-  virtual double getJacobian(double u, double v, double w, double jac[3][3]);
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1) const;
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1) const;
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1) const;
+  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1) const;
+  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]) const;
+  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])const;
+  virtual double getJacobian(double u, double v, double w, double jac[3][3]) const;
   virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
-  virtual int getNumShapeFunctions();
+  virtual int getNumShapeFunctions() const;
   virtual int getNumPrimaryShapeFunctions();
+  virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
   virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
@@ -100,16 +101,17 @@ class MSubTriangle : public MTriangle
   virtual const nodalBasis* getFunctionSpace(int order=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
   // the parametric coordinates are the coordinates in the local parent element
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
-  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
-  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
-  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
-  virtual double getJacobian(double u, double v, double w, double jac[3][3]);
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1) const;
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1) const;
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1) const;
+  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1) const;
+  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]) const;
+  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])const;
+  virtual double getJacobian(double u, double v, double w, double jac[3][3])const;
   virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
-  virtual int getNumShapeFunctions();
+  virtual int getNumShapeFunctions() const;
   virtual int getNumPrimaryShapeFunctions();
+  virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
   virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
@@ -155,16 +157,17 @@ class MSubLine : public MLine
   virtual const nodalBasis* getFunctionSpace(int order=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
   // the parametric coordinates are the coordinates in the local parent element
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
-  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
-  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
-  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
-  virtual double getJacobian(double u, double v, double w, double jac[3][3]);
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1) const;
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1) const;
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1) const;
+  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1) const;
+  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]) const;
+  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])const;
+  virtual double getJacobian(double u, double v, double w, double jac[3][3]) const;
   virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
-  virtual int getNumShapeFunctions();
+  virtual int getNumShapeFunctions() const;
   virtual int getNumPrimaryShapeFunctions();
+  virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
   virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
@@ -212,16 +215,17 @@ class MSubPoint : public MPoint
   virtual const nodalBasis* getFunctionSpace(int order=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
   // the parametric coordinates are the coordinates in the local parent element
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
-  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
-  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
-  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
-  virtual double getJacobian(double u, double v, double w, double jac[3][3]);
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1) const;
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1) const;
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1) const;
+  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1) const;
+  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]) const;
+  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])const;
+  virtual double getJacobian(double u, double v, double w, double jac[3][3]) const;
   virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
-  virtual int getNumShapeFunctions();
+  virtual int getNumShapeFunctions() const;
   virtual int getNumPrimaryShapeFunctions();
+  virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
   virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index f9adeda896..5d1d29a0b1 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -90,7 +90,7 @@ double MTetrahedron::getVolume()
   return det3x3(mat) / 6.;
 }
 
-void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3])
+void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3]) const
 {
   double mat[3][3], b[3], det;
   getMat(mat);
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 2dbb91c39a..220e4aaae3 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -60,9 +60,10 @@ class MTetrahedron : public MElement {
   virtual int getDim() const { return 3; }
   virtual int getNumVertices() const { return 4; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual const MVertex *getVertex(int num) const { return _v[num]; }  
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
   virtual int getNumEdges(){ return 6; }
-  virtual MEdge getEdge(int num)
+  virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_tetra(num, 0)], _v[edges_tetra(num, 1)]);
   }
@@ -109,7 +110,7 @@ class MTetrahedron : public MElement {
   {
     MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
   }
-  void getMat(double mat[3][3])
+  void getMat(double mat[3][3]) const
   {
     mat[0][0] = _v[1]->x() - _v[0]->x();
     mat[0][1] = _v[2]->x() - _v[0]->x();
@@ -127,7 +128,7 @@ class MTetrahedron : public MElement {
   virtual double getInnerRadius();
   virtual double getCircumRadius();
   virtual double etaShapeMeasure();
-  void xyz2uvw(double xyz[3], double uvw[3]);
+  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)
@@ -216,6 +217,7 @@ class MTetrahedron10 : public MTetrahedron {
   virtual int getPolynomialOrder() const { return 2; }
   virtual int getNumVertices() const { return 10; }
   virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+  virtual const MVertex *getVertex(int num) const { return num < 4 ? _v[num] : _vs[num - 4]; }  
   virtual MVertex *getVertexUNV(int num)
   {
     static const int map[10] = {0, 4, 1, 5, 2, 6, 7, 9, 8, 3};
@@ -330,6 +332,7 @@ class MTetrahedronN : public MTetrahedron {
   virtual int getPolynomialOrder() const { return _order; }
   virtual int getNumVertices() const { return 4 + _vs.size(); }
   virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+  virtual const MVertex *getVertex(int num) const{ return num < 4 ? _v[num] : _vs[num - 4]; }  
   virtual int getNumEdgeVertices() const { return 6 * (_order - 1); }
   virtual int getNumFaceVertices() const
   {
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index bc6c87c4bf..49d68dee7b 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -92,7 +92,7 @@ double MTriangle::gammaShapeMeasure()
 #endif
 }
 
- void MTriangle::xyz2uvw(double xyz[3], double uvw[3])
+ void MTriangle::xyz2uvw(double xyz[3], double uvw[3]) const
  {
    const double O[3] = {_v[0]->x(), _v[0]->y(), _v[0]->z()};
    const double d[3] = {xyz[0] - O[0], xyz[1] - O[1], xyz[2] - O[2]};
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 0adf7c219c..73b351f5ad 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -57,8 +57,9 @@ class MTriangle : public MElement {
   virtual double angleShapeMeasure();
   virtual int getNumVertices() const { return 3; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual const MVertex *getVertex(int num) const { return _v[num]; }
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
-  virtual void xyz2uvw(double xyz[3], double uvw[3]);
+  virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
   virtual MVertex *getOtherVertex(MVertex *v1, MVertex *v2)
   {
     if(_v[0] != v1 && _v[0] != v2) return _v[0];
@@ -67,7 +68,7 @@ class MTriangle : public MElement {
     return 0;
   }
   virtual int getNumEdges(){ return 3; }
-  virtual MEdge getEdge(int num)
+  virtual MEdge getEdge(int num) const
   {
     return MEdge(_v[edges_tri(num, 0)], _v[edges_tri(num, 1)]);
   }
@@ -193,12 +194,13 @@ class MTriangle6 : public MTriangle {
   virtual int getPolynomialOrder() const { return 2; }
   virtual int getNumVertices() const { return 6; }
   virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
+  virtual const MVertex *getVertex(int num) const { return num < 3 ? _v[num] : _vs[num - 3]; }  
   virtual MVertex *getVertexUNV(int num)
   {
     static const int map[6] = {0, 3, 1, 4, 2, 5};
     return getVertex(map[num]);
   }
-  virtual void xyz2uvw(double xyz[3], double uvw[3]){ MElement::xyz2uvw(xyz, uvw); }
+  virtual void xyz2uvw(double xyz[3], double uvw[3]) const{ MElement::xyz2uvw(xyz, uvw); }
   virtual int getNumEdgeVertices() const { return 3; }
   virtual int getNumEdgesRep();
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
@@ -274,6 +276,7 @@ class MTriangleN : public MTriangle {
   virtual int getPolynomialOrder() const { return _order; }
   virtual int getNumVertices() const { return 3 + _vs.size(); }
   virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
+  virtual const MVertex *getVertex(int num) const { return num < 3 ? _v[num] : _vs[num - 3]; }
   virtual int getNumFaceVertices() const
   {
     if(_order == 3 && _vs.size() == 6) return 0;
@@ -294,7 +297,7 @@ class MTriangleN : public MTriangle {
     if(_order == 10  && _vs.size() == 27) return 0;
     return 0;
   }
-  virtual void xyz2uvw(double xyz[3], double uvw[3]){ MElement::xyz2uvw(xyz, uvw); }
+  virtual void xyz2uvw(double xyz[3], double uvw[3]) const { MElement::xyz2uvw(xyz, uvw); }
   virtual int getNumEdgeVertices() const { return 3 * (_order - 1); }
   virtual int getNumEdgesRep();
   virtual int getNumFacesRep();
-- 
GitLab