diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000000000000000000000000000000000000..6232a0eab400a4a6d4bee60b6f22ddeb1a1b0e36
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,21 @@
+---
+# Install external dependencies
+before_script:
+  - duneci-install-module https://gitlab.dune-project.org/fufem/dune-matrix-vector.git
+  - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-solvers.git
+  - duneci-install-module https://gitlab.dune-project.org/fufem/dune-fufem.git
+
+dune:2.9 debian-11 gcc-10 C++20:
+  variables:
+    DUNECI_BRANCH: releases/2.9
+  image: registry.dune-project.org/docker/ci/dune:2.9-debian-11-gcc-10-20
+  before_script:
+  script: duneci-standard-test
+
+dune:git ubuntu-20-04 clang-10  C++20:
+  image: registry.dune-project.org/docker/ci/dune:git-ubuntu-20.04-clang-10-20
+  script: duneci-standard-test
+
+dune:git debian-11 gcc-10  C++20:
+  image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
+  script: duneci-standard-test
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 93c6973db5303856eb277d5168eb09568431884c..cdc6dd1a1ee9b791804f4bcdf8bd0a801302d365 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -12,7 +12,6 @@
 #include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number 
 #include <dune/solvers/solvers/umfpacksolver.hh> 
 
-using namespace Dune;
 using namespace MatrixOperations;
 
 using std::shared_ptr;
@@ -29,16 +28,16 @@ public:
 	
 	using GridView = typename Basis::GridView;
 	using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
-	using ScalarRT = FieldVector< double, 1>;
-	using VectorRT = FieldVector< double, dimworld>;
-	using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
+	using ScalarRT = Dune::FieldVector< double, 1>;
+	using VectorRT = Dune::FieldVector< double, dimworld>;
+	using MatrixRT = Dune::FieldMatrix< double, dimworld, dimworld>;
 	using FuncScalar = std::function< ScalarRT(const Domain&) >;
 	using FuncVector = std::function< VectorRT(const Domain&) >;
 	using Func2Tensor = std::function< MatrixRT(const Domain&) >;
   using Func2int = std::function< int(const Domain&) >;
-	using VectorCT = BlockVector<FieldVector<double,1> >;
-	using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
-	using ElementMatrixCT = Matrix<FieldMatrix<double,1,1> >;
+	using VectorCT = Dune::BlockVector<Dune::FieldVector<double,1> >;
+	using MatrixCT = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
+	using ElementMatrixCT = Dune::Matrix<Dune::FieldMatrix<double,1,1> >;
 
 	using HierarchicVectorView = Dune::Functions::HierarchicVectorWrapper< VectorCT, double>;
 	
@@ -50,7 +49,7 @@ protected:
   const Material& material_;
 
 	fstream& log_;      // Output-log
-	const ParameterTree& parameterSet_;
+	const Dune::ParameterTree& parameterSet_;
 
 	// const FuncScalar mu_; 
   // const FuncScalar lambda_; 
@@ -61,7 +60,7 @@ protected:
 
   VectorCT x_1_, x_2_, x_3_;            // (all) Corrector coefficient vectors 
   VectorCT phi_1_, phi_2_, phi_3_;      // Corrector phi_i coefficient vectors 
-  FieldVector<double,3> m_1_, m_2_, m_3_;  // Corrector m_i coefficient vectors 
+  Dune::FieldVector<double,3> m_1_, m_2_, m_3_;  // Corrector m_i coefficient vectors
 
   MatrixRT M1_, M2_, M3_;  // (assembled) corrector matrices M_i
   const std::array<MatrixRT*, 3 > mContainer = {&M1_ , &M2_, &M3_};
@@ -98,7 +97,7 @@ public:
             const Material& material,
             double gamma,
             std::fstream& log, 
-            const ParameterTree& parameterSet)
+            const Dune::ParameterTree& parameterSet)
           : basis_(basis), 
             material_(material),
             gamma_(gamma),
@@ -133,7 +132,7 @@ public:
     ///////////////////////////////
     const shared_ptr<Basis> getBasis() {return make_shared<Basis>(basis_);}
 
-    ParameterTree getParameterSet() const {return parameterSet_;}
+    Dune::ParameterTree getParameterSet() const {return parameterSet_;}
 
     fstream* getLog(){return &log_;}
 
@@ -171,7 +170,7 @@ public:
 
 
   // Get the occupation pattern of the stiffness matrix
