diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index a8229ea0219f4d4c73e2d4110818457d61836a1c..6581b74f101bae2aaca104f0138a8d485fb16072 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -145,7 +145,7 @@ class MPolyhedron : public MElement {
   {
     return (_orig ? _orig->getNumShapeFunctions() : 0);
   }
-  virtual int getNumPrimaryShapeFunctions()
+  virtual int getNumPrimaryShapeFunctions() const
   {
     return (_orig ? _orig->getNumPrimaryShapeFunctions() : 0);
   }
@@ -295,7 +295,7 @@ class MPolygon : public MElement {
   {
     return (_orig ? _orig->getNumShapeFunctions() : 0);
   }
-  virtual int getNumPrimaryShapeFunctions()
+  virtual int getNumPrimaryShapeFunctions() const
   {
     return (_orig ? _orig->getNumPrimaryShapeFunctions() : 0);
   }
diff --git a/Geo/MSubElement.cpp b/Geo/MSubElement.cpp
index ae58a393d8ecd4785ff69da95161a1a5ab411d53..2403623b9bf7f938f92e8460d904e970a47c5f54 100644
--- a/Geo/MSubElement.cpp
+++ b/Geo/MSubElement.cpp
@@ -64,7 +64,7 @@ double MSubTetrahedron::getJacobian(double u, double v, double w, double jac[3][
   if(_orig) return _orig->getJacobian(u, v, w, jac);
   return 0;
 }
-double MSubTetrahedron::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+double MSubTetrahedron::getPrimaryJacobian(double u, double v, double w, double jac[3][3]) const
 {
   if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
   return 0;
@@ -76,7 +76,7 @@ int MSubTetrahedron::getNumShapeFunctions() const
   return 0;
 }
 
-int MSubTetrahedron::getNumPrimaryShapeFunctions()
+int MSubTetrahedron::getNumPrimaryShapeFunctions() const
 {
   if(_orig) return _orig->getNumPrimaryShapeFunctions();
   return 0;
@@ -303,7 +303,7 @@ double MSubTriangle::getJacobian(double u, double v, double w, double jac[3][3])
   if(_orig) return _orig->getJacobian(u, v, w, jac);
   return 0;
 }
-double MSubTriangle::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+double MSubTriangle::getPrimaryJacobian(double u, double v, double w, double jac[3][3]) const
 {
   if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
   return 0;
@@ -315,7 +315,7 @@ int MSubTriangle::getNumShapeFunctions() const
   return 0;
 }
 
-int MSubTriangle::getNumPrimaryShapeFunctions()
+int MSubTriangle::getNumPrimaryShapeFunctions() const
 {
   if(_orig) return _orig->getNumPrimaryShapeFunctions();
   return 0;
@@ -531,7 +531,7 @@ double MSubLine::getJacobian(double u, double v, double w, double jac[3][3]) con
   if(_orig) return _orig->getJacobian(u, v, w, jac);
   return 0;
 }
-double MSubLine::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+double MSubLine::getPrimaryJacobian(double u, double v, double w, double jac[3][3]) const
 {
   if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
   return 0;
@@ -543,7 +543,7 @@ int MSubLine::getNumShapeFunctions() const
   return 0;
 }
 
-int MSubLine::getNumPrimaryShapeFunctions()
+int MSubLine::getNumPrimaryShapeFunctions() const
 {
   if(_orig) return _orig->getNumPrimaryShapeFunctions();
   return 0;
@@ -723,7 +723,7 @@ double MSubPoint::getJacobian(double u, double v, double w, double jac[3][3]) co
   if(_orig) return _orig->getJacobian(u, v, w, jac);
   return 0;
 }
-double MSubPoint::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+double MSubPoint::getPrimaryJacobian(double u, double v, double w, double jac[3][3]) const
 {
   if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
   return 0;
@@ -735,7 +735,7 @@ int MSubPoint::getNumShapeFunctions() const
   return 0;
 }
 
-int MSubPoint::getNumPrimaryShapeFunctions()
+int MSubPoint::getNumPrimaryShapeFunctions() const
 {
   if(_orig) return _orig->getNumPrimaryShapeFunctions();
   return 0;
diff --git a/Geo/MSubElement.h b/Geo/MSubElement.h
index 9415ead679716660ab1d0bf9db955717a69e9034..dc6da9b3c2c8495e356fa4ce4cce83ac3cb21681 100644
--- a/Geo/MSubElement.h
+++ b/Geo/MSubElement.h
@@ -51,9 +51,9 @@ class MSubTetrahedron : public MTetrahedron
   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 double getPrimaryJacobian(double u, double v, double w, double jac[3][3]) const;
   virtual int getNumShapeFunctions() const;
-  virtual int getNumPrimaryShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions() const;
   virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
@@ -109,9 +109,9 @@ class MSubTriangle : public MTriangle
   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 double getPrimaryJacobian(double u, double v, double w, double jac[3][3]) const;
   virtual int getNumShapeFunctions() const;
-  virtual int getNumPrimaryShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions() const;
   virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
@@ -167,9 +167,9 @@ class MSubLine : public MLine
   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 double getPrimaryJacobian(double u, double v, double w, double jac[3][3]) const;
   virtual int getNumShapeFunctions() const;
-  virtual int getNumPrimaryShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions() const;
   virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
@@ -225,9 +225,9 @@ class MSubPoint : public MPoint
   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 double getPrimaryJacobian(double u, double v, double w, double jac[3][3]) const;
   virtual int getNumShapeFunctions() const;
-  virtual int getNumPrimaryShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions() const;
   virtual const MVertex* getShapeFunctionNode(int i) const;
   virtual MVertex* getShapeFunctionNode(int i);
   virtual void xyz2uvw(double xyz[3], double uvw[3]) const;