diff --git a/AMDiS/src/Config.h b/AMDiS/src/Config.h
index a9f7415bf27bef60f5597a16af0fe483a98f7480..5e52b3bc5417810db44c8435cb11edf0335dbe4b 100644
--- a/AMDiS/src/Config.h
+++ b/AMDiS/src/Config.h
@@ -1,3 +1,27 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: 
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+
+
+/** \file Config.h */
+
 #pragma once
 
 /** \brief current AMDiS version */
@@ -35,4 +59,4 @@
   error: not supported compiler
 #endif
 
-#include "config/Config_defaults.h"
\ No newline at end of file
+#include "config/Config_defaults.h"
diff --git a/AMDiS/src/CreatorInterface.h b/AMDiS/src/CreatorInterface.h
index 00e008ed596653b54de676593aef8fc0a5d01744..dd6e436ad9390f4dd97135b9013806f8b52e2e4b 100644
--- a/AMDiS/src/CreatorInterface.h
+++ b/AMDiS/src/CreatorInterface.h
@@ -84,7 +84,7 @@ namespace AMDiS {
     virtual ~NullCreator() {}
 
     /// Implementation of \ref CreatorInterface::isNullCreator()
-    bool isNullCreator() override
+    virtual bool isNullCreator() override
     { 
       return true; 
     }
diff --git a/AMDiS/src/config/Config_clang.h b/AMDiS/src/config/Config_clang.h
index 1daa7f50061a6d7dd48a77e87d3f95d67010cefd..35b586e07741e6bd6c98e1373cc3397e154d9bee 100644
--- a/AMDiS/src/config/Config_clang.h
+++ b/AMDiS/src/config/Config_clang.h
@@ -1,7 +1,33 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: 
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+
+
+/** \file Config_clang.h */
+
 #pragma once
 
 #define CLANG_VERSION (__clang_major__ * 10000 + __clang_minor__ * 100 + __clang_patchlevel__)
-#define COMPILER_VERSION "Clang: CLANG_VERSION"
+
+#define COMPILER_NAME "clang"
+#define COMPILER_VERSION CLANG_VERSION
 
 // alignement specification
 // ------------------------
@@ -21,6 +47,7 @@ typedef size_t aligned_size_t   __attribute__ ((aligned(CACHE_LINE)));
 // --------------
 #if __cplusplus > 199711L
 
+// __has_feature(cxx_rvalue_references)
 #if CLANG_VERSION >= 20900
   #define HAS_VARIADIC_TEMPLATES 1
 #endif
@@ -57,4 +84,4 @@ typedef size_t aligned_size_t   __attribute__ ((aligned(CACHE_LINE)));
   #define HAS_TYPED_ENUMS 1
 #endif
 
-#endif
\ No newline at end of file
+#endif
diff --git a/AMDiS/src/config/Config_defaults.h b/AMDiS/src/config/Config_defaults.h
index 8341f43ad54b872b71b7826b56e1b4886b235a7b..64679fecf6825c81684864d55e6f450ac6a3f44f 100644
--- a/AMDiS/src/config/Config_defaults.h
+++ b/AMDiS/src/config/Config_defaults.h
@@ -1,3 +1,27 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: 
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+
+
+/** \file Config_defaults.h */
+
 #pragma once
 
 #ifndef COMPILER_NAME
@@ -66,4 +90,4 @@ typedef size_t aligned_size_t;
 
 #ifndef HAS_TYPED_ENUMS
   #define HAS_TYPED_ENUMS 0
-#endif
\ No newline at end of file
+#endif
diff --git a/AMDiS/src/config/Config_gcc.h b/AMDiS/src/config/Config_gcc.h
index 2247ed3beb5ba1dcad252e7558fc700bd5908539..a5aeb35efd62c714883d0c410a1dbe96bdcf25a3 100644
--- a/AMDiS/src/config/Config_gcc.h
+++ b/AMDiS/src/config/Config_gcc.h
@@ -1,8 +1,32 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: 
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+
+
+/** \file Config_gcc.h */
+
 #pragma once
 
 #define GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
 
-#define COMPIER_NAME "gcc"
+#define COMPILER_NAME "gcc"
 #define COMPILER_VERSION GCC_VERSION
 
 // alignement specification
@@ -31,7 +55,7 @@ typedef size_t aligned_size_t   __attribute__ ((aligned(CACHE_LINE)));
   #define HAS_ALIAS_TEMPLATES 1
 #endif
 
-#if __cpp_decltype >= 200707
+#if GCC_VERSION >= 40300
   #define HAS_DECLTYPE 1
 #endif
 
@@ -59,4 +83,4 @@ typedef size_t aligned_size_t   __attribute__ ((aligned(CACHE_LINE)));
   #define HAS_TYPED_ENUMS 1
 #endif
 
-#endif
\ No newline at end of file
+#endif
diff --git a/AMDiS/src/config/Config_intel.h b/AMDiS/src/config/Config_intel.h
index 4053bcb9017953ae8d048565c916221c3ff1e533..5b40615065f38df734b6ea2ee8724f891aa707de 100644
--- a/AMDiS/src/config/Config_intel.h
+++ b/AMDiS/src/config/Config_intel.h
@@ -1,9 +1,33 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: 
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+
+
+/** \file Config_intel.h */
+
 #pragma once
 
 // MMmm (M...major, m...minor)
 #define INTEL_VERSION __INTEL_COMPILER
 
-#define COMPIER_NAME "icc"
+#define COMPILER_NAME "icc"
 #define COMPILER_VERSION INTEL_VERSION
 
 // alignement specification
@@ -60,4 +84,4 @@ typedef __declspec(align(CACHE_LINE)) size_t aligned_size_t;
   #define HAS_TYPED_ENUMS 1
 #endif
 
-#endif
\ No newline at end of file
+#endif
diff --git a/AMDiS/src/config/Config_msc.h b/AMDiS/src/config/Config_msc.h
index 14edbf107821a100d63fd6e16ff39e153d776c00..9d2d41aba14d757eeb551f6b9bf9c03fb97100ee 100644
--- a/AMDiS/src/config/Config_msc.h
+++ b/AMDiS/src/config/Config_msc.h
@@ -1,3 +1,27 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: 
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+
+
+/** \file Config_msc.h */
+
 #pragma once
 
 #define MSC_VERSION _MSV_VER
@@ -6,7 +30,7 @@
 // MSC_VERSION == 1600 (Visual Studio 2010)
 // MSC_VERSION == 1500 (Visual Studio 2008)
 
-#define COMPIER_NAME "msc"
+#define COMPILER_NAME "msc"
 #define COMPILER_VERSION MSC_VERSION
 
 // alignement specification
@@ -63,4 +87,4 @@ typedef __declspec(align(CACHE_LINE)) size_t aligned_size_t;
   #define HAS_TYPED_ENUMS 1
 #endif
 
-#endif
\ No newline at end of file
+#endif
diff --git a/AMDiS/src/config/Config_pgi.h b/AMDiS/src/config/Config_pgi.h
index 377a985072d75cb9aa8cf824d528dd688ec00751..2b8c36806290213c513bbe85550810849dbd9e78 100644
--- a/AMDiS/src/config/Config_pgi.h
+++ b/AMDiS/src/config/Config_pgi.h
@@ -1,6 +1,30 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: 
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+
+
+/** \file Config_pgi.h */
+
 #pragma once
 
 #define PGI_VERSION (__PGIC__ * 10000 + __PGIC_MINOR * 100 + __PGIC_PATCHLEVEL__)
 
-#define COMPIER_NAME "pgi"
-#define COMPILER_VERSION PGI_VERSION
\ No newline at end of file
+#define COMPILER_NAME "pgi"
+#define COMPILER_VERSION PGI_VERSION
diff --git a/AMDiS/src/expressions/functorN_expr.hpp b/AMDiS/src/expressions/functorN_expr.hpp
index f333949f6c599cfb8f211e4f86018d4a315325f9..8c8d7debf2c81c4e87961cc0277203822ba1702b 100644
--- a/AMDiS/src/expressions/functorN_expr.hpp
+++ b/AMDiS/src/expressions/functorN_expr.hpp
@@ -66,7 +66,7 @@ namespace AMDiS
     };
     
     template<typename F>