-  void getOccupationPattern(MatrixIndexSet& nb)
+  void getOccupationPattern(Dune::MatrixIndexSet& nb)
   {
     //  OccupationPattern:
     // | phi*phi    m*phi |
@@ -215,7 +214,7 @@ public:
     // setOneBaseFunctionToZero
     //////////////////////////////////////////////////////////////////
     if(parameterSet_.get<bool>("set_oneBasisFunction_Zero ", true)){
-      FieldVector<int,3> row;
+      Dune::FieldVector<int,3> row;
       unsigned int arbitraryLeafIndex =  parameterSet_.get<unsigned int>("arbitraryLeafIndex", 0);
       unsigned int arbitraryElementNumber =  parameterSet_.get<unsigned int>("arbitraryElementNumber", 0);
 
@@ -287,7 +286,7 @@ public:
   //   int orderQR = 1;
   //   int orderQR = 2;
   //   int orderQR = 3;
-    const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+    const auto& quad = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
     
   //   double elementContribution = 0.0;
     
@@ -301,11 +300,11 @@ public:
       const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
       const auto integrationElement = geometry.integrationElement(quadPos);
 
-      std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
+      std::vector< Dune::FieldMatrix< double, 1, dim> > referenceGradients;
       localFiniteElement.localBasis().evaluateJacobian(quadPos,
                                                       referenceGradients);
       // Compute the shape function gradients on the grid element
-      std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+      std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());
 
       for (size_t i=0; i<gradients.size(); i++)
         jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
@@ -557,19 +556,19 @@ public:
   //   int orderQR = 2;
   //   int orderQR = 3;
     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);  
-    const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+    const auto& quad = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
   //   std::cout << "Quad-Rule order used: " << orderQR << std::endl;
 
     for (const auto& quadPoint : quad)
     {
-      const FieldVector<double,dim>& quadPos = quadPoint.position();
+      const Dune::FieldVector<double,dim>& quadPos = quadPoint.position();
       const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
       const double integrationElement = geometry.integrationElement(quadPos);
 
-      std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+      std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
       localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
 
-      std::vector< FieldVector< double, dim> > gradients(referenceGradients.size());  
+      std::vector< Dune::FieldVector< double, dim> > gradients(referenceGradients.size());
       for (size_t i=0; i< gradients.size(); i++)
         jacobian.mv(referenceGradients[i][0], gradients[i]);
       
@@ -658,7 +657,7 @@ public:
   {
     std::cout << "assemble Stiffness-Matrix begins." << std::endl;
 
-    MatrixIndexSet occupationPattern;
+    Dune::MatrixIndexSet occupationPattern;
     getOccupationPattern(occupationPattern);
     occupationPattern.exportIdx(matrix);
     matrix = 0;
@@ -679,7 +678,7 @@ public:
       const int localPhiOffset = localView.size();
   //     Dune::Timer Time;
       //std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
-      Matrix<FieldMatrix<double,1,1> > elementMatrix;
+      Dune::Matrix<Dune::FieldMatrix<double,1,1> > elementMatrix;
       // computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal);
       computeElementStiffnessMatrix(localView, elementMatrix);
       
@@ -792,7 +791,7 @@ public:
     // Output : determine all Component-Indices that correspond to that basis-function
     auto localView = basis_.localView();
 
-    FieldVector<int,3> row;
+    Dune::FieldVector<int,3> row;
     int elementCounter = 0;
 
     for (const auto& element : elements(basis_.gridView()))
@@ -823,7 +822,7 @@ public:
     unsigned int arbitraryLeafIndex =  parameterSet_.template get<unsigned int>("arbitraryLeafIndex", 0);
     unsigned int arbitraryElementNumber =  parameterSet_.template get<unsigned int>("arbitraryElementNumber", 0);
     //Determine 3 global indices (one for each component of an arbitrary local FE-function)
-    FieldVector<int,3> row = arbitraryComponentwiseIndices(arbitraryElementNumber,arbitraryLeafIndex);
+    Dune::FieldVector<int,3> row = arbitraryComponentwiseIndices(arbitraryElementNumber,arbitraryLeafIndex);
 
     for (int k = 0; k<dim; k++)
     {
@@ -876,12 +875,12 @@ public:
 
   auto integralMean(VectorCT& coeffVector)
   {
-    auto GVFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis_,coeffVector);
+    auto GVFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,dim>>(basis_,coeffVector);
     auto localfun = localFunction(GVFunction);
 
     auto localView = basis_.localView();
 
-    FieldVector<double,3> r = {0.0, 0.0, 0.0};
+    Dune::FieldVector<double,3> r = {0.0, 0.0, 0.0};
     double area = 0.0;
 
     // Compute Area integral & integral of FE-function
@@ -893,7 +892,7 @@ public:
       
   //     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; //TEST
       int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
-      const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
+      const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), orderQR);
 
       for(const auto& quadPoint : quad)
       {
@@ -968,20 +967,20 @@ public:
       if (Solvertype==1)  // CG - SOLVER
       {
           std::cout << "------------ CG - Solver ------------" << std::endl;
-          MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_);
+          Dune::MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_);
 
           // Sequential incomplete LU decomposition as the preconditioner
-          SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_,1.0);
+          Dune::SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_,1.0);
           int iter = parameterSet_.get<double>("iterations_CG", 999);
           // Preconditioned conjugate-gradient solver
