diff --git a/AMDiS/src/FiniteElemSpace.h b/AMDiS/src/FiniteElemSpace.h
index 5d662816a931f1bc50ed64ac2862d4876452b2a4..3e4ee0533de6b669f71488b197847a397cb22fd7 100644
--- a/AMDiS/src/FiniteElemSpace.h
+++ b/AMDiS/src/FiniteElemSpace.h
@@ -76,28 +76,28 @@ namespace AMDiS {
      */  
     inline std::string getName() const { 
       return name;
-    };
+    }
 
     /** \brief
      * Returns \ref admin
      */
     inline DOFAdmin* getAdmin() const { 
       return admin;
-    };
+    }
 
     /** \brief
      * Returns \ref basFcts
      */
     inline const BasisFunction* getBasisFcts() const { 
       return basFcts;
-    };
+    }
 
     /** \brief
      * Returns \ref mesh
      */
     inline Mesh* getMesh() const { 
       return mesh; 
-    };
+    }
   
   protected:
     /** \brief
@@ -107,7 +107,7 @@ namespace AMDiS {
     FiniteElemSpace(DOFAdmin* admin_,
 		    const BasisFunction* basisFcts, 
 		    Mesh* mesh,
-		    const std::string& name_="");
+		    const std::string& name_ = "");
 
   protected:
     /** \brief
diff --git a/AMDiS/src/SolutionDataStorage.h b/AMDiS/src/SolutionDataStorage.h
index 5d5ff0ea0d8840dc8b11e0c64440800806953be0..87c312d33976747c758b04d54817d4c0a8ca38aa 100644
--- a/AMDiS/src/SolutionDataStorage.h
+++ b/AMDiS/src/SolutionDataStorage.h
@@ -25,8 +25,14 @@ namespace AMDiS {
   class SolutionDataStorage 
   {
   public:
+    /** \brief
+     *
+     */
     SolutionDataStorage(std::string name);
 
+    /** \brief
+     *
+     */
     ~SolutionDataStorage();
 
     /** \brief
@@ -35,17 +41,60 @@ namespace AMDiS {
     void setFixFESpace(typename SolutionHelper<T>::type feSpace,
 		       bool cleanupFeSpace = false);
 
+    /** \brief
+     *
+     */
+    void push(T *solution,
+	      double timestamp,
+	      bool cleanupSolution = true);
+
+    /** \brief
+     *
+     */
     void push(T *solution,
 	      double timestamp,
-	      FiniteElemSpace *feSpace = NULL,
+	      typename SolutionHelper<T>::type feSpace,
 	      bool cleanupSolution = true,
 	      bool cleanupFeSpace = true);
 
-    bool pop(T *solution,
+    /** \brief
+     *
+     */
+    bool pop(T **solution,
+	     double *timestep);
+
+    bool pop(T **solution,
 	     double *timestep,
-	     FiniteElemSpace *feSpace = NULL);
+	     typename SolutionHelper<T>::type feSpace);
+
+    /** \brief
+     * Returns for a given solution number the corresponding fe Space. If the
+     * the fe Space is fixed, the fe Space for all solutions is stored at
+     * position 0.
+     */
+    typename SolutionHelper<T>::type getFeSpace(int i = 0) {
+      return feSpaces[i];
+    }
 
   protected:
+    /** \brief
+     * Auxiliary function for deleting either a single fe space, or a
+     * vector of fe spaces.
+     */
+    void deleteFeSpace(FiniteElemSpace* feSpace) {
+      DELETE feSpace;
+    }
+
+    void deleteFeSpace(std::vector<FiniteElemSpace*> feSpaces) {
+      for (int i = 0; i < static_cast<int>(feSpaces.size()); i++)
+	DELETE feSpaces[i];
+
+      feSpaces.empty();
+    }
+
+    /** \brief
+     * Deletes all pointers and empties all internal vectors.
+     */
     void cleanup();
 
     /** \brief
@@ -72,7 +121,7 @@ namespace AMDiS {
      * Here, all the solutions (i.e. either DOFVectors or SystemVectors)
      * are stored.
      */
