diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index a39c5a0ac4fcc505fac4e8f898b3a598fcdfd35e..2d9d4949c378b16bcae33bdf8439a3c7c2049cfb 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -1,6 +1,6 @@
 #ifndef AMDIS_H
 #define AMDIS_H
-#include "MathFunctions.h"
+
 #include "AbstractFunction.h"
 #include "AdaptInfo.h"
 #include "AdaptInstationary.h"
@@ -49,6 +49,7 @@
 #include "MacroElement.h"
 #include "MacroWriter.h"
 #include "Marker.h"
+#include "MathFunctions.h"
 #include "MatrixVector.h"
 #include "Mesh.h"
 #include "MeshStructure.h"
diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h
index 272c89046e90c3b7dfb47f8e880d17e1a9db02a4..686225f7d0b94aa305c50ee4a1f522c3a6983ebc 100644
--- a/AMDiS/src/AMDiS_fwd.h
+++ b/AMDiS/src/AMDiS_fwd.h
@@ -101,6 +101,19 @@ namespace AMDiS {
   struct AtomicBoundary;
 
   template<typename ReturnType, typename ArgumentType> class AbstractFunction;
+  template<typename ReturnType, 
+	   typename ArgumentType1, 
+	   typename ArgumentType2>                     class BinaryAbstractFunction;
+  template<typename ReturnType, 
+	   typename ArgumentType1, 
+	   typename ArgumentType2,
+	   typename ArgumentType3>                     class TertiaryAbstractFunction;
+  template<typename ReturnType, 
+	   typename ArgumentType1, 
+	   typename ArgumentType2, 
+	   typename ArgumentType3,
+	   typename ArgumentType4>                     class QuartAbstractFunction;
+
   template<typename T>                                 class DOFIndexed;
   template<typename T>                                 class DOFVectorBase;
   template<typename T>                                 class DOFVector;
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index cb978648331bbb4e5ce4ef25a05d365ba584be72..ab89835fc07bb09193c3c67459e373cd91269a8f 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -865,5 +865,97 @@ namespace AMDiS {
   }
 
 
+  double integrate(const DOFVector<double> &vec1,
+		   const DOFVector<double> &vec2,
+		   BinaryAbstractFunction<double, double, double> *fct)
+  {
+    if (vec1.getFESpace()->getMesh() == vec2.getFESpace()->getMesh())
+      return intSingleMesh(vec1, vec2, fct);
+    else
+      return intDualMesh(vec1, vec2, fct);
+  }
+
+
+  double intSingleMesh(const DOFVector<double> &vec1,
+		       const DOFVector<double> &vec2,
+		       BinaryAbstractFunction<double, double, double> *fct)
+  {
+    FUNCNAME("intDualmesh()");
+    
+    TEST_EXIT(fct)("No function defined!\n");
+
+    int deg = 2 * vec1.getFESpace()->getBasisFcts()->getDegree();
+    Quadrature* quad = 
+      Quadrature::provideQuadrature(vec1.getFESpace()->getMesh()->getDim(), deg);
+    FastQuadrature *fastQuad =
+      FastQuadrature::provideFastQuadrature(vec1.getFESpace()->getBasisFcts(), *quad, INIT_PHI);
+
+    std::vector<double> qp1(fastQuad->getNumPoints());
+    std::vector<double> qp2(fastQuad->getNumPoints());
+
+    double value = 0.0;
+    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
+    TraverseStack stack;    
+    ElInfo *elInfo = stack.traverseFirst(vec1.getFESpace()->getMesh(), -1, traverseFlag);
+    while (elInfo) {
+      double det = elInfo->getDet();
+      double tmp = 0.0;      
+      vec1.getVecAtQPs(elInfo, quad, fastQuad, &(qp1[0]));
+      vec2.getVecAtQPs(elInfo, quad, fastQuad, &(qp2[0]));      
+      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
+ 	tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
+      
+      value += tmp * det;
+      
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    return value;
+  }
+
+
+  double intDualMesh(const DOFVector<double> &vec1,
+		     const DOFVector<double> &vec2,
+		     BinaryAbstractFunction<double, double, double> *fct)
+  {
+    FUNCNAME("intDualmesh()");
+    
+    TEST_EXIT(fct)("No function defined!\n");
+
+    int deg = 2 * vec1.getFESpace()->getBasisFcts()->getDegree();
+    Quadrature* quad = 
+      Quadrature::provideQuadrature(vec1.getFESpace()->getMesh()->getDim(), deg);
+    FastQuadrature *fastQuad =
+      FastQuadrature::provideFastQuadrature(vec1.getFESpace()->getBasisFcts(), *quad, INIT_PHI);
+
+    std::vector<double> qp1(fastQuad->getNumPoints());
+    std::vector<double> qp2(fastQuad->getNumPoints());
+
+    double value = 0.0;
+    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
+    DualTraverse dualTraverse;    
+    DualElInfo dualElInfo;
+    bool cont = dualTraverse.traverseFirst(vec1.getFESpace()->getMesh(),
+					   vec2.getFESpace()->getMesh(),
+					   -1, -1, traverseFlag, traverseFlag,
+					   dualElInfo);
+    while (cont) {
+      double det = dualElInfo.smallElInfo->getDet();
+      double tmp = 0.0;      
+      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo,
+		       quad, fastQuad, &(qp1[0]));
+      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo,
+		       quad, fastQuad, &(qp2[0]));      
+      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
+ 	tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
+      
+      value += tmp * det;
+      
+      cont = dualTraverse.traverseNext(dualElInfo);
+    }
+
+    return value;
+  }
+
 }
 
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 4ba37da1adfe1dce7d2e8ebdb196ebf3701964f3..d8716809661a8bd50f1ca4f2a37dd4da74c58cd0 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -789,6 +789,26 @@ namespace AMDiS {
 
   WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
 					     WorldVector<DOFVector<double>*> *result);
+
+  /** \brief
+   * Computes the integral of a function that includes two different DOFVectors. This
+   * function works also for the case that the DOFVectors are defined on two different
+   * meshes.
+   */
+  double integrate(const DOFVector<double> &vec1,
+		   const DOFVector<double> &vec2,
+		   BinaryAbstractFunction<double, double, double> *fct);
+
+  /// Computes the integral of a function with two DOFVectors defined on the same mesh.
+  double intSingleMesh(const DOFVector<double> &vec1,
+		       const DOFVector<double> &vec2,
+		       BinaryAbstractFunction<double, double, double> *fct);
+
+  /// Computes the integral of a function with two DOFVectors defined on different meshes.
+  double intDualMesh(const DOFVector<double> &vec1,
+		     const DOFVector<double> &vec2,
+		     BinaryAbstractFunction<double, double, double> *fct);
+
 }
 
 #include "DOFVector.hh"