-          CGSolver<VectorCT> solver(op,
+          Dune::CGSolver<VectorCT> solver(op,
                                   ilu0, //NULL,
                                   1e-8, // desired residual reduction factorlack
                                   iter,   // maximum number of iterations
                                   Solver_verbosity,
                                   true    // Try to estimate condition_number
                                   );   // verbosity of the solver
-          InverseOperatorResult statistics;
+          Dune::InverseOperatorResult statistics;
           std::cout << "solve linear system for x_1.\n";
           solver.apply(x_1_, load_alpha1_, statistics);
           std::cout << "solve linear system for x_2.\n";
@@ -999,13 +998,13 @@ public:
       {
           std::cout << "------------ GMRES - Solver ------------" << std::endl;
           // Turn the matrix into a linear operator
-          MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_);
+          Dune::MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_);
 
           // Fancy (but only) way to not have a preconditioner at all
-          Richardson<VectorCT,VectorCT> preconditioner(1.0);
+          Dune::Richardson<VectorCT,VectorCT> preconditioner(1.0);
 
           // Construct the iterative solver
-          RestartedGMResSolver<VectorCT> solver(
+          Dune::RestartedGMResSolver<VectorCT> solver(
               stiffnessOperator, // Operator to invert
               preconditioner,    // Preconditioner
               1e-10,             // Desired residual reduction factor
@@ -1015,7 +1014,7 @@ public:
               Solver_verbosity);                // Verbosity of the solver
 
           // Object storing some statistics about the solving process
-          InverseOperatorResult statistics;
+          Dune::InverseOperatorResult statistics;
 
           // solve for different Correctors (alpha = 1,2,3)
           solver.apply(x_1_, load_alpha1_, statistics);     //load_alpha1 now contains the corresponding residual!!
@@ -1033,9 +1032,9 @@ public:
           std::cout << "------------ QR - Solver ------------" << std::endl;
           log_ << "solveLinearSystems: We use QR solver.\n";
           //TODO install suitesparse
-          SPQR<MatrixCT> sPQR(stiffnessMatrix_);
+          Dune::SPQR<MatrixCT> sPQR(stiffnessMatrix_);
           sPQR.setVerbosity(1);
-          InverseOperatorResult statisticsQR;
+          Dune::InverseOperatorResult statisticsQR;
 
           sPQR.apply(x_1_, load_alpha1_, statisticsQR);
           std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
@@ -1198,16 +1197,16 @@ public:
       std::string vtkOutputName = parameterSet_.get("outputPath", "../../outputs") + "/CellProblem-result"; 
       std::cout << "vtkOutputName:" << vtkOutputName << std::endl;
 
-      VTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView());
+      Dune::VTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView());
       vtkWriter.addVertexData(
-          Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_1_),
-          VTK::FieldInfo("Corrector phi_1 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));     
+          Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_1_),
+          Dune::VTK::FieldInfo("Corrector phi_1 level"+ std::to_string(level) , Dune::VTK::FieldInfo::Type::vector, dim));
       vtkWriter.addVertexData(
-          Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_2_),
-          VTK::FieldInfo("Corrector phi_2 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
+          Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_2_),
+          Dune::VTK::FieldInfo("Corrector phi_2 level"+ std::to_string(level) , Dune::VTK::FieldInfo::Type::vector, dim));
       vtkWriter.addVertexData(
-          Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_3_),
-          VTK::FieldInfo("Corrector phi_3 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
+          Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_3_),
+          Dune::VTK::FieldInfo("Corrector phi_3 level"+ std::to_string(level) , Dune::VTK::FieldInfo::Type::vector, dim));
       vtkWriter.write(vtkOutputName  + "-level"+ std::to_string(level));
       std::cout << "wrote Corrector-VTK data to file: " + vtkOutputName + "-level" + std::to_string(level) << std::endl;
   }
