diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index bb255e2cfe35b8133328150bab438d4ce95effc1..d1358583d44b496107409b7d720be998dffad72c 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -125,7 +125,6 @@ namespace AMDiS {
     for (int i = 0; i < mSize; i++) {
       matrix[i].resize(0);
     }
-    return;
   }
 
   bool DOFMatrix::symmetric()
@@ -772,4 +771,47 @@ namespace AMDiS {
     return(result);
   }
 
+  void addDOFMatrix(DOFMatrix *result, const DOFMatrix *a, const DOFMatrix *b)
+  {
+    result->clear();
+
+    DOFMatrix::Iterator resultIterator(result, USED_DOFS);
+    DOFMatrix::Iterator aIterator(const_cast<DOFMatrix*>(a), USED_DOFS);
+    DOFMatrix::Iterator bIterator(const_cast<DOFMatrix*>(b), USED_DOFS);
+
+    for (resultIterator.reset(), aIterator.reset(); 
+	 !aIterator.end(); 
+	 ++resultIterator, ++aIterator) {
+
+      *resultIterator = *aIterator;
+    }      
+
+    for (resultIterator.reset(), bIterator.reset(); 
+	 !bIterator.end(); 
+	 ++resultIterator, ++bIterator) {
+      std::vector<MatEntry>::iterator resultRowIt;
+      std::vector<MatEntry>::const_iterator bRowIt;
+
+      for (bRowIt = bIterator->begin();
+	   bRowIt != bIterator->end();
+	   ++bRowIt) {
+	bool added = false;
+	for (resultRowIt = resultIterator->begin();
+	     resultRowIt != resultIterator->end();
+	     ++resultRowIt) {
+	  if (bRowIt->col == resultRowIt->col) {
+	    resultRowIt->entry += bRowIt->entry;
+	    
+	    added = true;
+	    break;
+	  }
+	}
+
+	if (!added) {
+	  resultIterator->push_back(*bRowIt);
+	}
+      }	         
+    }      
+  }
+
 }
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index a83fe2dd9eb5489abb728f25113564b701accf32..77b51fe238c4c2b6939d296e912457d6e9c7d32e 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -739,6 +739,8 @@ namespace AMDiS {
 
   double max(std::vector<MatEntry> *row);
 
+  void addDOFMatrix(DOFMatrix *result, const DOFMatrix *a, const DOFMatrix *b);
+
 }
 
 #endif  // AMDIS_DOFMATRIX_H