diff --git a/AMDiS/compositeFEM/src/CompositeFEMOperator.cc b/AMDiS/compositeFEM/src/CompositeFEMOperator.cc
index 2806a5fcd3707b08c89dac0d85bb26f191eded55..8206a4592271442a452c5d735c109616d737e7b9 100644
--- a/AMDiS/compositeFEM/src/CompositeFEMOperator.cc
+++ b/AMDiS/compositeFEM/src/CompositeFEMOperator.cc
@@ -1,16 +1,13 @@
+#include <boost/numeric/mtl/mtl.hpp>
 #include "CompositeFEMOperator.h"
-
-#include "ElementMatrix.h"
-#include "ElementVector.h"
 #include "OpenMP.h"
-
 #include "SubElementAssembler.h"
 #include "SubElInfo.h"
 #include "SubPolytope.h"
 
 void 
 CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo, 
-				       ElementMatrix *userMat, 
+				       ElementMatrix& userMat, 
 				       double factor)
 {
   FUNCNAME("CompositeFEMOperator::getElementMatrix");
@@ -18,9 +15,6 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
   VectorOfFixVecs<DimVec<double> > *intersecPoints = NULL;
   SubPolytope *subPolytope = NULL;
   double levelSetSubPolytope;
-  ElementMatrix *elMat;
-  ElementMatrix *subPolMat1;
-  ElementMatrix *subPolMat2;
   DimVec<double> subElVertexBarCoords(elInfo->getMesh()->getDim());
 
   /**
@@ -99,9 +93,9 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
     ERROR_EXIT("cannot get position of subpolytope\n");
   }
 
-  subPolMat1 = NEW ElementMatrix(subElementAssembler->getNRow(),
-				 subElementAssembler->getNCol());
-  subPolMat1->set(0.0);
+  ElementMatrix subPolMat1(subElementAssembler->getNRow(),
+			   subElementAssembler->getNCol());
+  set_to_zero(subPolMat1);
   subElementAssembler->getSubPolytopeMatrix(subPolytope,
 					    subElementAssembler,
 					    elInfo,
@@ -110,12 +104,12 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
   /**
    * Integration on second subpolytope produced by the intersection.
    */
-  elMat = NEW ElementMatrix(subElementAssembler->getNRow(),
-			    subElementAssembler->getNCol());
-  elMat->set(0.0);
-  subPolMat2 = NEW ElementMatrix(subElementAssembler->getNRow(),
-				 subElementAssembler->getNCol());
-  subPolMat2->set(0.0);
+  ElementMatrix elMat(subElementAssembler->getNRow(),
+		      subElementAssembler->getNCol());
+  set_to_zero(elMat);
+  ElementMatrix subPolMat2(subElementAssembler->getNRow(),
+			   subElementAssembler->getNCol());
+  set_to_zero(subPolMat2);
 
   int myRank = omp_get_thread_num();
 
@@ -138,28 +132,21 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
 					    elInfo,
 					    subPolMat2);
 
-  axpy(-1.0, *subPolMat2, *elMat);
+  elMat -= subPolMat2;
 
   // Get integral on element as sum of the two integrals on subpolytopes.
-  axpy(1.0, *subPolMat1, *elMat);
+  elMat += subPolMat1;
 
   // Add integral to userMat.
-  axpy(factor, *elMat, *userMat);
-
-  /**
-   * Free data.
-   */
-  DELETE subPolytope;
-  DELETE elMat;
-  DELETE subPolMat1;
-  DELETE subPolMat2;
+  userMat += factor * elMat;
 
-  return;
+  // Free data
+  delete subPolytope;
 }
 
 void 
 CompositeFEMOperator::getElementVector(const ElInfo *elInfo, 
-				       ElementVector *userVec, 
+				       ElementVector& userVec, 
 				       double factor)
 {
   FUNCNAME("CompositeFEMOperator::getElementVector");
@@ -167,9 +154,6 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
   VectorOfFixVecs<DimVec<double> >*intersecPoints = NULL;
   SubPolytope *subPolytope = NULL;
   double levelSetSubPolytope;
-  ElementVector *elVec;
-  ElementVector *subPolVec1;
-  ElementVector *subPolVec2;
   DimVec<double> subElVertexBarCoords(elInfo->getMesh()->getDim());
 
   /**
@@ -247,8 +231,8 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
     ERROR_EXIT("cannot get position of subpolytope\n");
   }
 
-  subPolVec1 = NEW ElementVector(subElementAssembler->getNRow());
-  subPolVec1->set(0.0);
+  ElementVector subPolVec1(subElementAssembler->getNRow());
+  set_to_zero(subPolVec1);
   subElementAssembler->getSubPolytopeVector(subPolytope,
 					    subElementAssembler,
 					    elInfo,
@@ -257,10 +241,10 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
   /**
    * Integration on second subpolytope produced by the intersection.
    */
-  elVec = NEW ElementVector(subElementAssembler->getNRow());
-  elVec->set(0.0);
-  subPolVec2 = NEW ElementVector(subElementAssembler->getNRow());
-  subPolVec2->set(0.0);
+  ElementVector elVec(subElementAssembler->getNRow());
+  set_to_zero(elVec);
+  ElementVector subPolVec2(subElementAssembler->getNRow());
+  set_to_zero(subPolVec2);
 
   int myRank = omp_get_thread_num();
 
@@ -282,21 +266,14 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
 					    elInfo,
 					    subPolVec2);
 
-  axpy(-1.0, *subPolVec2, *elVec);
+  elVec -= subPolVec2;
 
   // Get integral on element as sum of the two integrals on subpolytopes.
-  axpy(1.0, *subPolVec1, *elVec);
+  elVec += subPolVec1;
 
   // Add integral to userVec.
-  axpy(factor, *elVec, *userVec);
-
-  /**
-   * Free data.
-   */
-  DELETE subPolytope;
-  DELETE elVec;
-  DELETE subPolVec1;
-  DELETE subPolVec2;
+  userVec += factor * elVec;
 
-  return;
+  // Free data
+  delete subPolytope;
 }
diff --git a/AMDiS/compositeFEM/src/CompositeFEMOperator.h b/AMDiS/compositeFEM/src/CompositeFEMOperator.h
index 7ef3c7a6bfa3d037a2102bdf95ffc0d708ba2d92..53521719144aae1b198d00409507ebeb822e7b02 100644
--- a/AMDiS/compositeFEM/src/CompositeFEMOperator.h
+++ b/AMDiS/compositeFEM/src/CompositeFEMOperator.h
@@ -6,7 +6,6 @@
 #include "FixVec.h"
 #include "Flag.h"
 #include "ElementLevelSet.h"
-#include "MemoryManager.h"
 #include "Operator.h"
 
 #include "ElementLevelSet.h"