-    std::vector<T> solutions;
+    std::vector<T*> solutions;
 
     /** \brief
      * Stores to every solution its FE Space. If \ref fixedFESpace is set
diff --git a/AMDiS/src/SolutionDataStorage.hh b/AMDiS/src/SolutionDataStorage.hh
index 7472667c52fa37f0106adab358d32ade74e3be3e..cbcb1709ffa6bc219669e97cc581e109f0103f28 100644
--- a/AMDiS/src/SolutionDataStorage.hh
+++ b/AMDiS/src/SolutionDataStorage.hh
@@ -44,9 +44,7 @@ namespace AMDiS {
   template<typename T>
   void SolutionDataStorage<T>::push(T *solution,
 				    double timestamp,
-				    FiniteElemSpace *feSpace,
-				    bool cleanupSolution,
-				    bool cleanupFeSpace)
+				    bool cleanupSolution)
   {
     // If pop was the last operation, cleanup and reset the data storage.
     if (poped) {
@@ -58,29 +56,35 @@ namespace AMDiS {
     timestamps.push_back(timestamp);
     cleanupSolutions.push_back(cleanupSolution);
 
+    lastPos++;
+  }
+
+  template<typename T>
+  void SolutionDataStorage<T>::push(T *solution,
+				    double timestamp,
+				    typename SolutionHelper<T>::type feSpace,
+				    bool cleanupSolution,
+				    bool cleanupFeSpace)
+  {
+    push(solution, timestamp, cleanupSolution, cleanupFeSpace);
+
     // Store fe space only, if we do not have a fixed fe space.
     if (!fixedFESpace) {
       feSpaces.push_back(feSpace);
       cleanupFeSpaces.push_back(cleanupFeSpace);
     }
-
-    lastPos++;
   }
 
   template<typename T>
-  bool SolutionDataStorage<T>::pop(T *solution,
-				   double *timestamp,
-				   FiniteElemSpace *feSpace)
+  bool SolutionDataStorage<T>::pop(T **solution,
+				   double *timestamp)
   {
     if (lastPos < 0) {
       return false;
     }
 
-    solution = solutions[lastPos];
+    *solution = solutions[lastPos];
     *timestamp = timestamps[lastPos];
-    if (!fixedFESpace) {
-      feSpace = feSpaces[lastPos];
-    }
 
     lastPos--;
     poped = true;
@@ -88,6 +92,24 @@ namespace AMDiS {
     return true;
   }
 
+  template<typename T>
+  bool SolutionDataStorage<T>::pop(T **solution,
+				   double *timestep,
+				   typename SolutionHelper<T>::type feSpace)
+  {
+    if (!pop(solution, timestep))
+      return false;
+
+    if (!fixedFESpace) {
+      feSpace = feSpaces[lastPos + 1]; // + 1, because lastPos was decremented in pop call above
+    } else {
+      feSpace = feSpaces[0];
+    }
+
+    return true;
+  }
+
+
   template<typename T>
   void SolutionDataStorage<T>::cleanup()
   {
@@ -99,7 +121,7 @@ namespace AMDiS {
 
     for (int i = 0; i < static_cast<int>(cleanupFeSpaces.size()); i++) {
       if (cleanupFeSpaces[i]) {
-	DELETE feSpaces[i];
+	deleteFeSpace(feSpaces[i]);
       }
     }
 
diff --git a/AMDiS/src/SystemVector.h b/AMDiS/src/SystemVector.h
index 9c372e5eb54d4a7272d1d3a7ee04db8c19000929..0117ba844a9e4c645a14a8298a47bb6a35793908 100644
--- a/AMDiS/src/SystemVector.h
+++ b/AMDiS/src/SystemVector.h
@@ -124,7 +124,7 @@ namespace AMDiS {
       
     {
       vectors.set(NULL);
-    };
+    }
 
     /** \brief
      * Copy Constructor.
@@ -136,11 +136,19 @@ namespace AMDiS {
       
     {
       for (int i = 0; i < vectors.getSize(); i++) {
-	vectors[i] = new DOFVector<double>(*rhs.vectors[i]);
+	vectors[i] = NEW DOFVector<double>(*rhs.vectors[i]);
       }
-    };
+    }
+
+    virtual ~SystemVector() 
+    {}
 
-    virtual ~SystemVector() {};
+    void createNewDOFVectors(std::string name)
+    {     
+      for (int i = 0; i < vectors.getSize(); i++) {
+	vectors[i] = NEW DOFVector<double>(feSpace[i], "tmp");
+      }
+    }
 
     /** \brief
      * Sets \ref vectors[index] = vec.
@@ -148,7 +156,7 @@ namespace AMDiS {
     inline void setDOFVector(int index, DOFVector<double> *vec) {
       TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
       vectors[index] = vec;
-    };
+    }
 
     /** \brief
      * Returns \ref vectors[index].
@@ -156,7 +164,7 @@ namespace AMDiS {
     inline DOFVector<double> *getDOFVector(int index) {      
       TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
       return vectors[index];
-    };
+    }
 
     /** \brief
      * Returns \ref vectors[index].
@@ -164,7 +172,7 @@ namespace AMDiS {
     inline const DOFVector<double> *getDOFVector(int index) const {
       TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
       return vectors[index];
-    };
+    }
 
     /** \brief
      * Returns sum of used vector sizes.
@@ -176,25 +184,25 @@ namespace AMDiS {
 	totalSize += vectors[i]->getUsedSize();
       }
       return totalSize;
-    };
+    }
 
     /** \brief
      * Returns number of contained vectors.
      */
     inline int getNumVectors() const { 
       return vectors.getSize(); 
-    };
+    }
 
     inline int getSize() const {
       return vectors.getSize();
-    };
+    }
 
     /** \brief
      * Returns \ref feSpace.
      */
     inline FiniteElemSpace *getFESpace(int i) const { 
       return feSpace[i]; 
-    };
+    }
 
     /** \brief
      * Here the system vector is interpreted as one large vector. The given
@@ -211,7 +219,7 @@ namespace AMDiS {
       }
 
       return (*(vectors[vectorIndex]))[localIndex];
-    };
+    }
 
     /** \brief
      * For const access.
@@ -225,7 +233,7 @@ namespace AMDiS {
       }
 
       return (*(vectors[vectorIndex]))[localIndex];
-    };
+    }
 
     /** \brief
      * Sets all entries in all vectors to value.
@@ -235,7 +243,7 @@ namespace AMDiS {
       for (int i = 0; i < size; i++) {
 	vectors[i]->set(value);
       }
-    };
+    }
 
     /** \brief
      * Sets all entries in all vectors to value.
@@ -246,7 +254,7 @@ namespace AMDiS {
 	(*(vectors[i])) = value;
       }
       return *this;
-    };
+    }
 
     /** \brief
      * Assignement operator.
@@ -258,7 +266,7 @@ namespace AMDiS {
 	(*(vectors[i])) = (*(rhs.getDOFVector(i)));
       }
       return *this;
-    };
+    }
 
     void serialize(std::ostream &out) {
       int size = vectors.getSize();
@@ -266,7 +274,7 @@ namespace AMDiS {
       for (int i = 0; i < size; i++) {
 	vectors[i]->serialize(out);
       }
-    };
+    }
 
     void deserialize(std::istream &in) {
       int size, oldSize = vectors.getSize();
@@ -278,7 +286,7 @@ namespace AMDiS {
       for (int i = 0; i < size; i++) {
 	vectors[i]->deserialize(in);
       }  
-    };
+    }
 
     void copy(const SystemVector& rhs) {
       int size = vectors.getSize();
@@ -286,21 +294,27 @@ namespace AMDiS {
       for (int i = 0; i < size; i++) {
 	vectors[i]->copy(*(const_cast<SystemVector&>(rhs).getDOFVector(i)));
       }
-    };
+    }
 
     void interpol(std::vector<AbstractFunction<double, WorldVector<double> >*> *f) {
       int size = vectors.getSize();
       for (int i = 0; i < size; i++) {
 	vectors[i]->interpol((*f)[i]);
       }
-    };
+    }
+
+    void interpol(SystemVector *v, double factor) {
+      for (int i = 0; i < v->getSize(); i++) {
+	vectors[i]->interpol(v->getDOFVector(i), factor);
+      }
+    }
 
     void print() {
       int size = vectors.getSize();
       for (int i = 0; i < size; i++) {
 	vectors[i]->print();
       }
-    };
+    }
 
 
 
@@ -330,7 +344,7 @@ namespace AMDiS {
       *(x.getDOFVector(i)) *= d;
     }
     return x;
-  };
+  }
 
   /** \brief
    * scalar product
@@ -342,9 +356,9 @@ namespace AMDiS {
     int size = x.getNumVectors();
     for (int i = 0; i < size; i++) {
       result += (*x.getDOFVector(i)) * (*y.getDOFVector(i));
-    };
+    }
     return result;
-  };
+  }
 
   /** \brief
    * addition of two system vectors
@@ -359,7 +373,7 @@ namespace AMDiS {
       (*(x.getDOFVector(i))) += (*(y.getDOFVector(i)));
     }
     return x;
-  };
+  }
 
   /**
    * subtraction of two system vectors.
@@ -374,7 +388,7 @@ namespace AMDiS {
       (*(x.getDOFVector(i))) -= (*(y.getDOFVector(i)));
     }
     return x;
-  };
+  }
 
   /** \brief
    * multiplication with a scalar
@@ -386,7 +400,7 @@ namespace AMDiS {
       (*(result.getDOFVector(i))) *= d;
     }
     return result;
-  };
+  }
 
   /** \brief
    * multiplication with a scalar
@@ -398,7 +412,7 @@ namespace AMDiS {
       (*(result.getDOFVector(i))) *= d;
     }
     return result;
-  };
+  }
 
   /** \brief
    * addition of two system vectors
@@ -414,21 +428,21 @@ namespace AMDiS {
       (*(result.getDOFVector(i))) += (*(y.getDOFVector(i)));
     }
     return result;
-  };
+  }
 
   /** \brief
    * Calls SystemVector::set(). Used for solving.
    */
   inline void set(SystemVector& x, double value) {
     x.set(value);
-  }; 
+  } 
 
   /** \brief
    * Calls SystemVector::set(). Used for solving.
    */
   inline void setValue(SystemVector& x, double value) {
     x.set(value);
-  }; 
+  }
 
   /** \brief
    * Norm of system vector.
@@ -440,7 +454,7 @@ namespace AMDiS {
       result += x->getDOFVector(i)->squareNrm2();
     }
     return sqrt(result);
-  };
+  }
 
   /** \brief
    * L2 norm of system vector.
@@ -452,7 +466,7 @@ namespace AMDiS {
       result += x->getDOFVector(i)->L2NormSquare();
     }
     return sqrt(result);
-  };
+  }
 
   /** \brief
    * H1 norm of system vector.
@@ -464,7 +478,7 @@ namespace AMDiS {
       result += x->getDOFVector(i)->H1NormSquare();
     }
     return sqrt(result);
-  };
+  }
 
   inline void mv(Matrix<DOFMatrix*> &matrix,
 		 const SystemVector &x,
@@ -495,7 +509,7 @@ namespace AMDiS {
 	}
       }
     }
-  };
+  }
 
   /** \brief
    * y = a*x + y;
@@ -514,7 +528,7 @@ namespace AMDiS {
     for (i = 0; i < size; i++) {
       axpy(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
     }
-  };
+  }
 
   /** \brief
    * y = x + a*y
@@ -535,11 +549,11 @@ namespace AMDiS {
    */
   inline int size(SystemVector* vec) {
     return vec->getUsedSize();
-  };
+  }
 
   inline void print(SystemVector* vec) {
     vec->print();
-  };
+  }
 
 }