@@ -1232,9 +1231,9 @@ public:
     double error_3 = 0.0;
 
     auto localView = basis_.localView();
-    auto GVFunc_1 = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_);
-    auto GVFunc_2 = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_);
-    auto GVFunc_3 = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_);
+    auto GVFunc_1 = Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_);
+    auto GVFunc_2 = Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_);
+    auto GVFunc_3 = Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_);
     auto localfun_1 = localFunction(GVFunc_1);
     auto localfun_2 = localFunction(GVFunc_2);
     auto localfun_3 = localFunction(GVFunc_3);
@@ -1248,7 +1247,7 @@ public:
       const auto& localFiniteElement = localView.tree().child(0).finiteElement();
 
       int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
-      const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
+      const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), orderQR);
 
       for(const auto& quadPoint : quad)
       {
@@ -1274,9 +1273,9 @@ public:
 
     auto localView = basis_.localView();
 
-    auto GVFunc_1 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_));
-    auto GVFunc_2 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_));
-    auto GVFunc_3 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_));
+    auto GVFunc_1 = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_));
+    auto GVFunc_2 = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_));
+    auto GVFunc_3 = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_));
     auto localfun_1 = localFunction(GVFunc_1);
     auto localfun_2 = localFunction(GVFunc_2);
     auto localfun_3 = localFunction(GVFunc_3);
@@ -1293,7 +1292,7 @@ public:
       const auto nSf = localFiniteElement.localBasis().size();
 
       int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 );  
-      const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+      const auto& quad = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
 
       for (const auto& quadPoint : quad)
       {
@@ -1325,9 +1324,9 @@ public:
 
     auto localView = basis_.localView();
 
-    auto GVFunc_1 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_));
-    auto GVFunc_2 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_));
-    auto GVFunc_3 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_));
+    auto GVFunc_1 = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_));
+    auto GVFunc_2 = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_));
+    auto GVFunc_3 = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_));
     auto localfun_1 = localFunction(GVFunc_1);
     auto localfun_2 = localFunction(GVFunc_2);
     auto localfun_3 = localFunction(GVFunc_3);
@@ -1387,7 +1386,7 @@ public:
       //     int orderQR = 1;
       //     int orderQR = 2;
       //     int orderQR = 3;
-          const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
+          const auto& quad = Dune::QuadratureRules<double, dim>::rule(e.type(), orderQR);
 
           for (const auto& quadPoint : quad) 
           {
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index 9329f1634fda83af82f0b28c2c896b99aa859505..0a281d4f232e0df9529b6ec3038ae706cc283239 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -9,7 +9,6 @@
 #include <dune/istl/matrix.hh>
 #include <dune/common/parametertree.hh>
 
-using namespace Dune;
 using namespace MatrixOperations;
 using std::shared_ptr;
 using std::make_shared;
@@ -91,7 +90,7 @@ public:
         // auto lambda_ = *correctorComputer_.getLambda();
         auto gamma = correctorComputer_.getGamma();
         auto basis = *correctorComputer_.getBasis();
-        ParameterTree parameterSet = correctorComputer_.getParameterSet();
+        Dune::ParameterTree parameterSet = correctorComputer_.getParameterSet();
 
 		shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_.getCorr_phi1(), 
                                             correctorComputer_.getCorr_phi2(),
@@ -114,7 +113,7 @@ public:
             double prestrain = 0.0;
             auto localView = basis.localView();
             // auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiContainer[a]));
-            auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[a]));
+            auto GVFunc_a = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[a]));
             //   auto GVFunc_b = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,phiContainer[b]));
             auto localfun_a = localFunction(GVFunc_a);
             //   auto localfun_b = localFunction(GVFunc_b);
@@ -157,7 +156,7 @@ public:
             //     int orderQR = 1;
             //     int orderQR = 2;
             //     int orderQR = 3;