@@ -35,10 +34,7 @@ using namespace std;
 class CompositeFEMOperator : public Operator
 {
 public:
-  MEMORY_MANAGED(CompositeFEMOperator);
-  /**
-   * Constructor.
-   */
+  /// Constructor.
   CompositeFEMOperator(Flag operatorType_,
 		       ElementLevelSet *elLS_,
 		       const FiniteElemSpace *rowFESpace_,
@@ -67,7 +63,7 @@ public:
    * integration domain.
    */
   void getElementMatrix(const ElInfo *elInfo, 
-			ElementMatrix *userMat, 
+			ElementMatrix& userMat, 
 			double factor = 1.0);
 
   /** 
@@ -78,7 +74,7 @@ public:
    * integration domain.
    */
   void getElementVector(const ElInfo *elInfo, 
-			ElementVector  *userVec, 
+			ElementVector& userVec, 
 			double factor = 1.0);
 
 protected:
diff --git a/AMDiS/compositeFEM/src/PenaltyOperator.cc b/AMDiS/compositeFEM/src/PenaltyOperator.cc
index 59daa976f9785eaab44d4e71d9ea580a3627b990..26ff5bdd6a820515c5214ea6141717db96caa5af 100644
--- a/AMDiS/compositeFEM/src/PenaltyOperator.cc
+++ b/AMDiS/compositeFEM/src/PenaltyOperator.cc
@@ -1,5 +1,4 @@
 #include "PenaltyOperator.h"
-
 #include "SurfaceOperator.h"
 
 double 
@@ -17,7 +16,7 @@ PenaltyOperator::getPenaltyCoeff(const ElInfo *elInfo)
 
 void 
 PenaltyOperator::getElementMatrix(const ElInfo *elInfo, 
-				  ElementMatrix *userMat, 
+				  ElementMatrix& userMat, 
 				  double factor)
 {
   FUNCNAME("PenaltyOperator::getElementMatrix");
@@ -135,7 +134,7 @@ PenaltyOperator::getElementMatrix(const ElInfo *elInfo,
 
 void 
 PenaltyOperator::getElementVector(const ElInfo *elInfo, 
-				  ElementVector *userVec, 
+				  ElementVector& userVec, 
 				  double factor)
 {
   FUNCNAME("PenaltyOperator::getElementVector");
@@ -189,12 +188,12 @@ PenaltyOperator::getElementVector(const ElInfo *elInfo,
       else {
 	surfaceOp->adaptSurfaceOperator((*tempCoords));
       }
-      surfaceOp->getElementVector(elInfo, userVec, factor*penaltyCoeff);
+      surfaceOp->getElementVector(elInfo, userVec, factor * penaltyCoeff);
       
       // Treat second simplex.
       (*tempCoords)[0] = (*intersecPoints)[3];
       surfaceOp->adaptSurfaceOperator((*tempCoords));
-      surfaceOp->getElementVector(elInfo, userVec, factor*penaltyCoeff);
+      surfaceOp->getElementVector(elInfo, userVec, factor * penaltyCoeff);
     }
 
     else {
diff --git a/AMDiS/compositeFEM/src/PenaltyOperator.h b/AMDiS/compositeFEM/src/PenaltyOperator.h
index a9de2fe556f9e2dadf459c9cfe653e31eacb9a7c..129af7722c81e63c39bc96c60a48707addda084c 100644
--- a/AMDiS/compositeFEM/src/PenaltyOperator.h
+++ b/AMDiS/compositeFEM/src/PenaltyOperator.h
@@ -5,10 +5,8 @@
 
 #include "ElementLevelSet.h"
 #include "Flag.h"
-#include "MemoryManager.h"
 #include "Operator.h"
 #include "SurfaceOperator.h"
-
 #include "ElementLevelSet.h"
 
 namespace AMDiS {
@@ -29,7 +27,6 @@ using namespace std;
 class PenaltyOperator : public Operator
 {
 public:
-  MEMORY_MANAGED(PenaltyOperator);
   /**
    * Constructor.
    */
@@ -76,7 +73,7 @@ public:
    * domain.
    */
   void getElementMatrix(const ElInfo *elInfo, 
-			ElementMatrix *userMat, 
+			ElementMatrix& userMat, 
 			double factor = 1.0);
 
   /** 
@@ -87,7 +84,7 @@ public:
    * the integration domain.
    */
   void getElementVector(const ElInfo *elInfo, 
-			ElementVector *userVec, 
+			ElementVector& userVec, 
 			double factor = 1.0);
 
 protected:
diff --git a/AMDiS/compositeFEM/src/SubElementAssembler.cc b/AMDiS/compositeFEM/src/SubElementAssembler.cc
index 95f495f5212cb385150f162adf2afd7959385272..569312d0e1604ee277018d72879e2cb0e6d9ec96 100644
--- a/AMDiS/compositeFEM/src/SubElementAssembler.cc
+++ b/AMDiS/compositeFEM/src/SubElementAssembler.cc
@@ -1,6 +1,5 @@
 #include <vector>
-#include "ElementMatrix.h"
-#include "ElementVector.h"
+#include <boost/numeric/mtl/mtl.hpp>
 #include "SubElementAssembler.h"
 #include "ScalableQuadrature.h"
 #include "SubPolytope.h"
@@ -72,7 +71,7 @@ namespace AMDiS {
 
   void SubElementAssembler::getSubElementVector(SubElInfo *subElInfo, 
 						const ElInfo *elInfo, 
-						ElementVector *userVec)
+						ElementVector& userVec)
   {
     /** 
      * Manipulate the quadratures of the SubAssemblers for subelement.
@@ -87,14 +86,13 @@ namespace AMDiS {
      * the result must be corrected with respect to subelement.
      */
     double corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
-    for (int i = 0; i < nRow; i++) {
-      (*userVec)[i] *= corrFactor;
-    }
+    for (int i = 0; i < nRow; i++)
+      userVec[i] *= corrFactor;
   }
 
   void SubElementAssembler::getSubElementMatrix(SubElInfo *subElInfo, 
 						const ElInfo *elInfo, 
-						ElementMatrix *userMat)
+						ElementMatrix& userMat)
   {
     /** 
      * Manipulate the quadratures of the SubAssemblers for subelement.
@@ -114,7 +112,7 @@ namespace AMDiS {
     double corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
     for (int i = 0; i < nRow; i++) {
       for (int j = 0; j < nCol; j++) {
-	(*userMat)[i][j] *= corrFactor;
+	userMat[i][j] *= corrFactor;
       }
     }
   }
@@ -122,44 +120,34 @@ namespace AMDiS {
   void SubElementAssembler::getSubPolytopeVector(SubPolytope *subPolytope,
 						 SubElementAssembler *subElementAssembler,
 						 const ElInfo *elInfo,
-						 ElementVector *subPolVec)
+						 ElementVector& subPolVec)
   {
-    /**
-     * Note: There is no reset of subPolVec.
-     */
+    /// Note: There is no reset of subPolVec.
     std::vector<SubElInfo *>::iterator it;
-    ElementVector *subElVec = NEW ElementVector(nRow);
+    ElementVector subElVec(nRow);
 
-    /**
-     * Assemble for each subelement of subpolytope.
-     */
+    /// Assemble for each subelement of subpolytope.
     for (it = subPolytope->getSubElementsBegin(); 
 	 it != subPolytope->getSubElementsEnd();
 	 it++) {
-      subElVec->set(0.0);
+      set_to_zero(subElVec);
       subElementAssembler->getSubElementVector(*it, elInfo, subElVec);
 
-      /**
-       * Add results for subelement to total result for subpolytope.
-       */
-      for (int i = 0; i < nRow; i++) {
-	(*subPolVec)[i] += (*subElVec)[i];
-      }
+      /// Add results for subelement to total result for subpolytope.
+      subPolVec += subElVec;
     }
-
-    DELETE subElVec;
   }
 
   void SubElementAssembler::getSubPolytopeMatrix(SubPolytope *subPolytope,
 						 SubElementAssembler *subElementAssembler,
 						 const ElInfo *elInfo,
-						 ElementMatrix *subPolMat)
+						 ElementMatrix& subPolMat)
   {
     /**
      * Note: There is no reset of subPolMat.
      */
     std::vector<SubElInfo *>::iterator it;
-    ElementMatrix *subElMat = NEW ElementMatrix(nRow, nCol);
+    ElementMatrix subElMat(nRow, nCol);
 
     /**
      * Assemble for each subelement of subpolytope.
@@ -167,20 +155,16 @@ namespace AMDiS {
     for (it = subPolytope->getSubElementsBegin(); 
 	 it != subPolytope->getSubElementsEnd();
 	 it++) {
-      subElMat->set(0.0);
+      set_to_zero(subElMat);
       subElementAssembler->getSubElementMatrix(*it, elInfo, subElMat);
 
       /**
        * Add results for subelement to total result for subpolytope.
        */
-      for (int i = 0; i < nRow; i++) {
-	for (int j = 0; j < nCol; j++) {
-	  (*subPolMat)[i][j] += (*subElMat)[i][j];
-	}
-      }
+      for (int i = 0; i < nRow; i++)
+	for (int j = 0; j < nCol; j++)
+	  subPolMat[i][j] += subElMat[i][j];
     }
-
-    DELETE subElMat;
   }
 
 }
diff --git a/AMDiS/compositeFEM/src/SubElementAssembler.h b/AMDiS/compositeFEM/src/SubElementAssembler.h
index 991b4b1875ded41f4d4e0c8d355e01709a55f207..a6c016c1ef4aaec691840830896e62038eeb5825 100644
--- a/AMDiS/compositeFEM/src/SubElementAssembler.h
+++ b/AMDiS/compositeFEM/src/SubElementAssembler.h
@@ -3,7 +3,6 @@
 #ifndef AMDIS_SUBELEMENTASSEMBLER_H
 #define AMDIS_SUBELEMENTASSEMBLER_H
 
-#include "MemoryManager.h"
 #include "Assembler.h"
 #include "SubElInfo.h"
 #include "ScalableQuadrature.h"
@@ -54,8 +53,6 @@ namespace AMDiS {
   class SubElementAssembler : public StandardAssembler
   {
   public:
-    MEMORY_MANAGED(SubElementAssembler);
-
     SubElementAssembler(Operator *op, 
 			const FiniteElemSpace *rowFESpace,
 			const FiniteElemSpace *colFESpace = NULL);
@@ -63,34 +60,34 @@ namespace AMDiS {
     virtual ~SubElementAssembler()
     {
       if (zeroOrderScalableQuadrature)
-	DELETE zeroOrderScalableQuadrature;
+	delete zeroOrderScalableQuadrature;
       if (firstOrderGrdPsiScalableQuadrature)
-	DELETE firstOrderGrdPsiScalableQuadrature;
+	delete firstOrderGrdPsiScalableQuadrature;
       if (firstOrderGrdPhiScalableQuadrature)
-	DELETE firstOrderGrdPhiScalableQuadrature;
+	delete firstOrderGrdPhiScalableQuadrature;
       if (secondOrderScalableQuadrature) 
-	DELETE secondOrderScalableQuadrature;
-    };
+	delete secondOrderScalableQuadrature;
+    }
 
     void scaleQuadratures(const SubElInfo& subElInfo);
 
     void getSubElementVector(SubElInfo *subElInfo, 
 			     const ElInfo *elInfo, 
-			     ElementVector *userVec);
+			     ElementVector& userVec);
 
     void getSubElementMatrix(SubElInfo *subElInfo, 
 			     const ElInfo *elInfo, 
-			     ElementMatrix *userMat);
+			     ElementMatrix& userMat);
 
     void getSubPolytopeVector(SubPolytope *subPolytope,
 			      SubElementAssembler *subElementAssembler,
 			      const ElInfo *elInfo,
-			      ElementVector *userVec);
+			      ElementVector& userVec);
 
     void getSubPolytopeMatrix(SubPolytope *subPolytope,
 			      SubElementAssembler *subElementAssembler,
 			      const ElInfo *elInfo,
-			      ElementMatrix *userMat);
+			      ElementMatrix& userMat);
 
   protected:
     ScalableQuadrature *zeroOrderScalableQuadrature;
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index 5c95bc5d89ec1cc6fc1cac1f4445e25a649d0c7e..a56671a53e6ca935605faaeef7c8d5ba46002468 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -29,8 +29,6 @@
 #include "ElInfo2d.h"
 #include "ElInfo3d.h"
 #include "Element.h"
-#include "ElementMatrix.h"
-#include "ElementVector.h"
 #include "Error.h"
 #include "Estimator.h"
 #include "FileWriter.h"
diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h
index 8e3c08e7a4134de9cd3a3b4ae6b9edbfc6af0b2b..ca613519d74248fa59ddeca0e089c7cc0ac42592 100644
--- a/AMDiS/src/AMDiS_fwd.h
+++ b/AMDiS/src/AMDiS_fwd.h
@@ -23,10 +23,13 @@
 #ifndef AMDIS_AMDIS_FWD_INCLUDE
 #define AMDIS_AMDIS_FWD_INCLUDE
 
+#include <boost/numeric/mtl/mtl.hpp>
+
 namespace AMDiS {
 
   class AdaptInfo;
   class AdaptStationary;
+  class Assembler;
   class BasisFunction; 
   class BoundaryManager;
   class CGSolver;
@@ -37,8 +40,6 @@ namespace AMDiS {
   class DOFIndexedBase;
   class DOFMatrix;
   class Element;
-  class ElementMatrix;
-  class ElementVector;
   class ElInfo; 
   class Estimator;
   class FastQuadrature;
@@ -54,6 +55,7 @@ namespace AMDiS {
   class Mesh; 
   class OEMSolver;
   class Operator;
+  class OperatorTerm;
   class ProblemInstat;
   class ProblemIterationInterface;
   class ProblemTimeInterface;
@@ -83,6 +85,12 @@ namespace AMDiS {
   template<typename T>                                 class WorldVector;
   template<typename T>                                 class WorldMatrix;
   template<typename T>                                 class VectorOfFixVecs;
+
+
+
+  typedef mtl::dense2D<double>              ElementMatrix;
+  typedef mtl::dense_vector<double>         ElementVector;
+
 } // namespace AMDiS
 
 #endif // AMDIS_AMDIS_FWD_INCLUDE
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 0df8b234bd2c4bb4f02d923a00a229740ab25a53..f8b164f9db897fcdbd162f9468b7c10ee13f0c92 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -1,13 +1,10 @@
 #include <vector>
 #include <algorithm>
 #include <boost/numeric/mtl/mtl.hpp>
-
 #include "Assembler.h"
 #include "Operator.h"
 #include "Element.h"
 #include "QPsiPhi.h"
-#include "ElementMatrix.h"
-#include "ElementVector.h"
 #include "DOFVector.h"
 #include "OpenMP.h"
 
@@ -24,28 +21,19 @@ namespace AMDiS {
       remember(true),
       rememberElMat(false),
       rememberElVec(false),
-      elementMatrix(NULL),
-      elementVector(NULL),
+      elementMatrix(nRow, nCol),
+      elementVector(nRow),
       lastMatEl(NULL),
       lastVecEl(NULL),
       lastTraverseId(-1)
   
-  {
-    elementMatrix = NEW ElementMatrix(nRow, nCol);
-    elementVector = NEW ElementVector(nRow);
-  }
+  {}
 
   Assembler::~Assembler()
-  {
-    if (elementMatrix)
-      DELETE elementMatrix;
-      
-    if (elementVector)
-      DELETE elementVector;
-  }
+  {}
 
   void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
-					 ElementMatrix *userMat,
+					 ElementMatrix& userMat,
 					 double factor)
   {
     FUNCNAME("Assembler::calculateElementMatrix()");
@@ -61,18 +49,18 @@ namespace AMDiS {
     }
 
     if (el != lastMatEl || !operat->isOptimized()) {
-      if (rememberElMat) {
-	elementMatrix->set(0.0);
-      }
+      if (rememberElMat)
+	set_to_zero(elementMatrix);
+
       lastMatEl = el;
     } else {
       if (rememberElMat) {
-	axpy(factor, *elementMatrix, *userMat);
+	userMat += factor * elementMatrix;
 	return;
       }
     }
  
-    ElementMatrix *mat = rememberElMat ? elementMatrix : userMat;
+    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
 
     if (secondOrderAssembler)
       secondOrderAssembler->calculateElementMatrix(elInfo, mat);
@@ -83,17 +71,15 @@ namespace AMDiS {
     if (zeroOrderAssembler)
       zeroOrderAssembler->calculateElementMatrix(elInfo, mat);
 
-    if (rememberElMat && userMat) {     
-      axpy(factor, *elementMatrix, *userMat);
-    }      
-
+    if (rememberElMat)
+      userMat += factor * elementMatrix;
   }
 
   void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
 					 const ElInfo *colElInfo,
 					 const ElInfo *smallElInfo,
 					 const ElInfo *largeElInfo,
-					 ElementMatrix *userMat,
+					 ElementMatrix& userMat,
 					 double factor)
   {
     FUNCNAME("Assembler::calculateElementMatrix()");
@@ -110,18 +96,18 @@ namespace AMDiS {
     }
 
     if (el != lastMatEl || !operat->isOptimized()) {
-      if (rememberElMat) {
-	elementMatrix->set(0.0);
-      }
+      if (rememberElMat)
+	set_to_zero(elementMatrix);
+
       lastMatEl = el;
     } else {
       if (rememberElMat) {
-	axpy(factor, *elementMatrix, *userMat);
+	userMat += factor * elementMatrix;
 	return;
       }
     }
   
-    ElementMatrix *mat = rememberElMat ? elementMatrix : userMat;
+    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
 
     if (secondOrderAssembler)
       secondOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, 
@@ -136,13 +122,12 @@ namespace AMDiS {
       zeroOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, 
 						 smallElInfo, largeElInfo, mat);
   
-    if (rememberElMat && userMat) {
-      axpy(factor, *elementMatrix, *userMat);
-    }      
+    if (rememberElMat)
+      userMat += factor * elementMatrix;
   }
 
   void Assembler::calculateElementVector(const ElInfo *elInfo, 
-					 ElementVector *userVec,
+					 ElementVector& userVec,
 					 double factor)
   {
     FUNCNAME("Assembler::calculateElementVector()");
@@ -158,40 +143,39 @@ namespace AMDiS {
     }
     
     if (el != lastVecEl || !operat->isOptimized()) {
-      if (rememberElVec) {
-	elementVector->set(0.0);
-      }
+      if (rememberElVec)
+	set_to_zero(elementVector);
+	
       lastVecEl = el;
     } else {
       if (rememberElVec) {
-	axpy(factor, *elementVector, *userVec);
+	userVec += factor * elementVector;
 	return;
       }
     }
-    ElementVector *vec = rememberElVec ? elementVector : userVec;
+
+    ElementVector& vec = rememberElVec ? elementVector : userVec;
     if (operat->uhOld && remember) {
       matVecAssemble(elInfo, vec);
       if (rememberElVec) {
-	axpy(factor, *elementVector, *userVec);
+	userVec += factor * elementVector;
       }
       return;
     } 
-    if (firstOrderAssemblerGrdPsi) {
+
+    if (firstOrderAssemblerGrdPsi)
       firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
-    }
-    if (zeroOrderAssembler) {
+    if (zeroOrderAssembler)
       zeroOrderAssembler->calculateElementVector(elInfo, vec);
-    }
-    if (rememberElVec) {
-      axpy(factor, *elementVector, *userVec);
-    }      
+    if (rememberElVec)
+      userVec += factor * elementVector;
   }
 
   void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
 					 const ElInfo *auxElInfo,
 					 const ElInfo *smallElInfo,
 					 const ElInfo *largeElInfo,
-					 ElementVector *userVec, 
+					 ElementVector& userVec, 
 					 double factor)
   {
     FUNCNAME("Assembler::calculateElementVector()");
@@ -207,17 +191,17 @@ namespace AMDiS {
     }
     
     if (el != lastVecEl || !operat->isOptimized()) {
-      if (rememberElVec) {
-	elementVector->set(0.0);
-      }
+      if (rememberElVec)
+	set_to_zero(elementVector);
+
       lastVecEl = el;
     } else {
       if (rememberElVec) {
-	axpy(factor, *elementVector, *userVec);
+	userVec += factor * elementVector;
 	return;
       }
     }
-    ElementVector *vec = rememberElVec ? elementVector : userVec;
+    ElementVector& vec = rememberElVec ? elementVector : userVec;
 
     if (operat->uhOld && remember) {
 
@@ -228,8 +212,9 @@ namespace AMDiS {
       }
 
       if (rememberElVec) {
-	axpy(factor, *elementVector, *userVec);
+	userVec += factor * elementVector;
       }
+
       return;
     } 
 
@@ -246,7 +231,7 @@ namespace AMDiS {
 
   }
 
-  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector *vec)
+  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
   {
     FUNCNAME("Assembler::matVecAssemble()");
 
@@ -257,16 +242,15 @@ namespace AMDiS {
 
     operat->uhOld->getLocalVector(el, uhOldLoc);
     
-    if (el != lastMatEl) {
-      calculateElementMatrix(elInfo, NULL);
-    }
+    if (el != lastMatEl)
+      calculateElementMatrix(elInfo, elementMatrix);
     
     for (int i = 0; i < nBasFcts; i++) {
       double val = 0.0;
       for (int j = 0; j < nBasFcts; j++) {
-	val += (*elementMatrix)[i][j] * uhOldLoc[j];
+	val += elementMatrix[i][j] * uhOldLoc[j];
       }
-      (*vec)[i] += val;
+      vec[i] += val;
     }
     
 
@@ -275,7 +259,7 @@ namespace AMDiS {
 
   void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
 				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
-				 ElementVector *vec)
+				 ElementVector& vec)
   {
     FUNCNAME("Assembler::matVecAssemble()");
 
@@ -306,15 +290,16 @@ namespace AMDiS {
 
     
     if (mainEl != lastMatEl) {
-      calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, NULL);
+      calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
+			     elementMatrix);
     }
     
     for (int i = 0; i < nBasFcts; i++) {
       double val = 0.0;
       for (int j = 0; j < nBasFcts; j++) {
-	val += (*elementMatrix)[i][j] * uhOldLoc2[j];
+	val += elementMatrix[i][j] * uhOldLoc2[j];
       }
-      (*vec)[i] += val;
+      vec[i] += val;
     }
     
 
@@ -336,45 +321,6 @@ namespace AMDiS {
       zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
   }
 
-  void Assembler::initElementMatrix(ElementMatrix *elMat, 
-				    const ElInfo *rowElInfo,
-				    const ElInfo *colElInfo)
-  {
-    TEST_EXIT_DBG(elMat)("No ElementMatrix!\n");
-
-    elMat->set(0.0);
-         
-    rowFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
-						   rowFESpace->getAdmin(),
-						   &(elMat->rowIndices));
-
-    if (rowFESpace == colFESpace) {
-      elMat->colIndices = elMat->rowIndices;
-    } else {
-      if (colElInfo) {
-	colFESpace->getBasisFcts()->getLocalIndicesVec(colElInfo->getElement(),
-						       colFESpace->getAdmin(),
-						       &(elMat->colIndices));
-      } else {
-	// If there is no colElInfo pointer, use rowElInfo the get the indices.
-	colFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
-						       colFESpace->getAdmin(),
-						       &(elMat->colIndices));
-      }
-    }
-  }
-
-  void Assembler::initElementVector(ElementVector *elVec, const ElInfo *elInfo)
-  {
-    TEST_EXIT_DBG(elVec)("No ElementVector!\n");
-
-    elVec->set(0.0);
- 
-    rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), 
-						   rowFESpace->getAdmin(), 
-						   &(elVec->dofIndices));
-  }
-
   void Assembler::checkQuadratures()
   { 
     if (secondOrderAssembler) {
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index 2f8e846da23d9e2817d868326ae05b7956071989..e7da738f49c577811130119b94f166dbaaa4d184 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -60,29 +60,21 @@ namespace AMDiS {
     /// Destructor
     ~Assembler();
 
-    void initElementMatrix(ElementMatrix *elMat, 
-			   const ElInfo *rowElInfo,
-			   const ElInfo *colElInfo = NULL);
-
-
-    void initElementVector(ElementVector *elVec,
-			   const ElInfo *elInfo);
-
     /// Assembles the element matrix for the given ElInfo
     void calculateElementMatrix(const ElInfo *elInfo, 
-				ElementMatrix *userMat, 
+				ElementMatrix& userMat, 
 				double factor = 1.0);
 
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *userMat,
+				ElementMatrix& userMat,
 				double factor = 1.0);
 
     /// Assembles the element vector for the given ElInfo
     void calculateElementVector(const ElInfo *elInfo, 
-				ElementVector *userVec, 
+				ElementVector& userVec, 
 				double factor = 1.0);
 
 
@@ -90,7 +82,7 @@ namespace AMDiS {
 				const ElInfo *auxElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementVector *userVec, 
+				ElementVector& userVec, 
 				double factor = 1.0);
 
     /// Returns \ref rowFESpace.
@@ -190,13 +182,13 @@ namespace AMDiS {
      * Vector assembling by element matrix-vector multiplication. 
      * Usefull if an element matrix was already calculated.
      */
-    void matVecAssemble(const ElInfo *elInfo, ElementVector *vec);
+    void matVecAssemble(const ElInfo *elInfo, ElementVector& vec);
 
 
     ///
     void matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
 			const ElInfo *smallElInfo, const ElInfo *largeElInfo,
-			ElementVector *vec);
+			ElementVector& vec);
 
     /** \brief
      * Checks whether quadratures for subassemblers are already set.
@@ -242,10 +234,10 @@ namespace AMDiS {
     bool rememberElVec;
 
     /// Locally stored element matrix
-    ElementMatrix* elementMatrix;
+    ElementMatrix elementMatrix;
 
     /// Locally stored element vector
-    ElementVector* elementVector;
+    ElementVector elementVector;
   
     /** \brief
      * Used to check whether \ref initElement() must be called, because
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index afa329059e83a67293b50806463a8a8b3a74758a..632551fbf0a7374d78dfb4c3f857705e7946deec 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -1,6 +1,6 @@
 #include <algorithm>
 #include <png.h>
-
+#include <boost/numeric/mtl/mtl.hpp>
 #include "DOFMatrix.h"
 #include "QPsiPhi.h"
 #include "BasisFunction.h"
@@ -13,12 +13,9 @@
 #include "Operator.h"
 #include "BoundaryCondition.h"
 #include "BoundaryManager.h"
-#include "ElementMatrix.h"
 #include "Assembler.h"
 #include "Utilities.h"
 
-#include <boost/numeric/mtl/mtl.hpp>
-
 namespace AMDiS {
 
   using namespace mtl;
@@ -28,7 +25,9 @@ namespace AMDiS {
   DOFMatrix::DOFMatrix()
     : rowFESpace(NULL),
       colFESpace(NULL),
-      elementMatrix(NULL),
+      elementMatrix(3, 3),
+      nRow(0),
+      nCol(0),
       inserter(NULL)
   {}
 
@@ -50,8 +49,12 @@ namespace AMDiS {
       (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->addDOFIndexed(this);
 
     boundaryManager = NEW BoundaryManager(rowFESpace_);
-    elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
-				      colFESpace->getBasisFcts()->getNumber());
+
+    nRow = rowFESpace->getBasisFcts()->getNumber();
+    nCol = colFESpace->getBasisFcts()->getNumber();
+    elementMatrix.change_dim(nRow, nCol);
+    rowIndices.resize(nRow);
+    colIndices.resize(nCol);
 
     applyDBCs.clear();
 
@@ -74,15 +77,10 @@ namespace AMDiS {
   DOFMatrix::~DOFMatrix()
   {
     FUNCNAME("DOFMatrix::~DOFMatrix()");
-    if (rowFESpace && rowFESpace->getAdmin()) {
+    if (rowFESpace && rowFESpace->getAdmin())
       (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
-    }  
 
-    if (boundaryManager)
-      DELETE boundaryManager;
-
-    if (elementMatrix)
-      DELETE elementMatrix;
+    if (boundaryManager) delete boundaryManager;
     if (inserter) delete inserter;
   }
 
@@ -164,19 +162,18 @@ namespace AMDiS {
       boundaryManager = NULL;
     }
 
-    if (rhs.elementMatrix) {
-      elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
-					colFESpace->getBasisFcts()->getNumber());
-    } else {
-      elementMatrix = NULL;
-    }
+    nRow = rhs.nRow;
+    nCol = rhs.nCol;
+    elementMatrix.change_dim(nRow, nCol);
 
     return *this;
   }
 
   void DOFMatrix::addElementMatrix(double sign, 
-				   const ElementMatrix &elMat, 
+				   const ElementMatrix& elMat, 
 				   const BoundaryType *bound,
+				   ElInfo* rowElInfo,
+				   ElInfo* colElInfo,
 				   bool add)
   {
     FUNCNAME("DOFMatrix::addElementMatrix()");
@@ -184,14 +181,33 @@ namespace AMDiS {
     TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode");
     inserter_type &ins= *inserter;
 
-    DegreeOfFreedom row, col;
-    double entry;
+    // === Get indices mapping from local to global matrix indices. ===
 
-    int n_row = elMat.rowIndices.getSize();
-    int n_col = elMat.colIndices.getSize();
+    rowFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
+						   rowFESpace->getAdmin(),
+						   &rowIndices);
+    if (rowFESpace == colFESpace) {
+      colIndices = rowIndices;
+    } else {
+      if (colElInfo) {
+	colFESpace->getBasisFcts()->getLocalIndicesVec(colElInfo->getElement(),
+						       colFESpace->getAdmin(),
+						       &colIndices);
+      } else {
+	// If there is no colElInfo pointer, use rowElInfo the get the indices.
+	colFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
+						       colFESpace->getAdmin(),
+						       &colIndices);
+      }
+    }
 
-    for (int i = 0; i < n_row; i++)  {   // for all rows of element matrix
-      row = elMat.rowIndices[i];
+
+    // === Add element matrix to the global matrix using the indices mapping. ===
+
+    DegreeOfFreedom row, col;
+    double entry;
+    for (int i = 0; i < nRow; i++)  {   // for all rows of element matrix
+      row = rowIndices[i];
 
       BoundaryCondition *condition = 
 	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
@@ -200,8 +216,8 @@ namespace AMDiS {
 	if (!coupleMatrix)
 	  applyDBCs.insert(static_cast<int>(row));
       } else
-	for (int j = 0; j < n_col; j++) {  // for all columns
-	  col = elMat.colIndices[j];
+	for (int j = 0; j < nCol; j++) {  // for all columns
+	  col = colIndices[j];
 	  entry = elMat[i][j];
 	  if (add)
 	    ins[row][col]+= sign * entry;
@@ -248,8 +264,7 @@ namespace AMDiS {
   {
     FUNCNAME("DOFMatrix::assemble()");
 
-    operators[0]->getAssembler(omp_get_thread_num())->
-      initElementMatrix(elementMatrix, elInfo);
+    set_to_zero(elementMatrix);
 
     std::vector<Operator*>::iterator it = operators.begin();
     std::vector<double*>::iterator factorIt = operatorFactor.begin();
@@ -261,7 +276,7 @@ namespace AMDiS {
       }
     }        
 
-    addElementMatrix(factor, *elementMatrix, bound);   
+    addElementMatrix(factor, elementMatrix, bound, elInfo, NULL);   
   }
 
   void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
@@ -271,9 +286,9 @@ namespace AMDiS {
 
       TEST_EXIT_DBG(op)("No operator!\n");
 
-      op->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
+      set_to_zero(elementMatrix);
       op->getElementMatrix(elInfo, elementMatrix, factor);
-      addElementMatrix(factor, *elementMatrix, bound);
+      addElementMatrix(factor, elementMatrix, bound, elInfo, NULL);
   }
 
   void DOFMatrix::assemble(double factor, 
@@ -283,16 +298,13 @@ namespace AMDiS {
   {
     FUNCNAME("DOFMatrix::assemble()");
 
-    if(!op && operators.size() == 0) {
+    if (!op && operators.size() == 0)
       return;
-    }
 
-    Operator *operat = op ? op : operators[0];
-    operat->getAssembler(omp_get_thread_num())->
-      initElementMatrix(elementMatrix, rowElInfo, colElInfo);
-    
+    set_to_zero(elementMatrix);
+
     if (op) {
-      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo,
+      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
 			   elementMatrix);
     } else {
       std::vector<Operator*>::iterator it = operators.begin();
@@ -305,7 +317,7 @@ namespace AMDiS {
       }      
     }
 
-    addElementMatrix(factor, *elementMatrix, bound);       
+    addElementMatrix(factor, elementMatrix, bound, rowElInfo, colElInfo);       
   }
 
   void DOFMatrix::assemble2(double factor, 
@@ -319,11 +331,9 @@ namespace AMDiS {
       return;
     }
 
-    Operator *operat = op ? op : operators[0];
-    operat->getAssembler(omp_get_thread_num())->
-      initElementMatrix(elementMatrix, mainElInfo);
+    set_to_zero(elementMatrix);
     
-    if(op) {
+    if (op) {
       ERROR_EXIT("TODO");
 //       op->getElementMatrix(rowElInfo, colElInfo, 
 // 			   smallElInfo, largeElInfo,
@@ -343,7 +353,7 @@ namespace AMDiS {
       }      
     }
 
-    addElementMatrix(factor, *elementMatrix, bound);       
+    addElementMatrix(factor, elementMatrix, bound, mainElInfo, NULL);       
   }
 
   void DOFMatrix::finishAssembling()
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index fcfc99a9768465b3badc409432c3be424de42936..04becace783d19870c03bc3cd52cf3ff4820ff1c 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -27,7 +27,6 @@
 #include <memory>
 #include <list>
 #include <boost/numeric/mtl/mtl.hpp>
-
 #include "AMDiS_fwd.h"
 #include "Global.h"
 #include "Flag.h"
@@ -210,6 +209,8 @@ namespace AMDiS {
     void addElementMatrix(double sign, 
 			  const ElementMatrix& elMat, 
 			  const BoundaryType *bound,
+			  ElInfo* rowElInfo,
+			  ElInfo* colElInfo,
 			  bool add = true);
 
     /* \brief
@@ -481,7 +482,19 @@ namespace AMDiS {
     bool coupleMatrix;
 
     /// Temporary variable used in assemble()
-    ElementMatrix *elementMatrix;
+    ElementMatrix elementMatrix;
+
+    /// Number of basis functions in the row fe space.
+    int nRow;
+
+    /// Number of basis functions in the col fe space.
+    int nCol;
+
+    /// Maps local row indices of an element to global matrix indices.
+    Vector<DegreeOfFreedom> rowIndices;
+
+    /// Maps local col indices of an element to global matrix indices.
+    Vector<DegreeOfFreedom> colIndices;
 
     /* \brief
      * A set of row indices. When assembling the DOF matrix, all rows, that
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 35fac1a3d9c00b1c339e04f91ecc3a95c09aba43..38be5c24d2d97aecc77567bc054270e71b2fc87b 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -827,14 +827,10 @@ namespace AMDiS {
   {
     FUNCNAME("DOFVector::assemble()");
 
-    TEST_EXIT_DBG(this->elementVector)("No ElementVector!\n");
-
     if (!(op || this->operators.size())) 
       return;
 
-    Operator *operat = op ? op : this->operators[0];
-
-    operat->getAssembler(omp_get_thread_num())->initElementVector(this->elementVector, elInfo);
+    set_to_zero(this->elementVector);
 
     if (op) {
       op->getElementVector(elInfo, this->elementVector);
@@ -852,7 +848,7 @@ namespace AMDiS {
       }
     }
 
-    addElementVector(factor, *this->elementVector, bound);
+    addElementVector(factor, this->elementVector, bound, elInfo);
   }
 
   template<>
@@ -867,10 +863,7 @@ namespace AMDiS {
     if (!(op || this->operators.size())) 
       return;
 
-    Operator *operat = op ? op : this->operators[0];
-
-    operat->getAssembler(omp_get_thread_num())->
-	initElementVector(this->elementVector, mainElInfo);
+    set_to_zero(this->elementVector);
 
     if (op) {
       ERROR_EXIT("TODO");
@@ -891,7 +884,7 @@ namespace AMDiS {
     }
   
 
-    addElementVector(factor, *this->elementVector, bound);
+    addElementVector(factor, this->elementVector, bound, mainElInfo);
   }
 
 
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 965f64257432eb8972d1203d34c087206d26ea91..ac16ebd991b0ea42633b21b342358914d3a57110 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -53,7 +53,7 @@ namespace AMDiS {
 
     DOFVectorBase() 
       : feSpace(NULL),
-	elementVector(NULL),
+	elementVector(3),
         boundaryManager(NULL),
         nBasFcts(0)
     {}
@@ -119,8 +119,9 @@ namespace AMDiS {
 		   Operator *op = NULL);
 
     void addElementVector(T sign,
-			  const ElementVector &elVec, 
+			  const ElementVector& elVec, 
 			  const BoundaryType *bound,
+			  ElInfo *elInfo,
 			  bool add = true); 
 
     /* \brief
@@ -207,7 +208,7 @@ namespace AMDiS {
     std::string name;
 
     ///
-    ElementVector *elementVector;
+    ElementVector elementVector;
 
     ///
     std::vector<Operator*> operators;
@@ -304,7 +305,6 @@ namespace AMDiS {
     /// Empty constructor. No initialization!
     DOFVector() 
       : DOFVectorBase<T>(), 
-	//elementVector(NULL), 
 	feSpace(NULL),
 	coarsenOperation(NO_OPERATION)
     {}
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 04cd7c624e0447cd7dfdca7d253f72f9d78d6239..bcfb853bb02e06991d41aabf5c5f45ebbd774328 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -2,6 +2,11 @@
 #include <algorithm>
 #include <math.h>
 
+#include <boost/numeric/mtl/mtl.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
 #include "FixVec.h"
 #include "Boundary.h"
 #include "DOFAdmin.h"
@@ -13,19 +18,12 @@
 #include "Quadrature.h"
 #include "AbstractFunction.h"
 #include "BoundaryManager.h"
-#include "ElementVector.h"
 #include "Assembler.h"
 #include "OpenMP.h"
 #include "Operator.h"
 #include "Parameters.h"
 #include "Traverse.h"
 
-#include <boost/numeric/mtl/mtl.hpp>
-#include <boost/numeric/mtl/utility/tag.hpp>
-#include <boost/numeric/mtl/utility/category.hpp>
-#include <boost/numeric/linear_algebra/identity.hpp>
-
-
 // Defining the interface for MTL4
 namespace mtl {
 
@@ -73,7 +71,7 @@ namespace AMDiS {
   DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
     : feSpace(f),
       name(n),
-      elementVector(NULL),
+      elementVector(f->getBasisFcts()->getNumber()),
       boundaryManager(NULL)
   {
     nBasFcts = feSpace->getBasisFcts()->getNumber();
@@ -92,8 +90,6 @@ namespace AMDiS {
       grdTmp[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
       D2Phis[i] = NEW DimMat<double>(dim, NO_INIT);
     }
-
-    elementVector = NEW ElementVector(this->nBasFcts);
   }
   
   template<typename T>
@@ -111,7 +107,6 @@ namespace AMDiS {
   template<typename T>
   DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
     : DOFVectorBase<T>(f, n),
-      //elementVector(NULL), 
       coarsenOperation(NO_OPERATION)
   {
     init(f, n);
@@ -157,18 +152,22 @@ namespace AMDiS {
   void DOFVectorBase<T>::addElementVector(T factor, 
 					  const ElementVector &elVec, 
 					  const BoundaryType *bound,
+					  ElInfo *elInfo,
 					  bool add)
   {
     FUNCNAME("DOFVector::addElementVector()");
 
-    int n_row = elVec.getSize();
+    Vector<DegreeOfFreedom> indices(nBasFcts);
+    feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), 
+						feSpace->getAdmin(),
+						&indices);
 
-    for (DegreeOfFreedom i = 0; i < n_row; i++) {
+    for (DegreeOfFreedom i = 0; i < nBasFcts; i++) {
       BoundaryCondition *condition = 
 	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;
 
       if (!(condition && condition->isDirichlet())) {
-	DegreeOfFreedom irow = elVec.dofIndices[i];
+	DegreeOfFreedom irow = indices[i];
 	(*this)[irow] = (add ? (*this)[irow] : 0.0);
 	(*this)[irow] += factor * elVec[i];
       }
@@ -700,7 +699,7 @@ namespace AMDiS {
     DOFVectorBase<T>::feSpace = rhs.feSpace;
     this->nBasFcts = rhs.nBasFcts;
     vec = rhs.vec;
-    this->elementVector = new ElementVector(this->nBasFcts);
+    this->elementVector.change_dim(this->nBasFcts);
     interFct = rhs.interFct;
     coarsenOperation = rhs.coarsenOperation;
     this->operators = rhs.operators;
@@ -1042,12 +1041,9 @@ namespace AMDiS {
     static T* localVec = NULL;
     const DOFAdmin* admin = feSpace->getAdmin();
 
-    int i;
-    int nBasFcts = feSpace->getBasisFcts()->getNumber();
-
     T *result;
   
-    if(d) {
+    if (d) {
       result = d;
     } else {
       if(localVec) delete [] localVec;
@@ -1058,9 +1054,8 @@ namespace AMDiS {
     const DegreeOfFreedom *localIndices = 
       feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);
 
-    for(i = 0; i < nBasFcts; i++) {
+    for (int i = 0; i < nBasFcts; i++)
       result[i] = (*this)[localIndices[i]];
-    }
 
     return result;
   }
diff --git a/AMDiS/src/ElementMatrix.h b/AMDiS/src/ElementMatrix.h
deleted file mode 100644
index f29b9d6b90c576552b63f3d3635ff334aa0edb1a..0000000000000000000000000000000000000000
--- a/AMDiS/src/ElementMatrix.h
+++ /dev/null
@@ -1,55 +0,0 @@
-// ============================================================================
-// ==                                                                        ==
-// == AMDiS - Adaptive multidimensional simulations                          ==
-// ==                                                                        ==
-// ============================================================================
-// ==                                                                        ==
-// ==  crystal growth group                                                  ==
-// ==                                                                        ==
-// ==  Stiftung caesar                                                       ==
-// ==  Ludwig-Erhard-Allee 2                                                 ==
-// ==  53175 Bonn                                                            ==
-// ==  germany                                                               ==
-// ==                                                                        ==
-// ============================================================================
-// ==                                                                        ==
-// ==  http://www.caesar.de/cg/AMDiS                                         ==
-// ==                                                                        ==
-// ============================================================================
-
-/** \file ElementMatrix.h */
-
-#ifndef AMDIS_ELEMENTMATRIX_H
-#define AMDIS_ELEMENTMATRIX_H
-
-#include "MatrixVector.h"
-
-namespace AMDiS {
-
-  /** \ingroup Assembler
-   *  
-   * \brief
-   * Stores a single element matrix and the corresponding row and column 
-   * dof indices.
-   */
-  class ElementMatrix : public Matrix<double>
-  {
-  public:
-    /// Constructor
-    ElementMatrix(int numRows, int numCols) 
-      : Matrix<double>(numRows, numCols),
-	rowIndices(numRows),
-	colIndices(numCols)
-    {}
-
-  public:
-    /// Row dof indices.
-    Vector<DegreeOfFreedom> rowIndices;
-
-    /// Column dof indices.
-    Vector<DegreeOfFreedom> colIndices;
-  };
-
-}
-
-#endif // AMDIS_ELEMENTMATRIX_H
diff --git a/AMDiS/src/ElementVector.h b/AMDiS/src/ElementVector.h
deleted file mode 100644
index 4383ac71e868ea254bbd65ad2c15b4cc9099e1ff..0000000000000000000000000000000000000000
--- a/AMDiS/src/ElementVector.h
+++ /dev/null
@@ -1,50 +0,0 @@
-// ============================================================================
-// ==                                                                        ==
-// == AMDiS - Adaptive multidimensional simulations                          ==
-// ==                                                                        ==
-// ============================================================================
-// ==                                                                        ==
-// ==  crystal growth group                                                  ==
-// ==                                                                        ==
-// ==  Stiftung caesar                                                       ==
-// ==  Ludwig-Erhard-Allee 2                                                 ==
-// ==  53175 Bonn                                                            ==
-// ==  germany                                                               ==
-// ==                                                                        ==
-// ============================================================================
-// ==                                                                        ==
-// ==  http://www.caesar.de/cg/AMDiS                                         ==
-// ==                                                                        ==
-// ============================================================================
-
-/** \file ElementVector.h */
-
-#ifndef AMDIS_ELEMENTVECTOR_H
-#define AMDIS_ELEMENTVECTOR_H
-
-#include "MatrixVector.h"
-
-namespace AMDiS {
-
-  /** \ingroup Assembler
-   *  
-   * \brief
-   * Stores a single element vector and the corresponding dof indices.
-   */
-  class ElementVector : public Vector<double>
-  {
-  public:
-    /// Constructor.
-    ElementVector(int size) 
-      : Vector<double>(size),
-	dofIndices(size)
-    {}
-
-  public:
-    /// dof indices.
-    Vector<DegreeOfFreedom> dofIndices;
-  };
-
-}
-
-#endif // AMDIS_ELEMENTVECTOR_H
diff --git a/AMDiS/src/FiniteElemSpace.h b/AMDiS/src/FiniteElemSpace.h
index f09f1384fba1153cd1deca9235a8bff538fa3d19..243d59a89fa408cf6e87ae15fe4bcbe74edb8d26 100644
--- a/AMDiS/src/FiniteElemSpace.h
+++ b/AMDiS/src/FiniteElemSpace.h
@@ -28,7 +28,7 @@
 
 #include <string>
 #include <vector>
-#include "MemoryManager.h"
+#include "AMDiS_fwd.h"
 
 namespace AMDiS {
 
@@ -46,8 +46,6 @@ namespace AMDiS {
   class FiniteElemSpace
   {
   public:
-    MEMORY_MANAGED(FiniteElemSpace);
-
     /// Create an empty fe space.
     FiniteElemSpace();
 
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index f6a951b6c67b1fcaa26b116e46803668cce2ec8d..03bb3c2c5f553b887fc910e074861d356d2ef111 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -1,6 +1,5 @@
 #include <vector>
 #include <boost/numeric/mtl/mtl.hpp>
-
 #include "Assembler.h"
 #include "FirstOrderAssembler.h"
 #include "Operator.h"
@@ -8,7 +7,6 @@
 #include "FiniteElemSpace.h"
 #include "Quadrature.h"
 #include "DOFVector.h"
-#include "ElementMatrix.h"
 #include "OpenMP.h"
 
 namespace AMDiS {
@@ -113,7 +111,7 @@ namespace AMDiS {
     : FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
   {}
 
-  void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
 
@@ -143,7 +141,7 @@ namespace AMDiS {
       for (int i = 0; i < nRow; i++) {
 	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
 	for (int j = 0; j < nCol; j++) {
-	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
+	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
 	}
       }
     }
@@ -153,7 +151,7 @@ namespace AMDiS {
 				       const ElInfo *colElInfo,
 				       const ElInfo *smallElInfo,
 				       const ElInfo *largeElInfo,
-				       ElementMatrix *mat)
+				       ElementMatrix& mat)
   {
     FUNCNAME("Stand10::calculateElementMatrix()");
 
@@ -195,14 +193,14 @@ namespace AMDiS {
 	    val1 += grdPsi;
 	  }
 
-	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * val1) *
+	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * val1) *
 	    (*(phi->getPhi(j)))(quadrature->getLambda(iq));
 	}
       }
     }
   }
 
-  void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
+  void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
   {
     DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
     const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
@@ -223,7 +221,7 @@ namespace AMDiS {
 
       for (int i = 0; i < nRow; i++) {
 	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
-	(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
+	vec[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
       }
     }
   }
@@ -234,7 +232,7 @@ namespace AMDiS {
   {}
 
 
-  void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     VectorOfFixVecs<DimVec<double> > *grdPsi;
 
@@ -272,12 +270,12 @@ namespace AMDiS {
 
       for (int i = 0; i < nRow; i++) {
 	for (int j = 0; j < nCol; j++)
-	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
+	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
       }
     }
   }
 
-  void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
+  void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
   {
     VectorOfFixVecs<DimVec<double> > *grdPsi;
 
@@ -312,7 +310,7 @@ namespace AMDiS {
       grdPsi = psiFast->getGradient(iq);
 
       for (int i = 0; i < nRow; i++) {
-	(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
+	vec[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
       }
     }
   }
@@ -323,7 +321,7 @@ namespace AMDiS {
   }
 
 
-  void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     const int *k;
     const double *values;
@@ -362,7 +360,7 @@ namespace AMDiS {
 	for (int m = 0; m < nEntries[i][j]; m++) {
 	  val += values[m] * Lb[0][k[m]];
 	}
-	(*mat)[i][j] += val;
+	mat[i][j] += val;
       }
     }
   }
@@ -372,7 +370,7 @@ namespace AMDiS {
     : FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
   {}
 
-  void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
 
@@ -401,7 +399,7 @@ namespace AMDiS {
       for (int i = 0; i < nRow; i++) {
 	double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
 	for (int j = 0; j < nCol; j++)
-	  (*mat)[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
+	  mat[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
       }
     } 
   }
@@ -410,7 +408,7 @@ namespace AMDiS {
 				       const ElInfo *colElInfo,
 				       const ElInfo *smallElInfo,
 				       const ElInfo *largeElInfo,
-				       ElementMatrix *mat)
+				       ElementMatrix& mat)
   {
     FUNCNAME("Stand01::calculateElementMatrix()");
 
@@ -448,7 +446,7 @@ namespace AMDiS {
 
 	  (*(phi->getGrdPhi(j)))(quadrature->getLambda(iq), grdPhi);
 
-	  (*mat)[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * val0) * grdPhi);
+	  mat[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * val0) * grdPhi);
 
 	}
       }
@@ -459,7 +457,7 @@ namespace AMDiS {
     : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
   {}
 
-  void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     VectorOfFixVecs<DimVec<double> > *grdPhi;
 
@@ -496,7 +494,7 @@ namespace AMDiS {
 
       for (int i = 0; i < nRow; i++) {
 	for (int j = 0; j < nCol; j++)
-	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
+	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
       }
     }
   }
@@ -506,7 +504,7 @@ namespace AMDiS {
   {
   }
 
-  void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     const int *l;
     const double *values;
@@ -545,12 +543,12 @@ namespace AMDiS {
 	for (int m = 0; m < nEntries[i][j]; m++) {
 	  val += values[m] * Lb[0][l[m]];
 	}
-	(*mat)[i][j] += val;
+	mat[i][j] += val;
       }
     }
   }
 
-  void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
+  void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
   {
     const int *k;
     const double *values;
@@ -588,7 +586,7 @@ namespace AMDiS {
       for (int m = 0; m < nEntries[i]; m++) {
 	val += values[m] * Lb[0][k[m]];
       }
-      (*vec)[i] += val;
+      vec[i] += val;
     }
   }
 
diff --git a/AMDiS/src/FirstOrderAssembler.h b/AMDiS/src/FirstOrderAssembler.h
index 0d8de4b9815dc37122ed45e77c8443a93a37e5a7..819bc77c9c97c36f784714309b4b8b11fcf3e391 100644
--- a/AMDiS/src/FirstOrderAssembler.h
+++ b/AMDiS/src/FirstOrderAssembler.h
@@ -18,8 +18,6 @@ namespace AMDiS {
   class FirstOrderAssembler : public SubAssembler
   {
   public:
-    MEMORY_MANAGED(FirstOrderAssembler);
-
     /** \brief
      * Creates and returns the FirstOrderAssembler for Operator op and
      * the given assembler. If all terms are piecewise constant precalculated 
@@ -81,17 +79,17 @@ namespace AMDiS {
     Stand10(Operator *op, Assembler *assembler, Quadrature *quad);
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     ///
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat);
+				ElementMatrix& mat);
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *, ElementVector */*vec*/);
+    void calculateElementVector(const ElInfo *, ElementVector&);
   };
 
 
@@ -110,17 +108,17 @@ namespace AMDiS {
     Stand01(Operator *op, Assembler *assembler, Quadrature *quad);
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     /// 
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat);
+				ElementMatrix& mat);
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *, ElementVector *) {
+    void calculateElementVector(const ElInfo*, ElementVector&) {
       ERROR_EXIT("should not be called\n");
     }
   };
