diff --git a/Geo/MElement.h b/Geo/MElement.h
index 7aff9cc4dc4bff9cfa5d0cd8933194408deffde2..b3d1e46fa0a90700a4f427747ab55cf2af181ad4 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -11,16 +11,18 @@
 #include <string>
 
 #include "GmshMessage.h"
+#include "ElementType.h"
 #include "MVertex.h"
 #include "MEdge.h"
 #include "MFace.h"
+#include "FuncSpaceData.h"
 #include "nodalBasis.h"
 #include "polynomialBasis.h"
-#include "JacobianBasis.h"
 #include "MetricBasis.h"
 #include "GaussIntegration.h"
 
 class GModel;
+class JacobianBasis;
 
 // A mesh element.
 class MElement
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 862ca1acc499f1cc436909d4e2e4f3981eb97471..b1084ca229919c5a89af7e0c25999184eb90fae9 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -13,6 +13,7 @@
 #include "MHexahedron.h"
 #include "Numeric.h"
 #include "polynomialBasis.h"
+#include "JacobianBasis.h"
 #include "GmshMessage.h"
 #include <limits>
 #include <string.h>
diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index 97ed21a7e176404cf6e2c5ad254b5c9b52ed767e..604ab7f0642d23951b86d3ed5688997065d35022 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -7,9 +7,11 @@
 #include "GmshDefines.h"
 #include "polynomialBasis.h"
 #include "pyramidalBasis.h"
+#include "bezierBasis.h"
 #include "miniBasis.h"
 #include "MetricBasis.h"
 #include "CondNumBasis.h"
+#include "JacobianBasis.h"
 #include <map>
 #include <cstddef>
 
@@ -82,6 +84,25 @@ const JacobianBasis* BasisFactory::getJacobianBasis(FuncSpaceData fsd)
   return J;
 }
 
+const JacobianBasis* BasisFactory::getJacobianBasis(int tag, int order)
+{
+  const int type = ElementType::ParentTypeFromTag(tag);
+  if (type != TYPE_PYR)
+    return getJacobianBasis(FuncSpaceData(true, tag, order));
+  else
+    return getJacobianBasis(FuncSpaceData(true, tag, false, order+1, order));
+}
+
+const JacobianBasis* BasisFactory::getJacobianBasis(int tag)
+{
+  const int order = JacobianBasis::jacobianOrder(tag);
+  const int type = ElementType::ParentTypeFromTag(tag);
+  if (type != TYPE_PYR)
+    return getJacobianBasis(FuncSpaceData(true, tag, order));
+  else
+    return getJacobianBasis(FuncSpaceData(true, tag, false, order+2, order));
+}
+
 const MetricBasis* BasisFactory::getMetricBasis(int tag)
 {
   std::map<int, MetricBasis*>::const_iterator it = ms.find(tag);
@@ -112,6 +133,16 @@ const GradientBasis* BasisFactory::getGradientBasis(FuncSpaceData data)
   return G;
 }
 
+const GradientBasis* BasisFactory::getGradientBasis(int tag, int order)
+{
+  return getGradientBasis(FuncSpaceData(true, tag, order));
+}
+
+const GradientBasis* BasisFactory::getGradientBasis(int tag)
+{
+  return getGradientBasis(FuncSpaceData(tag));
+}
+
 const bezierBasis* BasisFactory::getBezierBasis(FuncSpaceData fsd)
 {
   FuncSpaceData data = fsd.getForPrimaryElement();
@@ -124,6 +155,17 @@ const bezierBasis* BasisFactory::getBezierBasis(FuncSpaceData fsd)
   return B;
 }
 