-    struct functor_degree<F, typename std::enable_if<boost::is_base_of<FunctorBase, F>::value>::type >
+    struct functor_degree<F, typename enable_if< boost::is_base_of<FunctorBase, F> >::type >
     {
       template<typename... Int>
       static int eval(F f, Int... d) { return f.getDegree(d...); }
@@ -217,9 +217,9 @@ namespace AMDiS
   {
     // result of generator-functions (used for enable_if constructs)
     template<typename F, typename... Terms>
-    struct FunctionN : boost::enable_if
+    struct FunctionN : enable_if
       <
-	typename boost::mpl::and_< typename traits::is_valid_arg<Terms>::type... >::type,
+	typename and_< typename traits::is_valid_arg<Terms>::type... >::type,
 	expressions::FunctionN< F, typename traits::to_expr<Terms>::type...> 
       > {};
       
diff --git a/AMDiS/src/expressions/value_expr.hpp b/AMDiS/src/expressions/value_expr.hpp
index 4f6de8bde3379e157904b2af5158712bff4843f4..aa9277b33eb526f14cf473ab7a875b03c0e97b64 100644
--- a/AMDiS/src/expressions/value_expr.hpp
+++ b/AMDiS/src/expressions/value_expr.hpp
@@ -187,7 +187,7 @@ namespace AMDiS
       
       Eye() {
 	M = 0;
-	for (size_t i = 0; i < Global::getGeo(WORLD); ++i)
+	for (size_t i = 0; i < (size_t)Global::getGeo(WORLD); ++i)
 	  M[i][i] = 1;
       }
 		    
diff --git a/AMDiS/src/expressions/vec_functors.hpp b/AMDiS/src/expressions/vec_functors.hpp
index d67d870bc5e234c8adc2b5c8eb681bcdfe70f5e8..904e362307a341b7d5b2c0232b682d863230c9d0 100644
--- a/AMDiS/src/expressions/vec_functors.hpp
+++ b/AMDiS/src/expressions/vec_functors.hpp
@@ -134,14 +134,15 @@ namespace AMDiS
   inline typename result_of::UnaryExpr<Functor, Term>::type
   unary_expr(const Term& t) 
   { 
-    return function_< Functor< typename Term::value_type > >(t); 
+    typedef Functor< typename Term::value_type > F;
+    return function_(F(), t); 
   }
   
-  template<class Functor, typename Term>
-  inline typename result_of::UnaryExprFull<Functor, Term>::type
+  template<class F, typename Term>
+  inline typename result_of::UnaryExprFull<F, Term>::type
   unary_expr_full(const Term& t) 
   { 
-    return function_< Functor >(t); 
+    return function_(F(), t); 
   }
   
     
@@ -152,9 +153,8 @@ namespace AMDiS
   { 
     typedef typename traits::to_expr<Term1>::to Expr1;
     typedef typename traits::to_expr<Term2>::to Expr2;
-    return function_
-      < Functor< typename Expr1::type::value_type, typename Expr2::type::value_type > >
-      (Expr1::get(t1), Expr2::get(t2)); 
+    typedef Functor< typename Expr1::type::value_type, typename Expr2::type::value_type > F;
+    return function_(F(), Expr1::get(t1), Expr2::get(t2)); 
   }
     
 
diff --git a/AMDiS/src/io/Arh2Writer.cc b/AMDiS/src/io/Arh2Writer.cc
index fd39d5b68bc2bef5a7d63c3a2c91b12c6c94acad..b125ba4f60940354f2f56851e06ab805e82e58a9 100644
--- a/AMDiS/src/io/Arh2Writer.cc
+++ b/AMDiS/src/io/Arh2Writer.cc
@@ -107,4 +107,4 @@ namespace AMDiS { namespace io {
 #endif
 
   }
-} }
\ No newline at end of file
+} }
diff --git a/AMDiS/src/io/detail/VtkWriter.cc b/AMDiS/src/io/detail/VtkWriter.cc
index 50f71973d5018cb1075db517e3cb8f6632f994cc..1c0494312dab0cfabfa802e176fa74a212b36c66 100644
--- a/AMDiS/src/io/detail/VtkWriter.cc
+++ b/AMDiS/src/io/detail/VtkWriter.cc
@@ -204,7 +204,7 @@ namespace AMDiS  { namespace io {
         
 	void writeParallelFile(string name, int nRanks,
 			      string fnPrefix, string fnPostfix,
-			      vector<string> &componentNames,
+			      const vector<string> &componentNames,
 			      ::AMDiS::io::VtkWriter::Vtuformat format,
 			      bool highPrecision,
 			      bool writeAsVector
diff --git a/AMDiS/src/solver/BlockMTLMatrix.h b/AMDiS/src/solver/BlockMTLMatrix.h
index c25f58666a92d0e001ef5de8d6d179eee63e4b2b..48dad2a226cc1b2af525db0b3648a2dc62937b9b 100644
--- a/AMDiS/src/solver/BlockMTLMatrix.h
+++ b/AMDiS/src/solver/BlockMTLMatrix.h
@@ -176,4 +176,4 @@ namespace mtl
   
 } // end namespace mtl
 
-#endif // AMDIS_BLOCK_MTL_MATRIX_H
\ No newline at end of file
+#endif // AMDIS_BLOCK_MTL_MATRIX_H
diff --git a/AMDiS/src/solver/ITL_Preconditioner.h b/AMDiS/src/solver/ITL_Preconditioner.h
index 8c798fb2a9eb65ffff0aec89b98192854be4ebca..a22890d069434c00b16dc28da8b668dbe1718664 100644
--- a/AMDiS/src/solver/ITL_Preconditioner.h
+++ b/AMDiS/src/solver/ITL_Preconditioner.h
@@ -103,7 +103,7 @@ namespace AMDiS {
     }
     
     /// Implementation of \ref ITL_PreconditionerBase::init()
-    void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override
+    virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override
     {
       if (precon)
 	delete precon;
@@ -111,7 +111,7 @@ namespace AMDiS {
     }
     
     /// Implementation of \ref PreconditionerInterface::exit()
-    void exit() override
+    virtual void exit() override
     {
       if (precon) {
 	delete precon;
@@ -120,14 +120,14 @@ namespace AMDiS {
     }
     
     /// Implementation of \ref ITL_PreconditionerBase::solve()
-    void solve(const VectorType& vin, VectorType& vout) const override
+    virtual void solve(const VectorType& vin, VectorType& vout) const override
     {
       assert(precon != NULL);
       precon->solve(vin, vout);
     }
     
     /// Implementation of \ref ITL_PreconditionerBase::adjoint_solve()
-    void adjoint_solve(const VectorType& vin, VectorType& vout) const override
+    virtual void adjoint_solve(const VectorType& vin, VectorType& vout) const override
     {
       assert(precon != NULL);
       precon->adjoint_solve(vin, vout);
diff --git a/AMDiS/src/solver/ITL_Runner.h b/AMDiS/src/solver/ITL_Runner.h
index 00c8f550f57ff515d6f834536bfc123be19a752e..84669bca75eb6a9c81c8d08e87f002ce536e9731 100644
--- a/AMDiS/src/solver/ITL_Runner.h
+++ b/AMDiS/src/solver/ITL_Runner.h
@@ -83,7 +83,7 @@ namespace AMDiS {
     
     
     /// Implementation of \ref RunnerBase::init()
-    void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override
+    virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override
     {
       preconPair.l->init(A, fullMatrix);
       preconPair.r->init(A, fullMatrix);
@@ -91,7 +91,7 @@ namespace AMDiS {
     
     
     /// Implementation of \ref RunnerBase::solve()
-    int solve(const MatrixType& A , VectorType& x, const VectorType& b) override
+    virtual int solve(const MatrixType& A , VectorType& x, const VectorType& b) override
     { FUNCNAME("ITL_Runner::solve()");
     
       TEST_EXIT(preconPair.l != NULL)("there is no left preconditioner\n");
@@ -133,7 +133,7 @@ namespace AMDiS {
 
     
     /// Implementation of \ref RunnerBase::adjoint_solve()
-    int adjoint_solve(const MatrixType& A , VectorType& x, const VectorType& b) override
+    virtual int adjoint_solve(const MatrixType& A , VectorType& x, const VectorType& b) override
     { FUNCNAME("ITL_Runner::adjoint_solve()");
     
       TEST_EXIT(preconPair.l != NULL)("there is no left preconditioner\n");
@@ -170,7 +170,7 @@ namespace AMDiS {
       
       
     /// Implementation of \ref RunnerInterface::exit()
-    void exit() override
+    virtual void exit() override
     {
       preconPair.l->exit();
       preconPair.r->exit();
@@ -178,14 +178,14 @@ namespace AMDiS {
     
     
     /// Implementation of \ref RunnerInterface::getLeftPrecon()
-    PreconditionerInterface* getLeftPrecon() override
+    virtual PreconditionerInterface* getLeftPrecon() override
     {
       return preconPair.l;
     }    
     
     
     /// Implementation of \ref RunnerInterface::getRightPrecon()
-    PreconditionerInterface* getRightPrecon() override
+    virtual PreconditionerInterface* getRightPrecon() override
     {
       return preconPair.r;
     }
diff --git a/AMDiS/src/solver/KrylovPreconditioner.h b/AMDiS/src/solver/KrylovPreconditioner.h
index 2b314aecce555829d05dfe6db6d577275317c814..c9dcf5f56bc806baa83212372d3b356b505c44f6 100644
--- a/AMDiS/src/solver/KrylovPreconditioner.h
+++ b/AMDiS/src/solver/KrylovPreconditioner.h
@@ -108,27 +108,27 @@ namespace AMDiS {
     }
     
     /// Implementation of \ref ITL_PreconditionerBase::init()
-    void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix_) override
+    virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix_) override
     {
       fullMatrix = &fullMatrix_;
       runner->init(A, fullMatrix_);
     }
     
     /// Implementation of \ref PreconditionerInterface::init()
-    void exit() override
+    virtual void exit() override
     {
       runner->exit();
     }
     
     /// Implementation of \ref ITL_PreconditionerBase::solve()
-    void solve(const VectorType& b, VectorType& x) const override
+    virtual void solve(const VectorType& b, VectorType& x) const override
     {
       initVector(x);
       runner->solve(*fullMatrix, x, b);
     }
     
     /// Implementation of \ref ITL_PreconditionerBase::adjoint_solve()
-    void adjoint_solve(const VectorType& b, VectorType& x) const override
+    virtual void adjoint_solve(const VectorType& b, VectorType& x) const override
     {
       initVector(x);
       runner->adjoint_solve(*fullMatrix, x, b);
diff --git a/AMDiS/src/solver/LinearSolver.h b/AMDiS/src/solver/LinearSolver.h
index 5dece3930ed82fec5ec328970d181fe588ab53b8..a8cf7cbf18e986a62c23e918da6d10d671bd82ae 100644
--- a/AMDiS/src/solver/LinearSolver.h
+++ b/AMDiS/src/solver/LinearSolver.h
@@ -173,11 +173,11 @@ namespace AMDiS {
   protected:  
 
     /// Implementation of \ref LinearSolverInterface::solveLinearSystem()
-    int solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
-			  SystemVector& x,
-			  SystemVector& b,
-			  bool createMatrixData,
-			  bool storeMatrixData) override
+    virtual int solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+				  SystemVector& x,
+				  SystemVector& b,
+				  bool createMatrixData,
+				  bool storeMatrixData) override
     {    
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
       MPI::COMM_WORLD.Barrier();
diff --git a/AMDiS/src/solver/UmfPackSolver.h b/AMDiS/src/solver/UmfPackSolver.h
index f24538f880da1dbd6a165a0d38c22677fd189c38..8c51c8ebd7ca0f1417a610a340f1b66391fc856b 100644
--- a/AMDiS/src/solver/UmfPackSolver.h
+++ b/AMDiS/src/solver/UmfPackSolver.h
@@ -52,7 +52,7 @@ namespace AMDiS {
     }
 
     /// Implementation of \ref RunnerBase::init()
-    void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override
+    virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override
     {
       if (solver != NULL) {
 	delete solver;
@@ -67,7 +67,7 @@ namespace AMDiS {
     }
 
     /// Implementation of \ref RunnerBase::solve()
-    int solve(const MatrixType& A, VectorType& x, const VectorType& b) override
+    virtual int solve(const MatrixType& A, VectorType& x, const VectorType& b) override
     { 
       FUNCNAME("Umfpack_Runner::solve()");
       TEST_EXIT(solver != NULL)("The umfpack solver was not initialized\n");
@@ -89,7 +89,7 @@ namespace AMDiS {
     }
 
     /// Implementation of \ref RunnerInterface::adjoint_solve()
-    void exit() override {}
+    virtual void exit() override {}
 
     ~Umfpack_Runner() 
     {
diff --git a/extensions/ExtendedProblemStat.h b/extensions/ExtendedProblemStat.h
index fb58f1d516eaab09f01a06dfffe72000feba5a49..1e42e1aee89bfabbe534d01e498b5855893f6aa8 100644
--- a/extensions/ExtendedProblemStat.h
+++ b/extensions/ExtendedProblemStat.h
@@ -156,7 +156,9 @@ public:
       if (applyDirichletBC(singularDirichletBC[k], asmMatrix, asmVector))
 	++num_dbc;
     }
-    MSG_DBG("DBC applied at %d DOFs\n", num_dbc);
+    #if DEBUG
+    MSG("DBC applied at %d DOFs\n", num_dbc);
+    #endif
     
     // update solverMatrix
     if (asmMatrix && (singularDirichletBC.size() > 0 || manualPeriodicBC.size() > 0)) {
@@ -199,7 +201,6 @@ public:
   void addDirichletBC(BoundaryType type, int row, int col,
 		      ValueContainer *values)
   {
-  MSG("addDirichletBC(%d, %d)\n", row, col);
     BoundaryTypeContainer *bound = new BoundaryTypeContainer(type);
     DirichletBcDataList.push_back(
       new DirichletBcData<BoundaryTypeContainer, ValueContainer>(
diff --git a/extensions/base_problems/BaseProblem.h b/extensions/base_problems/BaseProblem.h
index 30df9b87c05075cfd9d8bdc8192f791f68533d81..67de4cec13ed78593292117d1152f390af3ef4f5 100644
--- a/extensions/base_problems/BaseProblem.h
+++ b/extensions/base_problems/BaseProblem.h
@@ -89,15 +89,15 @@ public:
   }
 
   /// This method is called before \ref beginIteration, \ref oneIteration and \ref endIteration.
-  virtual void initTimestep(AdaptInfo *adaptInfo) 
+  virtual void initTimestep(AdaptInfo *adaptInfo) override 
   {}
 
   /// calls \ref writeFiles
-  virtual void closeTimestep(AdaptInfo *adaptInfo);
+  virtual void closeTimestep(AdaptInfo *adaptInfo) override;
 
-  virtual void beginIteration(AdaptInfo *adaptInfo);
-  virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION);
-  virtual void endIteration(AdaptInfo *adaptInfo);
+  virtual void beginIteration(AdaptInfo *adaptInfo) override;
+  virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION) override;
+  virtual void endIteration(AdaptInfo *adaptInfo) override;
 
   /// Calls writeFiles of the problem
   virtual void writeFiles(AdaptInfo *adaptInfo, bool force)
diff --git a/extensions/base_problems/CahnHilliard.h b/extensions/base_problems/CahnHilliard.h
index a178877a83da36c69c563541d12870bd00b619e4..b9889e763d819350e35f588d9cf65f317f488823 100644
--- a/extensions/base_problems/CahnHilliard.h
+++ b/extensions/base_problems/CahnHilliard.h
@@ -42,15 +42,15 @@ namespace detail {
     CahnHilliard(const std::string &name_);
     ~CahnHilliard() {};
 
-    void solveInitialProblem(AMDiS::AdaptInfo *adaptInfo) override;
+    virtual void solveInitialProblem(AMDiS::AdaptInfo *adaptInfo) override;
 
     double getEpsilon() { return eps; }
     int getDoubleWellType() { return doubleWell; }
 
-    void fillOperators() override;
-    void fillBoundaryConditions() override {}
+    virtual void fillOperators() override;
+    virtual void fillBoundaryConditions() override {}
     
-    void finalizeData() override;
+    virtual void finalizeData() override;
 
   protected: // protected variables
 
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_RB.cc b/extensions/base_problems/CahnHilliardNavierStokes_RB.cc
index c6a0bf43f5f911e24483076eac3ad75b9f4fb7f6..9306d5c91a04c8e80214fb3d5b2da62e575900e5 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes_RB.cc
+++ b/extensions/base_problems/CahnHilliardNavierStokes_RB.cc
@@ -24,9 +24,18 @@ namespace AMDiS { namespace base_problems {
 
 CahnHilliardNavierStokes_RB::CahnHilliardNavierStokes_RB(const std::string &name_) :
   super(name_),
+  reinit(NULL),
+  velocity(NULL),
   useMobility(false),
   useConservationForm(false),
+  dim(2),
+  laplaceType(0),
+  nonLinTerm(3),
   doubleWell(1),
+  oldTimestep(0.0),
+  viscosity(1.0),
+  density(1.0),
+  c(0.0),
   gamma(1.0),
   eps(0.1),
   minusEps(-0.1),
@@ -34,16 +43,8 @@ CahnHilliardNavierStokes_RB::CahnHilliardNavierStokes_RB(const std::string &name
   minusEpsInv(-10.0),
   epsSqr(0.01),
   minusEpsSqr(-0.01),
-  laplaceType(0),
-  nonLinTerm(3),
-  oldTimestep(0.0),
-  viscosity(1.0),
-  density(1.0),
-  c(0.0),
   sigma(1.0),
   surfaceTension(1.0),
-  reinit(NULL),
-  velocity(NULL),
   fileWriter(NULL)
 {
   // parameters for CH
@@ -91,8 +92,7 @@ CahnHilliardNavierStokes_RB::CahnHilliardNavierStokes_RB(const std::string &name
 
 
 CahnHilliardNavierStokes_RB::~CahnHilliardNavierStokes_RB() 
-{ FUNCNAME("CahnHilliardNavierStokes_RB::~CahnHilliardNavierStokes_RB()");
-
+{ 
   if (reinit != NULL) {
     delete reinit;
     reinit = NULL;
@@ -107,8 +107,7 @@ CahnHilliardNavierStokes_RB::~CahnHilliardNavierStokes_RB()
 
 
 void CahnHilliardNavierStokes_RB::initData()
-{ FUNCNAME("CahnHilliardNavierStokes_RB::initData()");
-
+{ 
   dim = getMesh()->getDim();
   // create instance redistancing class
   reinit = new HL_SignedDistTraverse("reinit", dim);
@@ -122,8 +121,7 @@ void CahnHilliardNavierStokes_RB::initData()
 
 
 void CahnHilliardNavierStokes_RB::transferInitialSolution(AdaptInfo *adaptInfo)
-{ FUNCNAME("CahnHilliardNavierStokes_RB::transferInitialSolution()");
-
+{ 
   calcVelocity();
   
   for (int i = 0; i < 2+dow; i++)
@@ -367,8 +365,7 @@ void CahnHilliardNavierStokes_RB::fillOperators()
 
 
 void CahnHilliardNavierStokes_RB::addLaplaceTerm(int i)
-{ FUNCNAME("CahnHilliardNavierStokes_RB::addLaplaceTerm()");
-
+{
     /// < nu*[grad(u)+grad(u)^t] , grad(psi) >
     if (laplaceType == 1) {
       for (unsigned j = 0; j < dow; ++j) {
@@ -389,8 +386,7 @@ void CahnHilliardNavierStokes_RB::addLaplaceTerm(int i)
 }
 
 void CahnHilliardNavierStokes_RB::closeTimestep(AdaptInfo *adaptInfo)
-{ FUNCNAME("CahnHilliardNavierStokes_RB::closeTimestep()");
-
+{
   calcVelocity();
   fileWriter->writeFiles(adaptInfo, false);
   writeFiles(adaptInfo, false);
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_RB.h b/extensions/base_problems/CahnHilliardNavierStokes_RB.h
index a524247a2daf68f8cda0b824d1c13bdc6d57924d..0eee04e790a9bae454356ced42f691f2046fdb12 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes_RB.h
+++ b/extensions/base_problems/CahnHilliardNavierStokes_RB.h
@@ -52,11 +52,7 @@ public: // public methods
   virtual void fillBoundaryConditions() {}
   virtual void addLaplaceTerm(int i);
 
-protected: // protected variables
-
-  HL_SignedDistTraverse *reinit;
-
-  DOFVector<WorldVector<double> >* velocity;
+protected: // methods
 
   void calcVelocity()
   {
@@ -67,6 +63,12 @@ protected: // protected variables
     else if (dow == 3)
       transformDOF(prob->getSolution()->getDOFVector(2), prob->getSolution()->getDOFVector(3), prob->getSolution()->getDOFVector(4), velocity, new AMDiS::Vec3WorldVec<double>);
   }
+
+protected: // protected variables
+
+  HL_SignedDistTraverse *reinit;
+
+  DOFVector<WorldVector<double> >* velocity;
   
   bool useMobility;
   bool useConservationForm;
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.cc b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.cc
index b3a5959496aae52ec6b12f93d20e86717cae03a9..01de9a750bf7c31b3b53865c71b9cdb8db57d401 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.cc
+++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.cc
@@ -25,28 +25,29 @@ namespace AMDiS { namespace base_problems {
 
 CahnHilliardNavierStokes_TwoPhase::CahnHilliardNavierStokes_TwoPhase(const std::string &name_) :
   super(name_),
+  reinit(NULL),
+  velocity(NULL),
   useMobility(false),
   useConservationForm(false),
-  doubleWell(1),
-  gamma(1.0),
-  eps(0.1),
-  minusEps(-0.1),
-  epsInv(10.0),
-  minusEpsInv(-10.0),
-  epsSqr(0.01),
-  minusEpsSqr(-0.01),
+  dim(2),
   laplaceType(0),
   nonLinTerm(3),
+  doubleWell(1),
   oldTimestep(0.0),
   viscosity1(1.0),
   viscosity2(1.0),
   density1(1.0),
   density2(1.0),
   c(0.0),
+  gamma(1.0),
+  eps(0.1),
+  minusEps(-0.1),
+  epsInv(10.0),
+  minusEpsInv(-10.0),
+  epsSqr(0.01),
+  minusEpsSqr(-0.01),
   sigma(1.0),
   surfaceTension(1.0),
-  reinit(NULL),
-  velocity(NULL),
   fileWriter(NULL)
 {
   // parameters for CH
@@ -97,8 +98,7 @@ CahnHilliardNavierStokes_TwoPhase::CahnHilliardNavierStokes_TwoPhase(const std::
 
 
 CahnHilliardNavierStokes_TwoPhase::~CahnHilliardNavierStokes_TwoPhase() 
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::~CahnHilliardNavierStokes_TwoPhase()");
-
+{
   if (reinit != NULL) {
     delete reinit;
     reinit = NULL;
@@ -113,8 +113,7 @@ CahnHilliardNavierStokes_TwoPhase::~CahnHilliardNavierStokes_TwoPhase()
 
 
 void CahnHilliardNavierStokes_TwoPhase::initData()
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::initData()");
-
+{
   dim = getMesh()->getDim();
   // create instance redistancing class
   reinit = new HL_SignedDistTraverse("reinit", dim);
@@ -128,8 +127,7 @@ void CahnHilliardNavierStokes_TwoPhase::initData()
 
 
 void CahnHilliardNavierStokes_TwoPhase::transferInitialSolution(AdaptInfo *adaptInfo)
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::transferInitialSolution()");
-
+{
   calcVelocity();
   
   for (int i = 0; i < 2+dow; i++)
@@ -394,8 +392,7 @@ void CahnHilliardNavierStokes_TwoPhase::addMaterialDerivativeNonConst(int i)
 
 
 void CahnHilliardNavierStokes_TwoPhase::addLaplaceTermConst(int i)
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::addLaplaceTerm()");
-
+{
   /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
   if (laplaceType == 1) {
     for (size_t j = 0; j < dow; ++j) {
@@ -413,8 +410,7 @@ void CahnHilliardNavierStokes_TwoPhase::addLaplaceTermConst(int i)
 
 
 void CahnHilliardNavierStokes_TwoPhase::addLaplaceTermNonConst(int i)
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::addLaplaceTerm()");
-
+{
   /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
   if (laplaceType == 1) {
     for (size_t j = 0; j < dow; ++j) {
@@ -437,8 +433,7 @@ void CahnHilliardNavierStokes_TwoPhase::addLaplaceTermNonConst(int i)
 
 
 void CahnHilliardNavierStokes_TwoPhase::closeTimestep(AdaptInfo *adaptInfo)
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::closeTimestep()");
-
+{
   calcVelocity();
   fileWriter->writeFiles(adaptInfo, false);
   writeFiles(adaptInfo, false);
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.h b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.h
index a211c43e95d5a2cfcba72a0f176d23783f86b1ec..d317a546c7b88bb1b579c7a9f2044bd8a6cee639 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.h
+++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.h
@@ -37,10 +37,10 @@ public: // public methods
   CahnHilliardNavierStokes_TwoPhase(const std::string &name_);
   ~CahnHilliardNavierStokes_TwoPhase();
 
-  virtual void initData();
+  virtual void initData() override;
 
-  virtual void transferInitialSolution(AdaptInfo *adaptInfo);  
-  virtual void closeTimestep(AdaptInfo *adaptInfo);
+  virtual void transferInitialSolution(AdaptInfo *adaptInfo) override;  
+  virtual void closeTimestep(AdaptInfo *adaptInfo) override;
   
   double getEpsilon() { return eps; }
   int getDoubleWellType() { return doubleWell; }
@@ -50,20 +50,15 @@ public: // public methods
     return velocity;
   }
   
-  virtual void fillOperators();
-  virtual void fillBoundaryConditions() {}
+  virtual void fillOperators() override;
+  virtual void fillBoundaryConditions() override {}
   
   virtual void addMaterialDerivativeConst(int i);
   virtual void addMaterialDerivativeNonConst(int i);
   virtual void addLaplaceTermConst(int i);
   virtual void addLaplaceTermNonConst(int i);
 
-protected: // protected variables
-
-  HL_SignedDistTraverse *reinit;
-
-  DOFVector<WorldVector<double> >* velocity;
-
+protected: // methods
   void calcVelocity()
   {
     if (dow == 1)
@@ -73,6 +68,12 @@ protected: // protected variables
     else if (dow == 3)
       transformDOF(prob->getSolution()->getDOFVector(2), prob->getSolution()->getDOFVector(3), prob->getSolution()->getDOFVector(4), velocity, new AMDiS::Vec3WorldVec<double>);
   }
+
+protected: // protected variables
+
+  HL_SignedDistTraverse *reinit;
+
+  DOFVector<WorldVector<double> >* velocity;
   
   bool useMobility;
   bool useConservationForm;
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
index 65250f50f441e009e15c013c3f46598993c725e7..67cb978a019d9e921679104c941eec9da999b38b 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
+++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
@@ -100,8 +100,7 @@ CahnHilliardNavierStokes_TwoPhase_RB::CahnHilliardNavierStokes_TwoPhase_RB(const
 
 
 CahnHilliardNavierStokes_TwoPhase_RB::~CahnHilliardNavierStokes_TwoPhase_RB() 
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::~CahnHilliardNavierStokes_TwoPhase_RB()");
-
+{
   if (velocity != NULL) {
     delete velocity;
     velocity = NULL;
@@ -112,8 +111,7 @@ CahnHilliardNavierStokes_TwoPhase_RB::~CahnHilliardNavierStokes_TwoPhase_RB()
 
 
 void CahnHilliardNavierStokes_TwoPhase_RB::initData()
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::initData()");
-
+{
   dim = getMesh()->getDim();
   
   if (velocity == NULL)
@@ -125,8 +123,7 @@ void CahnHilliardNavierStokes_TwoPhase_RB::initData()
 
 
 void CahnHilliardNavierStokes_TwoPhase_RB::transferInitialSolution(AdaptInfo *adaptInfo)
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::transferInitialSolution()");
-
+{
   calcVelocity();
   
   for (int i = 0; i < 2+dow; i++)
@@ -446,8 +443,7 @@ void CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()
 
 
 void CahnHilliardNavierStokes_TwoPhase_RB::addStressTerm(int i)
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::addStressTerm()");
-
+{
     /// < nu*[grad(u)+grad(u)^t] , grad(psi) >
     if (laplaceType == 1) {
       for (unsigned j = 0; j < dow; ++j) {	
@@ -509,8 +505,7 @@ void CahnHilliardNavierStokes_TwoPhase_RB::addStressTerm(int i)
 }
 
 void CahnHilliardNavierStokes_TwoPhase_RB::closeTimestep(AdaptInfo *adaptInfo)
-{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::closeTimestep()");
-
+{
   calcVelocity();
   fileWriter->writeFiles(adaptInfo, false);
   writeFiles(adaptInfo, false);
diff --git a/extensions/base_problems/CouplingBaseProblem.h b/extensions/base_problems/CouplingBaseProblem.h
index d12de90a546124eee84ccf6a55f7c89cd115b525..303d67bfde5020e264f4387ae28cdae19acd37fc 100644
--- a/extensions/base_problems/CouplingBaseProblem.h
+++ b/extensions/base_problems/CouplingBaseProblem.h
@@ -45,9 +45,9 @@ public:
 
   ~CouplingBaseProblem() { }  
   
-  void initialize(Flag initFlag,
-		  ProblemStatSeq *adoptProblem = NULL,
-		  Flag adoptFlag = INIT_NOTHING) override
+  virtual void initialize(Flag initFlag,
+			  ProblemStatSeq *adoptProblem = NULL,
+			  Flag adoptFlag = INIT_NOTHING) override
   {  
     for (size_t i = 0; i < baseProblems.size(); i++) {
       for (size_t j = 0; j < baseProblems[i]->getNumProblems(); j++)
diff --git a/extensions/base_problems/CouplingBaseProblem2.h b/extensions/base_problems/CouplingBaseProblem2.h
index 35b802f5ee665f8e601ae39619c8508a3e9c194c..df2c140ac29dbf101efc843735f9b284db4dddba 100644
--- a/extensions/base_problems/CouplingBaseProblem2.h
+++ b/extensions/base_problems/CouplingBaseProblem2.h
@@ -138,9 +138,9 @@ public:
   
   ~CouplingBaseProblem() { }  
   
-  void initialize(Flag initFlag,
-		  ProblemStatSeq *adoptProblem = NULL,
-		  Flag adoptFlag = INIT_NOTHING) override
+  virtual void initialize(Flag initFlag,
+			  ProblemStatSeq *adoptProblem = NULL,
+			  Flag adoptFlag = INIT_NOTHING) override
   {  
     using namespace AMDiS::extensions;
     tools::FOR_EACH< detail::AddProblem >::loop2(*this, baseProblems);
diff --git a/extensions/base_problems/CouplingBaseProblem2_cxx11.h b/extensions/base_problems/CouplingBaseProblem2_cxx11.h
index ab3c9f839273bc6235c57ee711af079ab25a11dc..748fce8b572f6513aae2f4ddb75fae6b5f6385ba 100644
--- a/extensions/base_problems/CouplingBaseProblem2_cxx11.h
+++ b/extensions/base_problems/CouplingBaseProblem2_cxx11.h
@@ -130,9 +130,9 @@ public:
    * iterationInterface in this method. This order can be changed manually in the oneIteration
    * method.
    **/
-  void initialize(Flag initFlag,
-		  ProblemStatSeq *adoptProblem = NULL,
-		  Flag adoptFlag = INIT_NOTHING) override
+  virtual void initialize(Flag initFlag,
+			  ProblemStatSeq *adoptProblem = NULL,
+			  Flag adoptFlag = INIT_NOTHING) override
   {
     tools::FOR_EACH< detail::AddProblem >::loop2(*this, baseProblems);
     tools::FOR_EACH< detail::AddIterationInterface >::loop2(*this, baseProblems);
diff --git a/extensions/base_problems/DiffuseDomainFsi.cc b/extensions/base_problems/DiffuseDomainFsi.cc
index 873af228edf0ec9585589a55f7bbf9c73b4bf7a2..51ccfb8fed4b9a8fc02c0d020cea1fe9d7a4533f 100644
--- a/extensions/base_problems/DiffuseDomainFsi.cc
+++ b/extensions/base_problems/DiffuseDomainFsi.cc
@@ -23,22 +23,14 @@ using namespace std;
 
 DiffuseDomainFsi::DiffuseDomainFsi(const std::string &name_) :
   super(name_),
-  phase(NULL),
-  bcPhase(NULL),
-  phaseInv(NULL),
-  velocity(NULL),
-  displacement(NULL),
-  beta(1.0),
-  epsilon(1.e-1),
-  alpha(3.0),
-  mu(1.0),
-  viscosity(1.0),
-  density(1.0),
-  lambda(1.0),
-  rho(1.0),
-  nonLinTerm(3),
-  fileWriterPhase(NULL),
-  fileWriterVelocity(NULL),
+  dbc_factor(1.0), pressure_factor(1.0),
+  beta(1.0), epsilon(1.e-1), alpha(3.0),
+  nonLinTerm(3), oldTimestep(1.0),
+  viscosity(1.0), density(1.0),
+  mu(1.0), lambda(1.0), rho(1.0),  
+  phase(NULL), phaseInv(NULL), phaseInvDisturbed(NULL),
+  bcPhase(NULL), velocity(NULL), displacement(NULL),
+  fileWriterPhase(NULL), fileWriterVelocity(NULL),
   fileWriterDisplacement(NULL)
 {
   
@@ -61,7 +53,7 @@ DiffuseDomainFsi::DiffuseDomainFsi(const std::string &name_) :
   //							2... u^old*grad(u_i)
   //							3... u*grad(u_i^old) + u^old*grad(u_i) - u^old*grad(u_i^old)
   Parameters::get(name + "->non-linear term", nonLinTerm); 
-};
+}
 
 DiffuseDomainFsi::~DiffuseDomainFsi()
 {
@@ -79,7 +71,7 @@ DiffuseDomainFsi::~DiffuseDomainFsi()
 }
 
 void DiffuseDomainFsi::initData()
-{
+{ FUNCNAME("DiffuseDomainFsi::initData()");
   super::initData();
 
   phaseInv = new DOFVector<double>(getFeSpace(0), "phaseInv");
@@ -88,15 +80,14 @@ void DiffuseDomainFsi::initData()
   displacement = new DOFVector<WorldVector<double> >(getFeSpace(dow+1), "displacement");
   
   dbc_factor = beta/std::pow(epsilon, alpha);
-  MSG_DBG("dbc_factor = %f\n", dbc_factor);
+  MSG("dbc_factor = %f\n", dbc_factor);
   pressure_factor = -(lambda+(2.0/3.0)*mu);
   MSG("pressure_factor = %f\n", pressure_factor);
 }
 
 
 void DiffuseDomainFsi::transferInitialSolution(AdaptInfo *adaptInfo)
-{ FUNCNAME("DiffuseDomainFsi::transferInitialSolution()");
-
+{ 
   fileWriterPhase = new FileWriter(name + "->phase->output", getFeSpace()->getMesh(), phase);
   fileWriterVelocity = new FileVectorWriter(name + "->velocity->output", getFeSpace()->getMesh(), velocity);
   fileWriterDisplacement = new FileVectorWriter(name + "->displacement->output", getFeSpace()->getMesh(), displacement);
@@ -257,8 +248,7 @@ void DiffuseDomainFsi::fillOperators()
 
 
 void DiffuseDomainFsi::addLaplaceTerm(int i)
-{ FUNCNAME("DiffuseDomainFsi::addLaplaceTerm()");
-    
+{     
   for (size_t j = 0; j < dow; ++j) {
 	/// < nu*[grad(u)+grad(u)^t] , grad(psi) >
 	Operator *opLaplaceUi1 = new Operator(prob->getFeSpace(i), prob->getFeSpace(j));
diff --git a/extensions/base_problems/DiffuseDomainFsi.h b/extensions/base_problems/DiffuseDomainFsi.h
index c794ebcf72c97f2d0d70dad7877b17cc9ef65416..f49ca6aa3c661c968960e9e033ac6a48146f9280 100644
--- a/extensions/base_problems/DiffuseDomainFsi.h
+++ b/extensions/base_problems/DiffuseDomainFsi.h
@@ -43,12 +43,12 @@ public: // typedefs
 public: // methods
 
   DiffuseDomainFsi(const std::string &name_);
-  ~DiffuseDomainFsi();
+  virtual ~DiffuseDomainFsi();
   
-  virtual void initData();
-  virtual void transferInitialSolution(AdaptInfo *adaptInfo);  
-  virtual void initTimestep(AdaptInfo *adaptInfo);
-  virtual void closeTimestep(AdaptInfo *adaptInfo);
+  virtual void initData() override;
+  virtual void transferInitialSolution(AdaptInfo *adaptInfo) override;  
+  virtual void initTimestep(AdaptInfo *adaptInfo) override;
+  virtual void closeTimestep(AdaptInfo *adaptInfo) override;
   
   // === getting/setting methods ===
   
diff --git a/extensions/base_problems/NavierStokesPhase_TaylorHood.h b/extensions/base_problems/NavierStokesPhase_TaylorHood.h
index 377b81dc58f39e13cba6d1da8c965fea4bf32eae..39efc40b53958f2815f083a33a5eae85741f39b3 100644
--- a/extensions/base_problems/NavierStokesPhase_TaylorHood.h
+++ b/extensions/base_problems/NavierStokesPhase_TaylorHood.h
@@ -94,8 +94,8 @@ public: // methods
     epsilon = epsilon_;
   }
 
-  virtual void fillOperators();
-  virtual void fillBoundaryConditions();
+  virtual void fillOperators() override;
+  virtual void fillBoundaryConditions() override;
   virtual void addLaplaceTerm(int i);
 
 protected: // variables
diff --git a/extensions/base_problems/NavierStokes_TH_MultiPhase.h b/extensions/base_problems/NavierStokes_TH_MultiPhase.h
index ae267b49195bd439f7e074c6232b2faa623fb90b..8769c4cae4a48117fa42231aa749fc1b1ce94b1b 100644
--- a/extensions/base_problems/NavierStokes_TH_MultiPhase.h
+++ b/extensions/base_problems/NavierStokes_TH_MultiPhase.h
@@ -50,14 +50,14 @@ namespace detail {
     * viscosity/density DOFVectors,
     * initTimeInterface of super
     **/
-    void initData() override;
-    void transferInitialSolution(AdaptInfo *adaptInfo) override;
+    virtual void initData() override;
+    virtual void transferInitialSolution(AdaptInfo *adaptInfo) override;
 
     /** initTimestep of super and
     * initialization of densityPhase and viscosityPhase
     **/
-    void initTimestep(AdaptInfo *adaptInfo) override;
-    void closeTimestep(AdaptInfo *adaptInfo) override;
+    virtual void initTimestep(AdaptInfo *adaptInfo) override;
+    virtual void closeTimestep(AdaptInfo *adaptInfo) override;
 
     /** pointer to the multiPhase DOFVector, that
     * indicates the different phases. Will be initialized
@@ -65,9 +65,9 @@ namespace detail {
     **/
     void setMultiPhase(DOFVector<double>* p) { multiPhase=p; }
 
-    void fillOperators() override;
+    virtual void fillOperators() override;
 
-    void addLaplaceTerm(int i) override;
+    virtual void addLaplaceTerm(int i) override;
     
     // get-methods
     
diff --git a/extensions/base_problems/NavierStokes_TH_MultiPhase.hh b/extensions/base_problems/NavierStokes_TH_MultiPhase.hh
index 9beca9369970e4790ed30b784e5a60fae163b934..62c0dec97c61225fde202e91990a948fe627ce5d 100644
--- a/extensions/base_problems/NavierStokes_TH_MultiPhase.hh
+++ b/extensions/base_problems/NavierStokes_TH_MultiPhase.hh
@@ -32,8 +32,8 @@ NavierStokes_TH_MultiPhase<P>::NavierStokes_TH_MultiPhase(const std::string &nam
     densityPhase(NULL),
     viscosityPhase(NULL),
     multiPhase(NULL),
-    stepSolution(NULL),
-    oldSolution(NULL)
+    oldSolution(NULL),
+    stepSolution(NULL)
 { FUNCNAME("NavierStokes_TH_MultiPhase::_constructor()");
 
   Initfile::get(self::name + "->viscosity1", viscosity1); // viscosity of fluid 1
diff --git a/extensions/base_problems/NavierStokes_TaylorHood.h b/extensions/base_problems/NavierStokes_TaylorHood.h
index 0e0ea152b012b78f3aa5c6aef42f40da9fe467c9..740cc68b8d1155da3053e5cc57837ae215fc5a14 100644
--- a/extensions/base_problems/NavierStokes_TaylorHood.h
+++ b/extensions/base_problems/NavierStokes_TaylorHood.h
@@ -45,15 +45,15 @@ namespace detail
 
     NavierStokes_TaylorHood(const std::string &name_, bool createProblem = true);
 
-    ~NavierStokes_TaylorHood();
+    virtual ~NavierStokes_TaylorHood();
 
-    void initData() override;
+    virtual void initData() override;
 
-    void transferInitialSolution(AdaptInfo *adaptInfo) override;  
+    virtual void transferInitialSolution(AdaptInfo *adaptInfo) override;  
 
-    void closeTimestep(AdaptInfo *adaptInfo) override;
+    virtual void closeTimestep(AdaptInfo *adaptInfo) override;
     
-    void writeFiles(AdaptInfo *adaptInfo, bool force = false) override;
+    virtual void writeFiles(AdaptInfo *adaptInfo, bool force = false) override;
 
     // === getting/setting methods ===
 
@@ -74,7 +74,7 @@ namespace detail
     double getDensity() { return density; }
     double getViscosity() { return viscosity; }
 
-    void fillOperators() override;
+    virtual void fillOperators() override;
     
     virtual void addLaplaceTerm(int i);
 
diff --git a/extensions/base_problems/PhaseFieldCrystal.h b/extensions/base_problems/PhaseFieldCrystal.h
index e68333957d07fcec72d26bab7ddd487cec393b94..e72ecb5f596863ea83abec1e50408c55c52b9edc 100644
--- a/extensions/base_problems/PhaseFieldCrystal.h
+++ b/extensions/base_problems/PhaseFieldCrystal.h
@@ -44,13 +44,13 @@ namespace detail {
     DOFVector<double>* getDensity() { return self::prob->getSolution(1); }
     double calcEnergy();
 
-    void initData() override ;
-    void transferInitialSolution(AdaptInfo *adaptInfo) override;
+    virtual void initData() override ;
+    virtual void transferInitialSolution(AdaptInfo *adaptInfo) override;
     
-    void closeTimestep(AdaptInfo *adaptInfo) override;
+    virtual void closeTimestep(AdaptInfo *adaptInfo) override;
     
-    void fillOperators() override ;
-    void fillBoundaryConditions() override {};
+    virtual void fillOperators() override ;
+    virtual void fillBoundaryConditions() override {};
     
     void rejectTimestep()
     {
diff --git a/extensions/base_problems/PhaseFieldCrystal_Phase.h b/extensions/base_problems/PhaseFieldCrystal_Phase.h
index fffb75accb9043249b0ddc8f2e07c5f6399eea6b..bf17b241f5cfa9bb3d0fbeaa6591cfc577bf4f48 100644
--- a/extensions/base_problems/PhaseFieldCrystal_Phase.h
+++ b/extensions/base_problems/PhaseFieldCrystal_Phase.h
@@ -58,7 +58,7 @@ namespace AMDiS { namespace base_problems {
 	*phaseDisturbed << max(valueOf(phase), 1.e-6);
       };
 
-      void fillOperators() override;
+      virtual void fillOperators() override;
       
       double calcEnergy();
 
diff --git a/extensions/base_problems/PolarizationField.h b/extensions/base_problems/PolarizationField.h
index 7ca6159ab9444dd028bd144668eedaad7ed43d34..8661510a4591b30713a3acf46dbd7ae99950b1a3 100644
--- a/extensions/base_problems/PolarizationField.h
+++ b/extensions/base_problems/PolarizationField.h
@@ -62,16 +62,16 @@ namespace detail
     ~PolarizationField();
 
     /// initialize the vectorField and corresponding fileWriter
-    void initData() override;
+    virtual void initData() override;
 
     /// calls \ref calcVectorField and \ref super::transferInitialSolution
-    void transferInitialSolution(AdaptInfo *adaptInfo) override;  
+    virtual void transferInitialSolution(AdaptInfo *adaptInfo) override;  
 
     /// calls \ref calcVectorField and \ref super::closeTimestep
-    void closeTimestep(AdaptInfo *adaptInfo) override;
+    virtual void closeTimestep(AdaptInfo *adaptInfo) override;
     
     /// write the solution and the vectorField
-    void writeFiles(AdaptInfo *adaptInfo, bool force = false) override;
+    virtual void writeFiles(AdaptInfo *adaptInfo, bool force = false) override;
 
     // === getting/setting methods ===
 
@@ -92,7 +92,7 @@ namespace detail
     }
 
     /// implementation of BaseProblem::fillOperators
-    void fillOperators() override;
+    virtual void fillOperators() override;
     
     /// used in fillOperators to add the term (grad(P_i), grad(psi))
     virtual void fillLaplacian();
diff --git a/extensions/base_problems/QuasiCrystal.h b/extensions/base_problems/QuasiCrystal.h
index cf80ab139efba5bc4945611229e9a8d29a3ab0c7..8f02bf804fdddb364316ba8f4f51be47d47e2aeb 100644
--- a/extensions/base_problems/QuasiCrystal.h
+++ b/extensions/base_problems/QuasiCrystal.h
@@ -26,10 +26,10 @@ namespace detail {
     QuasiCrystal(const std::string &name_);
     ~QuasiCrystal() {};
 
-    void fillOperators() override;
-    void fillBoundaryConditions() override {};
+    virtual void fillOperators() override;
+    virtual void fillBoundaryConditions() override {};
     
-    void finalizeData() override;
+    virtual void finalizeData() override;
 
   protected:
 
diff --git a/extensions/base_problems/QuasiCrystal_RB.cc b/extensions/base_problems/QuasiCrystal_RB.cc
index 7967a31e5a2ff84446e7d139d8a307a517c2dc07..f95df93ff121b7b60c0f31905f6ce2c98b35c968 100644
--- a/extensions/base_problems/QuasiCrystal_RB.cc
+++ b/extensions/base_problems/QuasiCrystal_RB.cc
@@ -23,12 +23,11 @@ QuasiCrystal_RB::QuasiCrystal_RB(const std::string &name_) :
   Parameters::get(name + "->density", density);
   Parameters::get(name + "->alpha", alpha);
   minusC = -c;
-};
+}
 
 
 void QuasiCrystal_RB::fillOperators()
-{ FUNCNAME("QuasiCrystal_RB::fillOperators()");
-  
+{ 
   // mass matrices
   Operator *opM = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
   opM->addTerm(new Simple_ZOT);
@@ -150,11 +149,10 @@ void QuasiCrystal_RB::fillOperators()
   
   prob->addVectorOperator(opM_4, 4);
   prob->addMatrixOperator(opM, 4, 4, &minus1, &minus1); // theta
-};
+}
 
 
 void QuasiCrystal_RB::fillBoundaryConditions()
-{ FUNCNAME("QuasiCrystal_RB::fillBoundaryConditions()");
-};
+{ }
 
-} }
\ No newline at end of file
+} }
diff --git a/extensions/demo/pfc/src/pfc.cc b/extensions/demo/pfc/src/pfc.cc
index d57cd94910da7ce71a55bc3ee62b2db6cac56eeb..46a8e863e06a62abcf6aa5d4a05637a6998aabd2 100644
--- a/extensions/demo/pfc/src/pfc.cc
+++ b/extensions/demo/pfc/src/pfc.cc
@@ -12,7 +12,7 @@
 #include "preconditioner/MTLPreconPfc.h"
 #endif
 
-#include "GenericOperatorTerm.h"
+#include "Expressions.h"
 
 using namespace AMDiS;
 
@@ -25,39 +25,39 @@ public:
   PfcPC(std::string name) : PhaseFieldCrystal(name) { }
   
   /// initialize the preconditioners, i.e. set parameters
-  void initData() 
+  void finalizeData() override
   {
-    super::initData();
+    super::finalizeData();
 
     // sequential PFC preconditioner
 #if (defined HAVE_SEQ_PETSC) || (defined HAVE_PETSC)
-    PetscPreconPfc* runner = dynamic_cast<PetscPreconPfc*>(prob->getSolver()->getRunner());
-    if (runner) {
-      dynamic_cast<PetscSolver<PetscPreconPfc>*>(prob->getSolver())->setNested(true);
-      runner->setData(getTau());	
+    PetscPreconPfc* precon = dynamic_cast<PetscPreconPfc*>(prob->getSolver()->getRightPrecon());
+    if (precon) {
+      precon->setData(getTau(), M0);	
     }
-#endif
-
+//     PetscPreconPfcDiag* precon2 = dynamic_cast<PetscPreconPfcDiag*>(prob->getSolver()->getRightPrecon());
+//     if (precon2) {
+//       precon2->setData(getTau(), M0);	
+    }
+#elif !defined(HAVE_PARALLEL_DOMAIN_AMDIS)
+     MTLPreconPfc* precon = dynamic_cast<MTLPreconPfc*>(prob->getSolver()->getRightPrecon());
+     if (precon)
+	precon->setData(getTau(), M0);	
+#else
     // parallel PFC preconditioner
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
     Parallel::PetscSolverPfc* solver = dynamic_cast<Parallel::PetscSolverPfc*>(prob->getSolver());
     if (solver)
-      solver->setData(getTau());
-#endif
-      
-    // sequential PFC preconditioner using MTL    
-#if (!defined HAVE_SEQ_PETSC) && (!defined HAVE_PETSC) && (!defined HAVE_PARALLEL_DOMAIN_AMDIS)
-    using AMDiS::extensions::MTLPreconPfc;
-    MTLPreconPfc* precon = dynamic_cast<MTLPreconPfc*>(prob->getSolver()->getRightPrecon());
-    if (precon)
-      precon->setData(getTau());	
+      solver->setData(getTau(), M0);
+    
+//     Parallel::PetscSolverPfcDiag* solver2 = dynamic_cast<Parallel::PetscSolverPfcDiag*>(prob->getSolver());
+//     if (solver2)
+//       solver2->setData(getTau(), M0);
 #endif
   }
   
   /// generate initial solution for evolution equation
-  void solveInitialProblem(AdaptInfo *adaptInfo)
-  { FUNCNAME("PFC_Demo::solveInitialProblem()");
-
+  void solveInitialProblem(AdaptInfo *adaptInfo) override
+  { 
     Flag initFlag = initDataFromFile(adaptInfo);
     if (initFlag.isSet(DATA_ADOPTED))
       return;
@@ -79,16 +79,14 @@ int main(int argc, char** argv)
   
   // add preconditioner / solver to the parameter list. Must be added before problem is initialized.
 #if (defined HAVE_SEQ_PETSC) || (defined HAVE_PETSC)
-  CreatorMap<LinearSolver>::addCreator("petsc_pfc", new PetscSolver<PetscPreconPfc>::Creator);
-#endif
-  
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-  CreatorMap<LinearSolver>::addCreator("p_petsc_pfc", new Parallel::PetscSolverPfc::Creator);
-#endif
-  
-#if (!defined HAVE_SEQ_PETSC) && (!defined HAVE_PETSC) && (!defined HAVE_PARALLEL_DOMAIN_AMDIS)
-  using AMDiS::extensions::MTLPreconPfc;
-  CreatorMap<BasePreconditioner>::addCreator("pfc", new MTLPreconPfc::Creator);
+  CreatorMap<PetscPreconditionerNested>::addCreator("pfc", new PetscPreconPfc::Creator);
+//   CreatorMap<PetscPreconditionerNested>::addCreator("pfc_diag", new PetscPreconPfcDiag::Creator);
+#elif !defined(HAVE_PARALLEL_DOMAIN_AMDIS)
+  CreatorMap<typename MTLPreconPfc::precon_base>::addCreator("pfc", new MTLPreconPfc::Creator);
+//   CreatorMap<typename MTLPreconPfc_Diag::base_precon>::addCreator("pfc", new MTLPreconPfc::Creator);
+#else
+  CreatorMap<LinearSolverInterface>::addCreator("p_petsc_pfc", new Parallel::PetscSolverPfc::Creator);
+//   CreatorMap<LinearSolverInterface>::addCreator("p_petsc_pfc_diag", new Parallel::PetscSolverPfcDiag::Creator);
 #endif
   
   // create and initialize the PFC BaseProblem
diff --git a/extensions/preconditioner/MTLPreconPfc.h b/extensions/preconditioner/MTLPreconPfc.h
index 21c16a66dd7b42207405b9f86299d78a8ba07d77..508876b8183a8778b7aae82182bc4ebc6e49ac4d 100644
--- a/extensions/preconditioner/MTLPreconPfc.h
+++ b/extensions/preconditioner/MTLPreconPfc.h
@@ -30,9 +30,9 @@ using namespace AMDiS::MTLTypes;
 template<typename MatrixType>
 struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
 {
-  typedef BlockPreconditioner<MatrixType>                           super;
-  typedef ITL_BasePreconditioner<MatrixType, MTLTypes::MTLVector>   precon_base;
-  typedef MTLPreconPfcBase<MatrixType>                              self;
+  typedef BlockPreconditioner<MatrixType>   super;
+  typedef typename super::precon_base       precon_base;
+  typedef MTLPreconPfcBase<MatrixType>      self;
  
   class Creator : public CreatorInterfaceName<precon_base>
   {
@@ -48,6 +48,8 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
     : super(),
       name(name), 
       tau(NULL),
+      M0(1.0),
+      matM(NULL), matL(NULL), matL0(NULL),
       PId(NULL), PDiagM(NULL), PDiagMpL(NULL), PDiagMpL2(NULL),
       solverM(NULL), solverMpL(NULL), solverMpL2(NULL)
   { }
@@ -61,13 +63,22 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
   void setData(Data* data)
   {
     tau = data->tau;
+    M0 = data->M0;
     matM = &data->getM();
     matL = &data->getL();
+    matL0 = &data->getL0();
   }
+    
+  void setData(double* tau_, double M0_ = 1.0)
+  {
+    tau = tau_;
+    M0 = M0_;
+  }
+
   
-  virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix);
+  virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override;
   
-  virtual void exit()
+  virtual void exit() override
   {    
     if (PId) { delete PId; PId = NULL; }
     if (PDiagM) { delete PDiagM; PDiagM = NULL; }
@@ -78,7 +89,7 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
     if (solverMpL2) { delete solverMpL2; solverMpL2 = NULL; }
   }
   
-  virtual void solve(const MTLVector& b, MTLVector& x) const;  
+  virtual void solve(const MTLVector& b, MTLVector& x) const override;  
     
   template <typename VectorX, typename VectorB>
   void solveM(VectorX& x, const VectorB& b, int nIter = -1) const
@@ -137,14 +148,21 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
     }
   }
   
+  const MTLMatrix& getM() const { return matM ? *matM : super::A->getSubMatrix(2,2); }
+  const MTLMatrix& getL0() const { return matL0 ? *matL0 : super::A->getSubMatrix(1,0); }
+  const MTLMatrix& getL() const { return matL ? *matL : super::A->getSubMatrix(2,1); }
   
-  const MTLMatrix& getM() const { return *matM; }
-  const MTLMatrix& getL() const { return *matL; }
   double getTau() const { return *tau; }
   
 protected: 
+  std::string name;
+  
+  double* tau;
+  double M0;
+  
   MTLMatrix* matM;
   MTLMatrix* matL;
+  MTLMatrix* matL0;
   
   MTLMatrix MpL;
   MTLMatrix MpL2;
@@ -161,10 +179,6 @@ protected:
   mtl::matrix::umfpack::solver<MTLMatrix> *solverM;
   mtl::matrix::umfpack::solver<MTLMatrix> *solverMpL;
   mtl::matrix::umfpack::solver<MTLMatrix> *solverMpL2;
-
-  std::string name;
-  
-  double* tau;
 };
 
 typedef MTLPreconPfcBase<MTLTypes::MTLMatrix> MTLPreconPfc;
diff --git a/extensions/preconditioner/MTLPreconPfc.hh b/extensions/preconditioner/MTLPreconPfc.hh
index d10df58f22208f586f24d0398fa0ae7085876148..396f5ea97c43947ee408a797c54d557595d96763 100644
--- a/extensions/preconditioner/MTLPreconPfc.hh
+++ b/extensions/preconditioner/MTLPreconPfc.hh
@@ -27,11 +27,11 @@ void MTLPreconPfcBase<MatrixType>::init(const SolverMatrix<Matrix<DOFMatrix*> >&
   typedef typename mtl::Collection<MTLMatrix>::value_type value_type;
   size_type n = num_rows(getM());
   
-  double delta = std::sqrt(getTau());
+  double delta = std::sqrt(M0 * getTau());
   
   // helper-matrix MpL = M + delta*L
   MpL.change_dim(n, n);
-  MpL = getM() + delta*getL();
+  MpL = getM() + (1.0/delta)*getL0();
 
   // helper-matrix MpL = M + sqrt(delta)*L
   MpL2.change_dim(n, n);
@@ -81,11 +81,11 @@ void MTLPreconPfcBase<MatrixType>::solve(const MTLVector& b, MTLVector& x) const
   const MTLVector b1(b[super::getRowRange(1)]);
   const MTLVector b2(b[super::getRowRange(2)]); 
   
-  double delta = std::sqrt(getTau());
+  double delta = std::sqrt(M0 * getTau());
   
   solveM(y0, b0);		// M*y0 = b0
-  y1 = getL() * y0;		// y1 = K*y0
-  tmp = b1 - getTau()*y1;	// tmp := b1 - tau*y1
+  y1 = getL0() * y0;		// y1 = K*y0
+  tmp = b1 - y1;		// tmp := b1 - tau*y1
 
   solveMpL(y1, tmp);		// (M + delta*K) * y1 = tmp
   x0 = y0 + (1.0/delta)*y1;	// x0 = y0 + (1/delta)*y1
diff --git a/extensions/preconditioner/PetscPreconPfc.cc b/extensions/preconditioner/PetscPreconPfc.cc
index 22b06545a28dd439c6913f46ba26aca621133423..bb4fc21fe384a7352e0ce61ec407418292313b71 100644
--- a/extensions/preconditioner/PetscPreconPfc.cc
+++ b/extensions/preconditioner/PetscPreconPfc.cc
@@ -47,7 +47,7 @@ namespace AMDiS {
 
     KSPSolve(data->kspM, b1, y1); 			// M*y1 = b1    
     MatMult(data->matK, y1, tmp); 			// tmp := K*y1
-    VecAYPX(tmp, -(data->tau), b2);			// tmp := b2 - tau*tmp
+    VecAYPX(tmp, -sqr(data->delta), b2);			// tmp := b2 - tau*tmp
     
     KSPSolve(data->kspMpK, tmp, y2); 			// (M+sqrt(tau)K)*y2 = tmp
     
@@ -83,7 +83,7 @@ namespace AMDiS {
     PCShellSetContext(pc, &data);
     
   
-    double delta = sqrt(*tau);
+    double delta = std::sqrt((*tau) * M0);
     
     MatNestGetSubMat(A.matrix, 2, 2, &data.matM);
     MatNestGetSubMat(A.matrix, 2, 1, &data.matK);
@@ -92,7 +92,7 @@ namespace AMDiS {
     MatAXPY(MpK, delta, data.matK, SAME_NONZERO_PATTERN);
     
     MatDuplicate(data.matM, MAT_COPY_VALUES, &MpK2);
-    MatAXPY(MpK2, sqrt(delta), data.matK, SAME_NONZERO_PATTERN);
+    MatAXPY(MpK2, std::sqrt(delta), data.matK, SAME_NONZERO_PATTERN);
     
     // init sub-solvers
     Runner::createSubSolver(data.kspM, data.matM, "mass_");
diff --git a/extensions/preconditioner/PetscPreconPfc.h b/extensions/preconditioner/PetscPreconPfc.h
index a61ba554b2011950d58ab2b0ddbd7242bea04f91..a2036addc0a57acfdd91ae7f57d8edbb60502772 100644
--- a/extensions/preconditioner/PetscPreconPfc.h
+++ b/extensions/preconditioner/PetscPreconPfc.h
@@ -49,12 +49,14 @@ namespace AMDiS {
   public:
     PetscPreconPfc()
       : super("pfc_", "pfc"),
-	tau(NULL)
+	tau(NULL),
+	M0(1.0)
     { }
 
-    void setData(double *tauPtr)
+    void setData(double *tau_, double M0_ = 1.0)
     {
-      tau = tauPtr;
+      tau = tau_;
+      M0 = M0_;
     }
 
     void init(PC pc, const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override;
@@ -66,6 +68,7 @@ namespace AMDiS {
     PfcData data;
 
     double *tau;
+    double M0;
   };
 
 }
diff --git a/extensions/preconditioner/PetscSolverPfc.cc b/extensions/preconditioner/PetscSolverPfc.cc
index 6101c175092439fb047b295f3bf985e080366717..64931cd7429b13aeec01da1b0b98205f048e8dd6 100644
--- a/extensions/preconditioner/PetscSolverPfc.cc
+++ b/extensions/preconditioner/PetscSolverPfc.cc
@@ -18,6 +18,7 @@
 #include "PetscSolverPfc.h"
 #include "parallel/PetscHelper.h"
 #include "parallel/PetscSolverGlobalMatrix.h"
+#include "Expressions.h"
 
 namespace AMDiS { namespace Parallel {
 
@@ -47,24 +48,24 @@ namespace AMDiS { namespace Parallel {
     VecDuplicate(b1, &y2);
     VecDuplicate(b1, &tmp);
 
-    KSPSolve(data->kspM, b1, y1); 			// M*y1 = b1    
-    MatMult(data->matK, y1, tmp); 			// tmp := K*y1
-    VecAYPX(tmp, -(data->tau), b2);			// tmp := b2 - tau*tmp
+    KSPSolve(data->kspM, b1, y1); 		// M*y1 = b1    
+    MatMult(data->matK, y1, tmp); 		// tmp := K*y1
+    VecAYPX(tmp, -sqr(data->delta), b2);		// tmp := b2 - M0*tau*tmp
     
-    KSPSolve(data->kspMpK, tmp, y2); 			// (M+sqrt(tau)K)*y2 = tmp
+    KSPSolve(data->kspMpK, tmp, y2); 		// (M+sqrt(tau)K)*y2 = tmp
     
-    MatMult(data->matM, y2, tmp);			// tmp := M*y2
-    KSPSolve(data->kspMpK2, tmp, x2); 			// MpK2*x2 = tmp
-    MatMult(data->matM, x2, tmp);			// tmp := M*x2
-    KSPSolve(data->kspMpK2, tmp, x2); 			// MpK2*x2 = tmp
+    MatMult(data->matM, y2, tmp);		// tmp := M*y2
+    KSPSolve(data->kspMpK2, tmp, x2); 		// MpK2*x2 = tmp
+    MatMult(data->matM, x2, tmp);		// tmp := M*x2
+    KSPSolve(data->kspMpK2, tmp, x2); 		// MpK2*x2 = tmp
 
-    VecCopy(x2, x1); 					// x1 := x2
+    VecCopy(x2, x1); 				// x1 := x2
     VecAXPBYPCZ(x1, 1.0, 1.0/(data->delta), -1.0/(data->delta), y1, y2);	// x1 = 1*y1 + factor*y2 - factor*x1
 
-    MatMult(data->matK, x2, tmp); 			// tmp := K*x2
-    VecAYPX(tmp, -1.0, b3);				// tmp := b3 - tmp
+    MatMult(data->matK, x2, tmp); 		// tmp := K*x2
+    VecAYPX(tmp, -1.0, b3);			// tmp := b3 - tmp
     
-    KSPSolve(data->kspM, tmp, x3); 			// M*x3 = tmp    
+    KSPSolve(data->kspM, tmp, x3); 		// M*x3 = tmp    
 
     VecDestroy(&y1);
     VecDestroy(&y2);
@@ -105,16 +106,15 @@ namespace AMDiS { namespace Parallel {
     PCShellSetApply(getPc(), pcPfcShell);
     PCShellSetContext(getPc(), &data);
     
-    double delta = sqrt(*tau);
+    double delta = std::sqrt((*tau) * M0);
     
     const FiniteElemSpace *feSpace = componentSpaces[0];
     
     // create mass-matrix
     DOFMatrix matrixM(feSpace, feSpace);
-    Operator massOp(feSpace, feSpace);
-    Simple_ZOT zot;
-    massOp.addTerm(&zot);
-    matrixM.assembleOperator(massOp);
+    Operator opM(feSpace, feSpace);
+    addZOT(opM, 1.0);
+    matrixM.assembleOperator(opM);
     
     solverM = createSubSolver(0, "M_");
     solverM->fillPetscMatrix(&matrixM);
@@ -123,11 +123,10 @@ namespace AMDiS { namespace Parallel {
     
     // create MpK-matrix
     DOFMatrix matrixMpK(feSpace, feSpace);
-    Operator laplaceOp2(feSpace, feSpace);
-    Simple_SOT sot2(delta);
-    laplaceOp2.addTerm(&zot);
-    laplaceOp2.addTerm(&sot2);
-    matrixMpK.assembleOperator(laplaceOp2);
+    Operator opMpK(feSpace, feSpace);
+    addZOT(opMpK, 1.0);
+    addSOT(opMpK, delta);
+    matrixMpK.assembleOperator(opMpK);
     
     solverMpK = createSubSolver(0, "MpK_");
     solverMpK->fillPetscMatrix(&matrixMpK);
@@ -137,11 +136,10 @@ namespace AMDiS { namespace Parallel {
     
     // create MpK2-matrix
     DOFMatrix matrixMpK2(feSpace, feSpace);
-    Operator laplaceOp3(feSpace, feSpace);
-    Simple_SOT sot3(sqrt(delta));
-    laplaceOp3.addTerm(&zot);
-    laplaceOp3.addTerm(&sot3);
-    matrixMpK2.assembleOperator(laplaceOp3);
+    Operator opMpK2(feSpace, feSpace);
+    addZOT(opMpK2, 1.0);
+    addSOT(opMpK2, std::sqrt(delta));
+    matrixMpK2.assembleOperator(opMpK2);
     
     solverMpK2 = createSubSolver(0, "MpK2_");
     solverMpK2->fillPetscMatrix(&matrixMpK2);
@@ -150,10 +148,9 @@ namespace AMDiS { namespace Parallel {
     
     // create laplace-matrix
     DOFMatrix matrixK(feSpace, feSpace);
-    Operator laplaceOp(feSpace, feSpace);
-    Simple_SOT sot;
-    laplaceOp.addTerm(&sot);
-    matrixK.assembleOperator(laplaceOp);
+    Operator opK(feSpace, feSpace);
+    addSOT(opK, 1.0);
+    matrixK.assembleOperator(opK);
     
     solverK = createSubSolver(0, "K_");
     solverK->fillPetscMatrix(&matrixK);
diff --git a/extensions/preconditioner/PetscSolverPfc.h b/extensions/preconditioner/PetscSolverPfc.h
index b358410f3fc98d465fbdfea4ae874c300853ccf9..5eeb98324bcfde97c1d1e305e5ee61e8892127e3 100644
--- a/extensions/preconditioner/PetscSolverPfc.h
+++ b/extensions/preconditioner/PetscSolverPfc.h
@@ -48,14 +48,16 @@ namespace AMDiS { namespace Parallel {
 	solverMpK(NULL),
 	solverMpK2(NULL),
 	useOldInitialGuess(false),
-	tau(NULL)
+	tau(NULL),
+	M0(1.0)
     { 	
       Parameters::get(initFileStr + "->use old initial guess", useOldInitialGuess);
     }
 
-    void setData(double *tauPtr)
+    void setData(double *tau_, double M0_ = 1.0)
     {
-      tau = tauPtr;
+      tau = tau_;
+      M0 = M0_;
     }
 
   protected:
@@ -73,10 +75,10 @@ namespace AMDiS { namespace Parallel {
 
     PetscSolver *solverM, *solverK, *solverMpK, *solverMpK2;
 
-    bool useOldInitialGuess;    
-    SystemVector* solution;
+    bool useOldInitialGuess;
 
     double *tau;
+    double M0;
   };
 
 } }