diff --git a/Mesh/DivideAndConquer.h b/Mesh/DivideAndConquer.h
index 1907f2dcf18e3a96f462a7eb8c0261e362a9122b..30fc2ecb3e53360f4af51e3cdedb969e0e3a8dfb 100644
--- a/Mesh/DivideAndConquer.h
+++ b/Mesh/DivideAndConquer.h
@@ -14,6 +14,7 @@
 #include "SPoint2.h"
 #include "simpleFunction.h"
 #include "Octree.h"
+#include "MElement.h"
 
 class GFace;
 typedef struct _CDLIST DListRecord, *DListPeek;
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index c531b3d8aeeea05232d222e9c67031bb0685e6ac..097588743a4e03129ee3cd92913c1658aa06d904 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -13,6 +13,7 @@
 #endif
 
 #include "MVertex.h"
+#include "MElement.h"
 //#include "directions3D.h"
 #include "BackgroundMesh.h"
 #include "intersectCurveSurface.h"
diff --git a/Numeric/simpleFunction.h b/Numeric/simpleFunction.h
index 09fce8f4856e97c3d596ebe1e9720edf677de07e..4b56eb33c22c19d1023a0b30bfa9e2db8e20a11d 100644
--- a/Numeric/simpleFunction.h
+++ b/Numeric/simpleFunction.h
@@ -7,7 +7,7 @@
 #define _SIMPLE_FUNCTION_H_
 
 // FIXME: Numeric/ should not depend on Geo/
-#include "MElement.h"
+class MElement;
 
 template <class scalar>
 class simpleFunction {
diff --git a/Numeric/simpleFunctionPython.h b/Numeric/simpleFunctionPython.h
new file mode 100644
index 0000000000000000000000000000000000000000..f3f0549e6c1e09a9dce40a21993d2c66e41d6234
--- /dev/null
+++ b/Numeric/simpleFunctionPython.h
@@ -0,0 +1,37 @@
+#ifndef _SIMPLE_FUNCTION_PYTHON_H_
+#define _SIMPLE_FUNCTION_PYTHON_H_
+#include "Python.h"
+#include "simpleFunction.h"
+
+class simpleFunctionPython : public simpleFunction<double> {
+  PyObject *_pycallback;
+ public:
+  simpleFunctionPython(PyObject *callback):
+    _pycallback(callback)
+  {
+    Py_INCREF(_pycallback);
+  }
+  ~simpleFunctionPython()
+  {
+    Py_DECREF(_pycallback);
+  }
+  double operator()(double x, double y, double z) const
+  {
+    PyObject *pyargs = Py_BuildValue("(ddd)", x, y, z);
+    PyObject *result = PyEval_CallObject(_pycallback, pyargs);
+    double r = 0;
+    if (result) {
+      int ok = PyArg_Parse(result, "d", &r);
+      if (not ok)
+        Msg::Error("The python function did not return a double.");
+      Py_DECREF(result);
+    }
+    else {
+      PyErr_Print();
+      Msg::Error("An error occurs in the python simple function.");
+    }
+    Py_DECREF(pyargs);
+    return r;
+  }
+};
+#endif
diff --git a/contrib/bamg/bamglib/Mesh2.h b/contrib/bamg/bamglib/Mesh2.h
index 4d02159a19ecf2d45d12187d58696d6d5d464214..ce475c00063aa043b4bbc474c11fe77fe3309b34 100644
--- a/contrib/bamg/bamglib/Mesh2.h
+++ b/contrib/bamg/bamglib/Mesh2.h
@@ -30,6 +30,7 @@
 #include <math.h>
 #include <limits.h>
 #include <time.h>
+#include "MElement.h"
 #if  (defined(unix) || defined(__unix)) && !defined(__AIX)
 #define SYSTIMES
 #include <sys/times.h>
diff --git a/wrappers/gmshpy/gmshNumeric.i b/wrappers/gmshpy/gmshNumeric.i
index 11a968608f61360c2e7f75919066ee707080c886..b57c4182a70ba71a1cf54bd9276c221591b86606 100644
--- a/wrappers/gmshpy/gmshNumeric.i
+++ b/wrappers/gmshpy/gmshNumeric.i
@@ -11,6 +11,8 @@
   #include "JacobianBasis.h"
   #include "fullMatrix.h"
   #include "nodalBasis.h"
+  #include "simpleFunction.h"
+  #include "simpleFunctionPython.h"
   #include "polynomialBasis.h"
   #include "pyramidalBasis.h"
 %}
@@ -18,8 +20,11 @@
 %include "GaussIntegration.h"
 %include "JacobianBasis.h"
 %include "fullMatrix.h"
+%include "simpleFunction.h"
 %template(fullMatrixDouble) fullMatrix<double>;
 %template(fullVectorDouble) fullVector<double>;
+%template(simpleFunctionDouble) simpleFunction<double>;
+%include "simpleFunctionPython.h"
 %include "nodalBasis.h"
 %include "polynomialBasis.h"
 %include "pyramidalBasis.h"