+const bezierBasis* BasisFactory::getBezierBasis(int parentTag, int order)
+{
+  int primaryTag = ElementType::getTag(parentTag, 1);
+  return getBezierBasis(FuncSpaceData(true, primaryTag, order));
+}
+
+const bezierBasis* BasisFactory::getBezierBasis(int tag)
+{
+  return getBezierBasis(FuncSpaceData(tag));
+}
+
 void BasisFactory::clearAll()
 {
   std::map<int, nodalBasis*>::iterator itF = fs.begin();
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index c75bcd2e4dba6b05eb2ed7d8dfc1e784cc913463..1ca2457dd2da7f341b1ce491e85f4e6d13540ca9 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -6,13 +6,14 @@
 #ifndef BASISFACTORY_H
 #define BASISFACTORY_H
 
-#include "JacobianBasis.h"
-#include "FuncSpaceData.h"
+#include <map>
 class nodalBasis;
 class MetricBasis;
 class GradientBasis;
 class bezierBasis;
 class CondNumBasis;
+class JacobianBasis;
+class FuncSpaceData;
 
 class BasisFactory
 {
@@ -34,21 +35,8 @@ class BasisFactory
   // Warning: bases returned by BasisFactory::getJacobianBasis(int tag) are the
   // only safe bases for using Bezier on the jacobian determinant!
   static const JacobianBasis* getJacobianBasis(FuncSpaceData);
-  static const JacobianBasis* getJacobianBasis(int tag, int order) {
-    const int type = ElementType::ParentTypeFromTag(tag);
-    if (type != TYPE_PYR)
-      return getJacobianBasis(FuncSpaceData(true, tag, order));
-    else
-      return getJacobianBasis(FuncSpaceData(true, tag, false, order+1, order));
-  }
-  static const JacobianBasis* getJacobianBasis(int tag) {
-    const int order = JacobianBasis::jacobianOrder(tag);
-    const int type = ElementType::ParentTypeFromTag(tag);
-    if (type != TYPE_PYR)
-      return getJacobianBasis(FuncSpaceData(true, tag, order));
-    else
-      return getJacobianBasis(FuncSpaceData(true, tag, false, order+2, order));
-  }
+  static const JacobianBasis* getJacobianBasis(int tag, int order);
+  static const JacobianBasis* getJacobianBasis(int tag);
 
   // Metric
   static const MetricBasis* getMetricBasis(int tag);
@@ -58,22 +46,13 @@ class BasisFactory
 
   // Gradients
   static const GradientBasis* getGradientBasis(FuncSpaceData);
-  static const GradientBasis* getGradientBasis(int tag, int order) {
-    return getGradientBasis(FuncSpaceData(true, tag, order));
-  }
-  static const GradientBasis* getGradientBasis(int tag) {
-    return getGradientBasis(FuncSpaceData(tag));
-  }
+  static const GradientBasis* getGradientBasis(int tag, int order);
+  static const GradientBasis* getGradientBasis(int tag);
 
   // Bezier
   static const bezierBasis* getBezierBasis(FuncSpaceData);
-  static const bezierBasis* getBezierBasis(int parentTag, int order) {
-    int primaryTag = ElementType::getTag(parentTag, 1);
-    return getBezierBasis(FuncSpaceData(true, primaryTag, order));
-  }
-  static const bezierBasis* getBezierBasis(int tag) {
-    return getBezierBasis(FuncSpaceData(tag));
-  }
+  static const bezierBasis* getBezierBasis(int parentTag, int order);
+  static const bezierBasis* getBezierBasis(int tag);
 
   static void clearAll();
 };
diff --git a/Numeric/CondNumBasis.cpp b/Numeric/CondNumBasis.cpp
index 2a514f962c9b63671528de6f003fa889d56ba660..eb4f72e94b0db6a644447c59d3147e15160de300 100644
--- a/Numeric/CondNumBasis.cpp
+++ b/Numeric/CondNumBasis.cpp
@@ -312,7 +312,7 @@ inline void calcGradInvCondNum3D(double dxdX, double dxdY, double dxdZ,
 
 
 CondNumBasis::CondNumBasis(int tag, int cnOrder) :
-    _tag(tag), _dim(ElementType::DimensionFromTag(tag)),
+    _dim(ElementType::DimensionFromTag(tag)),
     _condNumOrder(cnOrder >= 0 ? cnOrder : condNumOrder(tag))
 {
   if ( ElementType::ParentTypeFromTag(tag) == TYPE_TRIH){
diff --git a/Numeric/CondNumBasis.h b/Numeric/CondNumBasis.h
index 82a735996017370eaa132ec7d48b4f4285e0c475..ab598583152232be8e5da12437ff8c83a4b281be 100644
--- a/Numeric/CondNumBasis.h
+++ b/Numeric/CondNumBasis.h
@@ -8,15 +8,14 @@
 
 #include <map>
 #include <vector>
-#include "JacobianBasis.h"
 #include "fullMatrix.h"
-
+#include "JacobianBasis.h"
 
 class CondNumBasis {
  private:
   const GradientBasis *_gradBasis;
 
-  const int _tag, _dim, _condNumOrder;
+  const int _dim, _condNumOrder;
 
   fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY,
                      primGradShapeBarycenterZ;
diff --git a/Numeric/FuncSpaceData.h b/Numeric/FuncSpaceData.h
index fbc55f6d0f18a722c51ee6dc6329bf9a4755c623..546deb052c35a8898a941ee86b445077f44fd327 100644
--- a/Numeric/FuncSpaceData.h
+++ b/Numeric/FuncSpaceData.h
@@ -35,7 +35,6 @@ private:
   // otherwise,
   //   the space is {X^i Y^j Z^k | i,j <= '_nij', k <= '_nk'}, (hex-like space)
   // where X = xi/(1-zeta), Y = eta/(1-zeta) and Z = (1-zeta).
-  // Note that (x, y, z) are here the reference coordinates.
 
 public:
 
diff --git a/contrib/HighOrderMeshOptimizer/CADDistances.cpp b/contrib/HighOrderMeshOptimizer/CADDistances.cpp
index 242b0584a8d78dd3a7231042bcca2592ca7506d2..1969f5793c4d3ebea74c485a1b0fe23c6b35dd9e 100644
--- a/contrib/HighOrderMeshOptimizer/CADDistances.cpp
+++ b/contrib/HighOrderMeshOptimizer/CADDistances.cpp
@@ -10,6 +10,7 @@
 #include "MTriangle.h"
 #include "MQuadrangle.h"
 #include "BasisFactory.h"
+#include "JacobianBasis.h"
 #include "GModel.h"
 #if defined(HAVE_ANN)
 #include "ANN/ANN.h"
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 9ffdebb4e08322731d073e53568550d3db4700be..6c8534f6ebacaca6c784587f9dc97bddc402cc4c 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -35,6 +35,7 @@
 #include "ParamCoord.h"
 #include "OptHomMesh.h"
 #include "BasisFactory.h"
+#include "JacobianBasis.h"
 #include "OptHomIntegralBoundaryDist.h"
 
 Mesh::Mesh(const std::map<MElement*,GEntity*> &element2entity,