@@ -141,20 +139,20 @@ namespace AMDiS {
     Quad10(Operator *op, Assembler *assembler, Quadrature *quad);
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     ///
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat) 
+				ElementMatrix& mat) 
     {
       ERROR_EXIT("CEM not yet\n");
     }
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *, ElementVector *);
+    void calculateElementVector(const ElInfo *, ElementVector&);
   };
 
   /**
@@ -172,20 +170,20 @@ namespace AMDiS {
     Quad01(Operator *op, Assembler *assembler, Quadrature *quad);
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     ///
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat) 
+				ElementMatrix& mat) 
     {
       ERROR_EXIT("CEM not yet\n");
     }
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *, ElementVector *) {
+    void calculateElementVector(const ElInfo*, ElementVector&) {
       ERROR_EXIT("should not be called\n");
     }
   };
@@ -206,20 +204,20 @@ namespace AMDiS {
     Pre10(Operator *op, Assembler *assembler, Quadrature *quad);
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     ///
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat) 
+				ElementMatrix& mat) 
     {
       ERROR_EXIT("CEM not yet\n");
     }
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *, ElementVector *);
+    void calculateElementVector(const ElInfo*, ElementVector&);
 
   protected:
     /// Integral of the product of the derivative of psi and phi.
@@ -247,20 +245,20 @@ namespace AMDiS {
     Pre01(Operator *op, Assembler *assembler, Quadrature *quad);
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     ///
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat) 
+				ElementMatrix& mat)
     {
       ERROR_EXIT("CEM not yet\n");
     }
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *, ElementVector *) {
+    void calculateElementVector(const ElInfo*, ElementVector&) {
       ERROR_EXIT("should not be called\n");
     }
 
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 35fb49f653915aea763d5372bad63ee58ca944a7..427bf07341f55bd135de71f7a10df9b813b86362 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -3,8 +3,6 @@
 #include "Assembler.h"
 #include "FixVec.h"
 #include "DOFVector.h"
-#include "ElementMatrix.h"
-#include "ElementVector.h"
 #include "Quadrature.h"
 #include "OpenMP.h"
 
@@ -304,7 +302,7 @@ namespace AMDiS {
   }
 
   void Operator::getElementMatrix(const ElInfo *elInfo, 
-				  ElementMatrix *userMat, 
+				  ElementMatrix& userMat, 
 				  double factor)
   {
     int myRank = omp_get_thread_num();
@@ -318,7 +316,7 @@ namespace AMDiS {
 
   void Operator::getElementMatrix(const ElInfo *rowElInfo, const ElInfo *colElInfo,
 				  const ElInfo *smallElInfo, const ElInfo *largeElInfo,
-				  ElementMatrix *userMat, 
+				  ElementMatrix& userMat, 
 				  double factor)
   {
     int myRank = omp_get_thread_num();
@@ -333,7 +331,7 @@ namespace AMDiS {
   }
 
   void Operator::getElementVector(const ElInfo *elInfo, 
-				  ElementVector *userVec, 
+				  ElementVector& userVec, 
 				  double factor)
   {
     int myRank = omp_get_thread_num();
@@ -347,7 +345,7 @@ namespace AMDiS {
 
   void Operator::getElementVector(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
 				  const ElInfo *smallElInfo, const ElInfo *largeElInfo,
-				  ElementVector *userVec,
+				  ElementVector& userVec,
 				  double factor)
   {
     int myRank = omp_get_thread_num();
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 236872019c820362a275270cce1b754bae145090..530a2744da89e6ed1346b916984bb99eb84b8108 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -23,9 +23,9 @@
 #define AMDIS_OPERATOR_H
 
 #include <vector>
+#include "AMDiS_fwd.h"
 #include "FixVec.h"
 #include "Flag.h"
-#include "MemoryManager.h"
 #include "MatrixVector.h"
 #include "ElInfo.h"
 #include "AbstractFunction.h"
@@ -34,20 +34,6 @@
 
 namespace AMDiS {
 
-  class Assembler;
-  class ElInfo;
-  class FiniteElemSpace;
-  class Operator;
-  class SubAssembler;
-  class ElementMatrix;
-  class ElementVector;
-  class Quadrature;
-  template<typename T> class DOFVectorBase;
-
-  // ============================================================================
-  // ===== class OperatorTerm ===================================================
-  // ============================================================================
-
   /** 
    * \ingroup Assembler
    * 
@@ -60,8 +46,6 @@ namespace AMDiS {
   class OperatorTerm
   {
   public:
-    MEMORY_MANAGED(OperatorTerm);
-
     /** \brief
      * Constructs an OperatorTerm with initially no properties.
      * degree_ is used to determine the degree of the needed quadrature
@@ -3292,10 +3276,6 @@ namespace AMDiS {
   };
 
 
-  // ============================================================================
-  // ===== class Operator =======================================================
-  // ============================================================================
-
   /** \brief
    * An Operator holds all information needed to assemble the system matrix
    * and the right hand side. It consists of four OperatorTerm lists each storing
@@ -3315,8 +3295,6 @@ namespace AMDiS {
   class Operator
   {
   public:
-    MEMORY_MANAGED(Operator);
-
     /** \brief
      * Constructs an empty Operator of type operatorType for the given 
      * FiniteElemSpace. 
@@ -3372,14 +3350,14 @@ namespace AMDiS {
      * factor to userMat.
      */
     virtual void getElementMatrix(const ElInfo *elInfo, 
-				  ElementMatrix *userMat, 
+				  ElementMatrix& userMat, 
 				  double factor = 1.0);
 
     virtual void getElementMatrix(const ElInfo *rowElInfo,
 				  const ElInfo *colElInfo,
 				  const ElInfo *smallElInfo,
 				  const ElInfo *largeElInfo,
-				  ElementMatrix *userMat,
+				  ElementMatrix& userMat,
 				  double factor = 1.0);
 
     /** \brief
@@ -3387,14 +3365,14 @@ namespace AMDiS {
      * factor to userVec.
      */
     virtual void getElementVector(const ElInfo *elInfo, 
-				  ElementVector *userVec, 
+				  ElementVector& userVec, 
 				  double factor = 1.0);
 
     virtual void getElementVector(const ElInfo *mainElInfo, 
 				  const ElInfo *auxElInfo,
 				  const ElInfo *smallElInfo,
 				  const ElInfo *largeElInfo,
-				  ElementVector *userVec,
+				  ElementVector& userVec,
 				  double factor = 1.0);
 
     /// That function must be called after one assembling cycle has been finished.
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index fd932bb7d6665ef2974c448354f9824adf1c9a7b..eea1e692c7e4f7643c82393ea233eac4b5fd71c0 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -1,6 +1,5 @@
 #include <vector>
 #include <boost/numeric/mtl/mtl.hpp>