-                const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
+                const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), orderQR);
 
                 for (const auto& quadPoint : quad) 
                 {
@@ -246,7 +245,7 @@ public:
         //writeMatrixToMatlab(Q, "../../Matlab-Programs/QMatrix.txt");
         writeMatrixToMatlab(Qeff_, outputPath + "/QMatrix.txt");
         // write effective Prestrain in Matrix for Output
-        FieldMatrix<double,1,3> BeffMat;
+        Dune::FieldMatrix<double,1,3> BeffMat;
         BeffMat[0] = Beff_;
         writeMatrixToMatlab(BeffMat, outputPath + "/BMatrix.txt");
         return;
@@ -288,7 +287,7 @@ public:
             const auto& localFiniteElement = localView.tree().child(0).finiteElement();
 
             int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
-            const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
+            const auto& quad = Dune::QuadratureRules<double, dim>::rule(e.type(), orderQR);
             for (const auto& quadPoint : quad) 
             {
                 const auto& quadPos = quadPoint.position();
@@ -431,4 +430,4 @@ public:
 
 
 
-#endif
\ No newline at end of file
+#endif
diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index 577a151cde9d8306b1071f129ac1832d2023c8ef..fc1670e5aed0a87f2acdb6e54f92b5c786609aa7 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -3,10 +3,8 @@
 
 namespace MatrixOperations {
 
-	using namespace Dune;
-
-	using MatrixRT = FieldMatrix< double, 3, 3>;
-	using VectorRT = FieldVector< double, 3>;	
+	using MatrixRT = Dune::FieldMatrix< double, 3, 3>;
+	using VectorRT = Dune::FieldVector< double, 3>;
 
 	using std::sin;
 	using std::cos;
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index d0fd35b0ca8571a7570338c340a573ef91658856..904b1e581898883caa65f8a1b025f146ead4ee05 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -12,10 +12,9 @@
 // #include <dune/microstructure/materialDefinitions.hh> //deprecated
 
 
-using namespace Dune;
 using namespace MatrixOperations;
-using namespace Functions;
-using namespace Functions::BasisFactory;
+using namespace Dune::Functions;
+using namespace Dune::Functions::BasisFactory;
 using std::pow;
 using std::abs;
 using std::sqrt;
@@ -87,9 +86,9 @@ public:
   using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
   using LocalDomain = typename GridView::template Codim<0>::Geometry::LocalCoordinate;
   using Element = typename GridViewEntitySet<GridView, 0>::Element;
-  using ScalarRT = FieldVector< double, 1>;
-  using VectorRT = FieldVector< double, dimworld>;
-  using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
+  using ScalarRT = Dune::FieldVector< double, 1>;
+  using VectorRT = Dune::FieldVector< double, dimworld>;
+  using MatrixRT = Dune::FieldMatrix< double, dimworld, dimworld>;
   using FuncScalar = std::function< double(const Domain&) >;
   using Func2int = std::function< int(const Domain&) >;
   using Func2Tensor = std::function< MatrixRT(const Domain&) >;
@@ -103,7 +102,7 @@ public:
 protected:
 
   const GridView& gridView_;
-  const ParameterTree& parameterSet_;
+  const Dune::ParameterTree& parameterSet_;
   std::string materialFunctionName_;
 
   // MatrixFunc L1_;
@@ -111,10 +110,10 @@ protected:
   // MatrixFunc L3_;
 
   // Elasticity tensors (in voigt notation)
-  FieldMatrix<double,6,6> L1_;
-  FieldMatrix<double,6,6> L2_;
-  FieldMatrix<double,6,6> L3_;
-  FieldMatrix<double,6,6> L4_; //not used yet
+  Dune::FieldMatrix<double,6,6> L1_;
+  Dune::FieldMatrix<double,6,6> L2_;
+  Dune::FieldMatrix<double,6,6> L3_;
+  Dune::FieldMatrix<double,6,6> L4_; //not used yet
 
   Func2Tensor prestrain1_;
   Func2Tensor prestrain2_; 
@@ -138,14 +137,14 @@ protected:
 
 
   //Transformation matrix for Voigt notation
-  FieldMatrix<double,6,6> D_ = {{1.0, 0.0, 0.0, 0.0 , 0.0, 0.0},
+  Dune::FieldMatrix<double,6,6> D_ = {{1.0, 0.0, 0.0, 0.0 , 0.0, 0.0},
                               {0.0, 1.0, 0.0, 0.0 , 0.0, 0.0},
                               {0.0, 0.0, 1.0, 0.0 , 0.0, 0.0},
                               {0.0, 0.0, 0.0, sqrt(2.0) , 0.0, 0.0},
                               {0.0, 0.0, 0.0, 0.0 , sqrt(2.0), 0.0},
                               {0.0, 0.0, 0.0, 0.0 , 0.0, sqrt(2.0)}
                             };
-  FieldMatrix<double,6,6> D_inv = {{1.0, 0.0, 0.0, 0.0 , 0.0, 0.0},
+  Dune::FieldMatrix<double,6,6> D_inv = {{1.0, 0.0, 0.0, 0.0 , 0.0, 0.0},
                                   {0.0, 1.0, 0.0, 0.0 , 0.0, 0.0},
                                   {0.0, 0.0, 1.0, 0.0 , 0.0, 0.0},
                                   {0.0, 0.0, 0.0, 1.0/sqrt(2.0) , 0.0, 0.0},
@@ -160,7 +159,7 @@ public:
     // constructor
     ///////////////////////////////
     prestrainedMaterial(const GridView gridView,
-                        const ParameterTree& parameterSet)  
+                        const Dune::ParameterTree& parameterSet)
                        : gridView_(gridView), 
                          parameterSet_(parameterSet)
     {
@@ -266,7 +265,7 @@ public:
 
 
 
-FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module module, int phase)
+Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module module, int phase)
 {
   std::string materialParameters_string;
   std::cout << "Setup phase "<< phase << " as " << phaseType << " material." << std::endl;
@@ -289,14 +288,14 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
       break;
     }
     default:
-    DUNE_THROW(Exception, "Phase number not implemented."); 
+    DUNE_THROW(Dune::Exception, "Phase number not implemented.");
         break;
   }
 
   if(phaseType == "isotropic")                           //TODO SetupPHASE (phase_type)
   {
-      FieldVector<double,2> materialParameters(0);
-      module.get(materialParameters_string).toC<FieldVector<double,2>>(materialParameters);
+      Dune::FieldVector<double,2> materialParameters(0);
+      module.get(materialParameters_string).toC<Dune::FieldVector<double,2>>(materialParameters);
       return isotropic(materialParameters[0],materialParameters[1]);
   }
   if(phaseType == "orthotropic")                           //TODO SetupPHASE (phase_type)
@@ -306,8 +305,8 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
       VectorRT e1 = {1.0,0.0,0.0};
       VectorRT e2 = {0.0,1.0,0.0};
       VectorRT e3 = {0.0,0.0,1.0};
-      FieldVector<double,9> materialParameters(0);
-      module.get(materialParameters_string).toC<FieldVector<double,9>>(materialParameters);
+      Dune::FieldVector<double,9> materialParameters(0);
+      module.get(materialParameters_string).toC<Dune::FieldVector<double,9>>(materialParameters);
       return orthotropic(e1,e2,e3,materialParameters);
   }
   if(phaseType == "transversely_isotropic")                           //TODO SetupPHASE (phase_type)
@@ -317,8 +316,8 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
       VectorRT e1 = {1.0,0.0,0.0};
       VectorRT e2 = {0.0,1.0,0.0};
       VectorRT e3 = {0.0,0.0,1.0};
-      FieldVector<double,5> materialParameters(0);
-      module.get(materialParameters_string).toC<FieldVector<double,5>>(materialParameters);
+      Dune::FieldVector<double,5> materialParameters(0);
+      module.get(materialParameters_string).toC<Dune::FieldVector<double,5>>(materialParameters);
       return transversely_isotropic(e1,e2,e3,materialParameters);
   }
   if(phaseType == "general_anisotropic")                           //TODO SetupPHASE (phase_type)
@@ -328,13 +327,13 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
       VectorRT e1 = {1.0,0.0,0.0};
       VectorRT e2 = {0.0,1.0,0.0};
       VectorRT e3 = {0.0,0.0,1.0};
-      FieldMatrix<double,6,6> materialParameters(0); //compliance matric
-      module.get(materialParameters_string).toC<FieldMatrix<double,6,6>>(materialParameters);
+      Dune::FieldMatrix<double,6,6> materialParameters(0); //compliance matric
+      module.get(materialParameters_string).toC<Dune::FieldMatrix<double,6,6>>(materialParameters);
       printmatrix(std::cout, materialParameters, "Compliance matrix used:", "--");
       return general_anisotropic(e1,e2,e3,materialParameters);
   }
   else
