diff --git a/Numeric/approximationError.cpp b/Numeric/approximationError.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..407aaa5d5ce5a343abec8f8c1a65241862c96e28
--- /dev/null
+++ b/Numeric/approximationError.cpp
@@ -0,0 +1,26 @@
+#include "approximationError.h"
+#include "MElement.h"
+double approximationError (simpleFunction<double> &f, MElement *e) 
+{
+  double VALS [e->getNumVertices()];
+  for (int i=0;i<e->getNumVertices();i++){
+    MVertex *v = e->getVertex(i);
+    VALS[i] = f(v->x(),v->y(),v->z());
+  }
+  int npts; IntPt *pts;
+  e->getIntegrationPoints (2*e->getPolynomialOrder() + 2 , &npts, &pts);
+  double errSqr = 0.0;
+  for (int k=0;k<npts;k++){
+    const double u = pts[k].pt[0];
+    const double v = pts[k].pt[1];
+    const double w = pts[k].pt[2];
+    SPoint3 p; e->pnt(u, v, w, p);
+    const double Jac = e->getJacobianDeterminant(u,v,w);
+    const double C = e->interpolate(VALS,u,v,w);
+    const double F = f(p.x(),p.y(),p.z());
+    const double DIFF = C-F;
+    errSqr += pts[k].weight * Jac * (DIFF*DIFF);
+  }
+  return sqrt(errSqr);
+}
+
diff --git a/Numeric/approximationError.h b/Numeric/approximationError.h
new file mode 100644
index 0000000000000000000000000000000000000000..53caed44370384122721bb843bb7e9863319e222
--- /dev/null
+++ b/Numeric/approximationError.h
@@ -0,0 +1,10 @@
+#ifndef _APPROXIMATION_ERROR_
+#define _APPROXIMATION_ERROR_
+#include <map>
+#include "simpleFunction.h"
+class MElement;
+// computes E such as
+// E^2 = \int_e [C_e(f) - f]^2 de
+// where C_e(f) is clement's interpolation operator of f on e
+double approximationError (simpleFunction<double> &f, MElement *e) ;
+#endif