-
 #include "Assembler.h"
 #include "SecondOrderAssembler.h"
 #include "Operator.h"
@@ -8,7 +7,6 @@
 #include "FiniteElemSpace.h"
 #include "Quadrature.h"
 #include "DOFVector.h"
-#include "ElementMatrix.h"
 #include "OpenMP.h"
 
 namespace AMDiS {
@@ -99,12 +97,12 @@ namespace AMDiS {
   Pre2::~Pre2()
   {
     for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
-      DELETE *(tmpLALt[i]);
-      DELETE tmpLALt[i];
+      delete *(tmpLALt[i]);
+      delete tmpLALt[i];
     }
   }
 
-  void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     const int **nEntries;
     const int *k, *l;
@@ -133,7 +131,7 @@ namespace AMDiS {
 	for (int m = 0; m < nEntries[i][i]; m++) {
 	  val += values[m] * tmpMat[k[m]][l[m]];
 	}
-	(*mat)[i][i] += val;
+	mat[i][i] += val;
 
 	for (int j = i + 1; j < nCol; j++) {
 	  k = q11->getKVec(i, j);
@@ -143,8 +141,8 @@ namespace AMDiS {
 	  for (int m = 0; m < nEntries[i][j]; m++) {
 	    val += values[m] * tmpMat[k[m]][l[m]];
 	  }
-	  (*mat)[i][j] += val;
-	  (*mat)[j][i] += val;
+	  mat[i][j] += val;
+	  mat[j][i] += val;
 	}
       }
     } else {  /*  A not symmetric or psi != phi        */
@@ -157,7 +155,7 @@ namespace AMDiS {
 	  for (int m = 0; m < nEntries[i][j]; m++) {
 	    val += values[m] * tmpMat[k[m]][l[m]];
 	  }
-	  (*mat)[i][j] += val;
+	  mat[i][j] += val;
 	}
       }
     }
