diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 428a11c2082209b1d0e44922b5787f54bcdc6f8e..a934ce6670a549efe6ab1abf06df558c4a9d435c 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -450,7 +450,7 @@ const {
   return _computeDeterminantAndRegularize(this, jac);
 }
 
-double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+double MElement::getPrimaryJacobian(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.;
@@ -471,7 +471,7 @@ double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][
   return _computeDeterminantAndRegularize(this, jac);
 }
 
-void MElement::getSignedJacobian(fullVector<double> &jacobian, int o)
+void MElement::getSignedJacobian(fullVector<double> &jacobian, int o) const
 {
   const int numNodes = getNumVertices();
   fullMatrix<double> nodesXYZ(numNodes,3);
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 480b468221d6b8a7a663f7bf98976f26ffb95c66..b4d59b9d9b78fe58bc0ecd77c5cfad3e7cce23fa 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -268,16 +268,16 @@ class MElement
     }
     return detJ;
   }
-  virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
-  double getJacobianDeterminant(double u, double v, double w)
+  virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]) const;
+  double getJacobianDeterminant(double u, double v, double w) const
   {
     double jac[3][3]; return getJacobian(u, v, w, jac);
   }
-  void getSignedJacobian(fullVector<double> &jacobian, int o = -1);
+  void getSignedJacobian(fullVector<double> &jacobian, int o = -1) const;
   void getNodesCoord(fullMatrix<double> &nodesXYZ) const;
-  virtual int getNumShapeFunctions() const{ return getNumVertices(); }
-  virtual int getNumPrimaryShapeFunctions() { return getNumPrimaryVertices(); }
-  virtual const MVertex *getShapeFunctionNode(int i) const{ return getVertex(i); }
+  virtual int getNumShapeFunctions() const { return getNumVertices(); }
+  virtual int getNumPrimaryShapeFunctions() const { 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