-    DUNE_THROW(Exception, " Unknown material class for phase "); 
+    DUNE_THROW(Dune::Exception, " Unknown material class for phase ");
 }
 
 
@@ -463,7 +462,7 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
   ////////////////////////////////////////////////////////
 
   //--- Definition of isotropic elasticity tensor (in voigt notation)
-  FieldMatrix<double,6,6> isotropic(const double& mu, const double& lambda)
+  Dune::FieldMatrix<double,6,6> isotropic(const double& mu, const double& lambda)
   {
       return {{lambda+2.0*mu, lambda       , lambda       , 0.0   , 0.0   , 0.0   },
               {lambda       , lambda+2.0*mu, lambda       , 0.0   , 0.0   , 0.0   },
@@ -478,10 +477,10 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
   //--- Definition of orthotropic elasticity tensor (in voigt notation)
   // - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,E3,G12,G23,G31,nu12,nu13,nu23}
   // - Output: compliance Matrix
-  FieldMatrix<double,6,6> orthotropic(const VectorRT& framevector1,
+  Dune::FieldMatrix<double,6,6> orthotropic(const VectorRT& framevector1,
                                       const VectorRT& framevector2,
                                       const VectorRT& framevector3,
-                                      const FieldVector<double,9>& p //elastic moduli {E1,E2,E3,G12,G23,G31,nu12,nu13,nu23}
+                                      const Dune::FieldVector<double,9>& p //elastic moduli {E1,E2,E3,G12,G23,G31,nu12,nu13,nu23}
                                       ) 
   {
     // (see https://en.wikipedia.org/wiki/Orthotropic_material)
@@ -500,7 +499,7 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
       auto nu_32 = (p[8]/p[1])*p[2]; 
 
       //compliance matrix
-      FieldMatrix<double,6,6> C = {{1.0/E_1      , -nu_21/E_2   , -nu_31/E_3   , 0.0     , 0.0     , 0.0   },
+      Dune::FieldMatrix<double,6,6> C = {{1.0/E_1      , -nu_21/E_2   , -nu_31/E_3   , 0.0     , 0.0     , 0.0   },
                                        {-nu_12/E_1   , 1.0/E_2      , -nu_32/E_3   , 0.0     , 0.0     , 0.0   },
                                        {-nu_13/E_1   , -nu_23/E_2   , 1.0/E_3      , 0.0     , 0.0     , 0.0   },
                                        {  0.0        ,  0.0         , 0.0          , 1.0/G_23, 0.0     , 0.0   },
@@ -522,10 +521,10 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
   //--- Definition of transversely isotropic elasticity tensor (in voigt notation)
   // - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,G12,nu12,nu23}
   // - Output: compliance Matrix
-  FieldMatrix<double,6,6> transversely_isotropic(const VectorRT& framevector1,
+  Dune::FieldMatrix<double,6,6> transversely_isotropic(const VectorRT& framevector1,
                                                  const VectorRT& framevector2,
                                                  const VectorRT& framevector3,
-                                                 const FieldVector<double,5>& p //elastic moduli {E1,E2,G12,nu12,nu23}
+                                                 const Dune::FieldVector<double,5>& p //elastic moduli {E1,E2,G12,nu12,nu23}
                                                  ) 
   {
     // (see https://en.wikipedia.org/wiki/Poisson%27s_ratio)
@@ -544,7 +543,7 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
       //other three poisson ratios are not independent 
       
       //---compliance matrix
-      FieldMatrix<double,6,6> C = {{1.0/E_1      , -nu_21/E_2   , -nu_31/E_3   , 0.0     , 0.0     , 0.0   },
+      Dune::FieldMatrix<double,6,6> C = {{1.0/E_1      , -nu_21/E_2   , -nu_31/E_3   , 0.0     , 0.0     , 0.0   },
                                    {-nu_12/E_1   , 1.0/E_2      , -nu_32/E_3   , 0.0     , 0.0     , 0.0   },
                                    {-nu_13/E_1   , -nu_23/E_2   , 1.0/E_3      , 0.0     , 0.0     , 0.0   },
                                    {  0.0        ,  0.0         , 0.0          , 1.0/G_23, 0.0     , 0.0   },
@@ -565,10 +564,10 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
   // --- Definition of transversely isotropic elasticity tensor (in voigt notation)
   // - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,G12,nu12,nu23}
   // - Output: inverse compliance Matrix
-  FieldMatrix<double,6,6> general_anisotropic(const VectorRT& framevector1,
+  Dune::FieldMatrix<double,6,6> general_anisotropic(const VectorRT& framevector1,
                                               const VectorRT& framevector2,
                                               const VectorRT& framevector3,
-                                              const FieldMatrix<double,6,6>& C //compliance matrix
+                                              const Dune::FieldMatrix<double,6,6>& C //compliance matrix
                                               ) 
   {
       //--- return elasticity tensor (in Voigt notation)
@@ -581,14 +580,14 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
 
 
   // --- Helpers
-  static FieldVector<double,6> MatrixToVoigt(const MatrixRT& G) 
+  static Dune::FieldVector<double,6> MatrixToVoigt(const MatrixRT& G)
   {
     // return {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
     return {G[0][0], G[1][1], G[2][2], 2.0*G[1][2], 2.0*G[0][2], 2.0*G[0][1]};
 
   }
 
-  static MatrixRT VoigtToMatrix(const FieldVector<double,6>& v) 
+  static MatrixRT VoigtToMatrix(const Dune::FieldVector<double,6>& v)
   {
     // return {{v[0], v[5], v[4]}, {v[5], v[1], v[3]}, {v[4], v[3], v[2]}};
     return {{v[0], v[5]/2.0, v[4]/2.0}, {v[5]/2.0, v[1], v[3]/2.0}, {v[4]/2.0, v[3]/2.0, v[2]}};
@@ -600,8 +599,8 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
   {
     // FieldVector<double,6> G_tmp = {G[0][0], G[1][1], G[2][2], sqrt(2.0)*G[1][2], sqrt(2.0)*G[0][2], sqrt(2.0)*G[0][1]};
     // FieldVector<double,6> G_tmp = {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
-    FieldVector<double,6> G_tmp = MatrixToVoigt(G);
-    FieldVector<double,6> out(0);
+    Dune::FieldVector<double,6> G_tmp = MatrixToVoigt(G);
+    Dune::FieldVector<double,6> out(0);
     // printvector(std::cout, G_tmp, "G_tmp", "--");
 
     if (phase == 1)
@@ -640,7 +639,7 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
     else if (phase ==3)
       return prestrain3_(x);
     else 
-      DUNE_THROW(Exception, "Phase number not implemented."); 
+      DUNE_THROW(Dune::Exception, "Phase number not implemented.");
      
   }
 
@@ -667,7 +666,7 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
     ///////////////////////////////
     // Getter
     ///////////////////////////////
-    ParameterTree getParameterSet() const {return parameterSet_;}
+    Dune::ParameterTree getParameterSet() const {return parameterSet_;}
     // Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
     // Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
     Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;}
@@ -687,32 +686,32 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
   void writeVTKMaterialFunctions(std::array<int, dim> nElements, const int level)
   {
         std::string outputPath = parameterSet_.get("outputPath", "../../outputs");
-        using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+        using VTKGridType = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
         // VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
         VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},nElements);
         
-        using GridViewVTK = VTKGridType::LeafGridView;
+        using GridViewVTK = typename VTKGridType::LeafGridView;
         const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
         
         auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
         auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
 
         std::vector<double> indicatorFunction_CoeffP0;
-        Functions::interpolate(scalarP0FeBasis, indicatorFunction_CoeffP0, indicatorFunction_);
-        auto indicatorFunction_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, indicatorFunction_CoeffP0);
+        Dune::Functions::interpolate(scalarP0FeBasis, indicatorFunction_CoeffP0, indicatorFunction_);
+        auto indicatorFunction_P0 = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, indicatorFunction_CoeffP0);
         
         std::vector<double> indicatorFunction_CoeffP1;
-        Functions::interpolate(scalarP1FeBasis, indicatorFunction_CoeffP1, indicatorFunction_);
-        auto indicatorFunction_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, indicatorFunction_CoeffP1);
+        Dune::Functions::interpolate(scalarP1FeBasis, indicatorFunction_CoeffP1, indicatorFunction_);
+        auto indicatorFunction_P1 = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, indicatorFunction_CoeffP1);
         
-        VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
+        Dune::VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
         
         MaterialVtkWriter.addCellData(
             indicatorFunction_P0,
-            VTK::FieldInfo("indicatorFunction_P0", VTK::FieldInfo::Type::scalar, 1));    
+            Dune::VTK::FieldInfo("indicatorFunction_P0", Dune::VTK::FieldInfo::Type::scalar, 1));
         MaterialVtkWriter.addVertexData(
             indicatorFunction_P1,
-            VTK::FieldInfo("indicatorFunction_P1", VTK::FieldInfo::Type::scalar, 1));    
+            Dune::VTK::FieldInfo("indicatorFunction_P1", Dune::VTK::FieldInfo::Type::scalar, 1));
         
         MaterialVtkWriter.write(outputPath + "/MaterialFunctions-level"+ std::to_string(level) );
         std::cout << "wrote data to file:" + outputPath +"/MaterialFunctions-level" + std::to_string(level) << std::endl;  
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index 9899928ebf8160562e8741fd0b7e37565869a10b..93e0ffc4078a18cbb5db3605bf2d0a3af67943e5 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -56,7 +56,6 @@
 
 // #include <dune/fufem/discretizationerror.hh>
 #include <dune/fufem/dunepython.hh>
-#include <python2.7/Python.h>
 
 // #include <dune/fufem/functions/virtualgridfunction.hh> //TEST