@@ -175,14 +173,14 @@ namespace AMDiS {
       int nPoints = quadrature->getNumPoints();
       for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
 	for (int j = 0; j < nPoints; j++) {
-	  DELETE tmpLALt[i][j];
+	  delete tmpLALt[i][j];
 	}
-	DELETE [] tmpLALt[i];
+	delete [] tmpLALt[i];
       }
     }
   }
 
-  void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     int nPoints = quadrature->getNumPoints();
     int myRank = omp_get_thread_num();
@@ -222,14 +220,14 @@ namespace AMDiS {
 	grdPhi = phiFast->getGradient(iq);
 
 	for (int i = 0; i < nRow; i++) {
-	  (*mat)[i][i] += quadrature->getWeight(iq) * 
+	  mat[i][i] += quadrature->getWeight(iq) * 
 	    xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[i]);
 
 	  for (int j = i + 1; j < nCol; j++) {
 	    double val = quadrature->getWeight(iq) * 
 	      xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
-	    (*mat)[i][j] += val;
-	    (*mat)[j][i] += val;
+	    mat[i][j] += val;
+	    mat[j][i] += val;
 	  }
 	}
       }
@@ -242,7 +240,7 @@ namespace AMDiS {
 
 	for (int i = 0; i < nRow; i++) {
 	  for (int j = 0; j < nCol; j++) {
-	    (*mat)[i][j] += quadrature->getWeight(iq) * 
+	    mat[i][j] += quadrature->getWeight(iq) * 
 	      xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
 	  }
 	}
@@ -254,7 +252,7 @@ namespace AMDiS {
     : SecondOrderAssembler(op, assembler, quad, false)
   {}
 
-  void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     DimVec<double> grdPsi(dim, NO_INIT);
     VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
@@ -288,14 +286,14 @@ namespace AMDiS {
 
 	for (int i = 0; i < nRow; i++) {
 	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
-	  (*mat)[i][i] += quadrature->getWeight(iq) * 
+	  mat[i][i] += quadrature->getWeight(iq) * 
 	    (grdPsi * ((*LALt[iq]) * grdPhi[i]));
 
 	  for (int j = i + 1; j < nCol; j++) {
 	    double val = quadrature->getWeight(iq) * 
 	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
-	    (*mat)[i][j] += val;
-	    (*mat)[j][i] += val;
+	    mat[i][j] += val;
+	    mat[j][i] += val;
 	  }
 	}
       }
@@ -310,7 +308,7 @@ namespace AMDiS {
 	for (int i = 0; i < nRow; i++) {
 	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
 	  for (int j = 0; j < nCol; j++) {
-	    (*mat)[i][j] += quadrature->getWeight(iq) *
+	    mat[i][j] += quadrature->getWeight(iq) *
 	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
 	  }
 	}
@@ -318,16 +316,16 @@ namespace AMDiS {
     }
 
     for (int iq = 0; iq < nPoints; iq++) 
-      DELETE LALt[iq];
+      delete LALt[iq];
 
-    DELETE [] LALt;
+    delete [] LALt;
   }
 
   void Stand2::calculateElementMatrix(const ElInfo *rowElInfo, 
 				      const ElInfo *colElInfo,
 				      const ElInfo *smallElInfo,
 				      const ElInfo *largeElInfo,
-				      ElementMatrix *mat)
+				      ElementMatrix& mat)
   {
     FUNCNAME("Stand2::calculateElementMatrix()");
 
@@ -403,7 +401,7 @@ namespace AMDiS {
 
 	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
 
-	  (*mat)[i][j] += quadrature->getWeight(iq) *
+	  mat[i][j] += quadrature->getWeight(iq) *
 	      (grdPsi * ((*LALt[iq]) * val1));
 	}
       }
@@ -421,22 +419,16 @@ namespace AMDiS {
 	for (int i = 0; i < nRow; i++) {
 	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
 	  for (int j = 0; j < nCol; j++) {
-	    (*mat)[i][j] += quadrature->getWeight(iq) *
+	    mat[i][j] += quadrature->getWeight(iq) * 
 	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
 	  }
 	}
       }
 
-      if (p) {
-	mat->print();
-	WAIT_REALLY;
-      }
-
-
     for (int iq = 0; iq < nPoints; iq++) 
-      DELETE LALt[iq];
+      delete LALt[iq];
 
-    DELETE [] LALt;
+    delete [] LALt;
   }
 
 }
diff --git a/AMDiS/src/SecondOrderAssembler.h b/AMDiS/src/SecondOrderAssembler.h
index 0e2f113bb21af189b55bde3012f3f574916e3bb6..fed24e08260febb36f24acc3b649ec79d17b253a 100644
--- a/AMDiS/src/SecondOrderAssembler.h
+++ b/AMDiS/src/SecondOrderAssembler.h
@@ -16,8 +16,6 @@ namespace AMDiS {
   class SecondOrderAssembler : public SubAssembler
   {
   public:
-    MEMORY_MANAGED(SecondOrderAssembler);
-
     /** \brief
      * Creates and returns the SecondOrderAssembler for Operator op and
      * the given assembler. If all terms are piecewise constant precalculated 
@@ -58,23 +56,21 @@ namespace AMDiS {
   class Stand2 : public SecondOrderAssembler 
   {
   public:
-    MEMORY_MANAGED(Stand2);
-
     /// Constructor.
     Stand2(Operator *op, Assembler *assembler, Quadrature *quad);
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     ///
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat);
+				ElementMatrix& mat);
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *, ElementVector *) {
+    void calculateElementVector(const ElInfo *, ElementVector&) {
       ERROR_EXIT("should not be called\n");
     }
   };
@@ -88,28 +84,26 @@ namespace AMDiS {
   class Quad2 : public SecondOrderAssembler 
   {
   public:
-    MEMORY_MANAGED(Quad2);
-
     /// Constructor.
     Quad2(Operator *op, Assembler *assembler, Quadrature *quad);
 
     ~Quad2();
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     ///
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat) 
+				ElementMatrix& mat) 
     {
       ERROR_EXIT("CEM not yet\n");
     }
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *, ElementVector */*vec*/) {
+    void calculateElementVector(const ElInfo *, ElementVector&) {
       ERROR_EXIT("should not be called\n");
     }
 
@@ -130,8 +124,6 @@ namespace AMDiS {
   class Pre2 : public SecondOrderAssembler 
   {
   public:
-    MEMORY_MANAGED(Pre2);
-
     /// Constructor.
     Pre2(Operator *op, Assembler *assembler, Quadrature *quad);
 
@@ -139,20 +131,20 @@ namespace AMDiS {
     ~Pre2();
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     ///
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat) 
+				ElementMatrix& mat) 
     {
       ERROR_EXIT("Not yet implemented!\n");
     }
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *, ElementVector *) {
+    void calculateElementVector(const ElInfo *, ElementVector&) {
       ERROR_EXIT("should not be called\n");
     }
 
diff --git a/AMDiS/src/SubAssembler.h b/AMDiS/src/SubAssembler.h
index a56dbd4444e2ca2cc3eb11720d7c90311e7d44fe..4eaf222dccb93ea56e6207e5f194be0a17ac9372 100644
--- a/AMDiS/src/SubAssembler.h
+++ b/AMDiS/src/SubAssembler.h
@@ -3,26 +3,13 @@
 
 #include <map>
 #include <vector>
+#include "AMDiS_fwd.h"
 #include "FixVec.h"
 #include "Flag.h"
-#include "MemoryManager.h"
 #include "OpenMP.h"
 
 namespace AMDiS {
 
-  class Assembler;
-  class Quadrature;
-  class FastQuadrature;
-  class Operator;
-  class OperatorTerm;
-  class BasisFunction;
-  class FiniteElemSpace;
-  class ElInfo;
-  class ElementMatrix;
-  class ElementVector;
-
-  template<typename T> class DOFVectorBase;
-
   /**
    * \ingroup Assembler
    * 
@@ -39,8 +26,6 @@ namespace AMDiS {
   class SubAssembler
   {
   public:
-    MEMORY_MANAGED(SubAssembler);
-
     /** \brief
      * Creates a SubAssembler belonging to assembler for the terms of given 
      * order of Operator op. If order is equal to one, type spezifies what kind 
@@ -62,20 +47,20 @@ namespace AMDiS {
      * for mat must be provided by the caller.
      */
     virtual void calculateElementMatrix(const ElInfo *elInfo,
-					ElementMatrix *mat) = 0;
+					ElementMatrix& mat) = 0;
 
     virtual void calculateElementMatrix(const ElInfo *rowElInfo,
 					const ElInfo *colElInfo,
 					const ElInfo *smallElInfo,
 					const ElInfo *largeElInfo,
-					ElementMatrix *mat) = 0;
+					ElementMatrix& mat) = 0;
 
     /** \brief
      * Calculates the element vector for elInfo and adds it to vec. Memory
      * for vec must be provided by the caller.
      */
     virtual void calculateElementVector(const ElInfo *elInfo,
-					ElementVector *vec) = 0;
+					ElementVector& vec) = 0;
 
     /// Returns \ref terms
     inline std::vector<OperatorTerm*> *getTerms() { 
diff --git a/AMDiS/src/SurfaceAssembler.h b/AMDiS/src/SurfaceAssembler.h
index 9c95a1f9ec4238aa77e324678a414480650de1b8..35fa9b17d417a8e6ca6020f2f4e2ba3133034dd8 100644
--- a/AMDiS/src/SurfaceAssembler.h
+++ b/AMDiS/src/SurfaceAssembler.h
@@ -26,15 +26,10 @@
 #include "Mesh.h"
 #include "Assembler.h"
 #include "SubQuadrature.h"
-#include "ElementMatrix.h"
 #include "QPInfo.h"
 
 namespace AMDiS {
 
-  // ============================================================================
-  // ===== class SurfaceAssembler ================================================
-  // ============================================================================
-
   /** 
    * \ingroup Integration
    *
@@ -43,8 +38,6 @@ namespace AMDiS {
   class SurfaceAssembler : public Assembler
   {
   public:
-    MEMORY_MANAGED(SurfaceAssembler);
-
     /// Creates a SurfaceAssembler conforming to operat for the given \ref coords.
     SurfaceAssembler(Operator *operat,
 		     const FiniteElemSpace *rowFESpace,
@@ -136,9 +129,8 @@ namespace AMDiS {
 	FixVec<WorldVector<double>, VERTEX> worldCoords(rowDim_ - 1, NO_INIT);
 
 	// transform barycentric coords to world coords
-	for (int i = 0; i < rowDim_; i++) {
+	for (int i = 0; i < rowDim_; i++)
 	  elInfo->coordToWorld(coords_[i], &worldCoords[i]);
-	}
 
 	// set determinant for world coords of the side
 	det_ = elInfo->calcDet(worldCoords);
diff --git a/AMDiS/src/SurfaceOperator.h b/AMDiS/src/SurfaceOperator.h
index 915537d3fdbd15b67e036e0d617a263faf905923..c939d1c53e43bfb37b4b5d5eb475a593505b1bd4 100644
--- a/AMDiS/src/SurfaceOperator.h
+++ b/AMDiS/src/SurfaceOperator.h
@@ -26,15 +26,10 @@
 #include "Mesh.h"
 #include "Operator.h"
 #include "SurfaceQuadrature.h"
-#include "ElementMatrix.h"
 #include "OpenMP.h"
 
 namespace AMDiS {
 
-  // ============================================================================
-  // ===== class SurfaceOperator ================================================
-  // ============================================================================
-
   /** 
    * \ingroup Integration
    *
@@ -50,8 +45,6 @@ namespace AMDiS {
   class SurfaceOperator : public Operator
   {
   public:
-    MEMORY_MANAGED(SurfaceOperator);
-
     /// Creates a SurfaceOperator conforming to operat for the given \ref coords.
     SurfaceOperator(Operator *operat, 
 		    VectorOfFixVecs<DimVec<double> > &coords) 
@@ -130,7 +123,7 @@ namespace AMDiS {
      * the base class function.
      */
     virtual void getElementMatrix(const ElInfo *elInfo, 
-				  ElementMatrix *userMat, 
+				  ElementMatrix& userMat, 
 				  double factor = 1.0)
     {
       int i;
@@ -160,19 +153,17 @@ namespace AMDiS {
      * the base class function.
      */
     virtual void getElementVector(const ElInfo *elInfo, 
-				  ElementVector *userVec, 
+				  ElementVector& userVec, 
 				  double factor = 1.0) 
     {
-      int i;
       int dim = rowFESpace->getMesh()->getDim();
       double origDet = elInfo->getDet();
 
       FixVec<WorldVector<double>, VERTEX> worldCoords(dim-1, NO_INIT);
 
       // transform barycentric coords to world coords
-      for(i = 0; i < dim; i++) {
+      for (int i = 0; i < dim; i++)
 	elInfo->coordToWorld(coords_[i], worldCoords[i]);
-      }
 
       // set determinant for world coords of the side
       const_cast<ElInfo*>(elInfo)->setDet(elInfo->calcDet(worldCoords));
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index b5b5ebd61c7493c90fc8b5bcd4879f2b3ec11818..8d040afe6be1a9b8afdc8d66df669c7e1fe36abc 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -1,6 +1,5 @@
 #include <vector>
 #include <boost/numeric/mtl/mtl.hpp>
-
 #include "Assembler.h"
 #include "ZeroOrderAssembler.h"
 #include "Operator.h"
@@ -8,7 +7,6 @@
 #include "FiniteElemSpace.h"
 #include "Quadrature.h"
 #include "DOFVector.h"
-#include "ElementMatrix.h"
 #include "OpenMP.h"
 
 namespace AMDiS {
@@ -84,7 +82,7 @@ namespace AMDiS {
     : ZeroOrderAssembler(op, assembler, quad, false)
   {}
 
-  void StandardZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void StandardZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
     const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
@@ -112,11 +110,11 @@ namespace AMDiS {
 
 	for (int i = 0; i < nRow; i++) {
 	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
-	  (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i];
+	  mat[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i];
 	  for (int j = i + 1; j < nCol; j++) {
 	    double val = quadrature->getWeight(iq) * c[iq] * psival * phival[j];
-	    (*mat)[i][j] += val;
-	    (*mat)[j][i] += val;
+	    mat[i][j] += val;
+	    mat[j][i] += val;
 	  }
 	}
       }
@@ -132,7 +130,7 @@ namespace AMDiS {
 	for (int i = 0; i < nRow; i++) {
 	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
 	  for (int j = 0; j < nCol; j++) {
-	    (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
+	    mat[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
 	  }
 	}
       }
@@ -143,7 +141,7 @@ namespace AMDiS {
 					   const ElInfo *colElInfo,
 					   const ElInfo *smallElInfo,
 					   const ElInfo *largeElInfo,
-					   ElementMatrix *mat) 
+					   ElementMatrix& mat) 
   {
     FUNCNAME("StandardZOA::calculateElementMatrix()");
 
@@ -160,14 +158,14 @@ namespace AMDiS {
 	tmpMat[i][j] = 0.0;
 
 	for (int k = 0; k < nRow; k++) {
-	  tmpMat[i][j] += m[i][j] * (*mat)[j][i];
+	  tmpMat[i][j] += m[i][j] * mat[j][i];
 	}
       }
     }
 
     for (int i = 0; i < nRow; i++) {
       for (int j = 0; j < nRow; j++) {
-	(*mat)[i][j] = tmpMat[i][j];
+	mat[i][j] = tmpMat[i][j];
       }
     }
 
@@ -197,7 +195,7 @@ namespace AMDiS {
 
 	  val0 *= val1;
 
-	  (*mat)[i][j] += val0;
+	  mat[i][j] += val0;
 	}
       }
     }
@@ -205,7 +203,7 @@ namespace AMDiS {
 
   }
 
-  void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
+  void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
   {
     int nPoints = quadrature->getNumPoints();
     std::vector<double> c(nPoints, 0.0);
@@ -222,7 +220,7 @@ namespace AMDiS {
       for (int i = 0; i < nRow; i++) {
 	double psi = (*(owner->getRowFESpace()->getBasisFcts()->getPhi(i)))
 	  (quadrature->getLambda(iq));
-	(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi;
+	vec[i] += quadrature->getWeight(iq) * c[iq] * psi;
       }
     }
   }
@@ -231,7 +229,7 @@ namespace AMDiS {
     : ZeroOrderAssembler(op, assembler, quad, true)
   {}
 
-  void FastQuadZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void FastQuadZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     int nPoints = quadrature->getNumPoints();
     int myRank = omp_get_thread_num();
@@ -264,11 +262,11 @@ namespace AMDiS {
 	const double *psi = psiFast->getPhi(iq);
 	const double *phi = phiFast->getPhi(iq);
 	for (int i = 0; i < nRow; i++) {
-	  (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
+	  mat[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
 	  for (int j = i + 1; j < nCol; j++) {
 	    double val = quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
-	    (*mat)[i][j] += val;
-	    (*mat)[j][i] += val;
+	    mat[i][j] += val;
+	    mat[j][i] += val;
 	  }
 	}
       }
@@ -280,7 +278,7 @@ namespace AMDiS {
 	const double *phi = phiFast->getPhi(iq);
 	for (int i = 0; i < nRow; i++) {
 	  for (int j = 0; j < nCol; j++) {
-	    (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
+	    mat[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
 	  }
 	}
       }
@@ -291,7 +289,7 @@ namespace AMDiS {
 					   const ElInfo *colElInfo,
 					   const ElInfo *smallElInfo,
 					   const ElInfo *largeElInfo,
-					   ElementMatrix *mat) 
+					   ElementMatrix& mat) 
   {
     FUNCNAME("FastQuadZOA::calculateElementMatrix()");
 
@@ -336,14 +334,14 @@ namespace AMDiS {
 	  }
 	  val *= tmpval;
 
-	  (*mat)[j][i] += val;
+	  mat[j][i] += val;
 	}
       }
 
     }    
   }
 
-  void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
+  void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
   {
     int myRank = omp_get_thread_num();
     int nPoints = quadrature->getNumPoints();
@@ -372,7 +370,7 @@ namespace AMDiS {
 
       const double *psi = psiFast->getPhi(iq);
       for (int i = 0; i < nRow; i++) {
-	(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
+	vec[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
       }
     }
   }
@@ -382,7 +380,7 @@ namespace AMDiS {
   {
   }
 
-  void PrecalcZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
+  void PrecalcZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
   {
     if (firstCall) {
 #ifdef _OPENMP
@@ -410,23 +408,23 @@ namespace AMDiS {
 
     if (symmetric) {
       for (int i = 0; i < nRow; i++) {
-	(*mat)[i][i] += c[0] * q00->getValue(i,i);
+	mat[i][i] += c[0] * q00->getValue(i,i);
 	for (int j = i + 1; j < nCol; j++) {
 	  double val = c[0] * q00->getValue(i, j);
-	  (*mat)[i][j] += val;
-	  (*mat)[j][i] += val;
+	  mat[i][j] += val;
+	  mat[j][i] += val;
 	}
       }
     } else {
       for (int i = 0; i < nRow; i++) {
 	for (int j = 0; j < nCol; j++) {
-	  (*mat)[i][j] += c[0] * q00->getValue(i, j);
+	  mat[i][j] += c[0] * q00->getValue(i, j);
 	}
       }
     }
   }
 
-  void PrecalcZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
+  void PrecalcZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
   {
     if (firstCall) {
 #ifdef _OPENMP
@@ -452,7 +450,7 @@ namespace AMDiS {
     c[0] *= elInfo->getDet();
 
     for (int i = 0; i < nRow; i++)
-      (*vec)[i] += c[0] * q0->getValue(i);
+      vec[i] += c[0] * q0->getValue(i);
   }
 
 }
diff --git a/AMDiS/src/ZeroOrderAssembler.h b/AMDiS/src/ZeroOrderAssembler.h
index 2d89690ef37d3c64d625b15cc6260cc6c5670c3b..99bf4578c388c9fd03930707af07bad654ac3f38 100644
--- a/AMDiS/src/ZeroOrderAssembler.h
+++ b/AMDiS/src/ZeroOrderAssembler.h
@@ -18,8 +18,6 @@ namespace AMDiS {
   class ZeroOrderAssembler : public SubAssembler
   {
   public:
-    MEMORY_MANAGED(ZeroOrderAssembler);
-
     /** \brief
      * Creates and returns the ZeroOrderAssembler for Operator op and
      * the given assembler. If all terms are piecewise constant precalculated 
@@ -61,22 +59,20 @@ namespace AMDiS {
   class StandardZOA :  public ZeroOrderAssembler 
   {
   public:
-    MEMORY_MANAGED(StandardZOA);
-
     /// Constructor.
     StandardZOA(Operator *op, Assembler *assembler, Quadrature *quad);
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat);
+				ElementMatrix& mat);
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *elInfo, ElementVector *vec);
+    void calculateElementVector(const ElInfo *elInfo, ElementVector& vec);
   };
 
 
@@ -89,22 +85,20 @@ namespace AMDiS {
   class FastQuadZOA :  public ZeroOrderAssembler
   {
   public:
-    MEMORY_MANAGED(FastQuadZOA);
-
     /// Constructor.
     FastQuadZOA(Operator *op, Assembler *assembler, Quadrature *quad);
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat);
+				ElementMatrix& mat);
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *elInfo, ElementVector *vec);
+    void calculateElementVector(const ElInfo *elInfo, ElementVector& vec);
   };
 
 
@@ -117,19 +111,17 @@ namespace AMDiS {
   class PrecalcZOA : public ZeroOrderAssembler
   {
   public:
-    MEMORY_MANAGED(PrecalcZOA);
-
     /// Constructor.
     PrecalcZOA(Operator *op, Assembler *assembler, Quadrature *quad);
 
     /// Implements SubAssembler::calculateElementMatrix().
-    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat);
+    void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
 
     void calculateElementMatrix(const ElInfo *rowElInfo,
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
-				ElementMatrix *mat) 
+				ElementMatrix& mat) 
     {
       FUNCNAME("PrecalcZOA::calculateElementMatrix()");
       
@@ -137,7 +129,7 @@ namespace AMDiS {
     }
 
     /// Implements SubAssembler::calculateElementVector().
-    void calculateElementVector(const ElInfo *elInfo, ElementVector *vec);
+    void calculateElementVector(const ElInfo *elInfo, ElementVector& vec);
 
   protected:
     /// Integral of the product of psi and phi.