diff --git a/dumux/multidomain/common/multidomainassembler.hh b/dumux/multidomain/common/multidomainassembler.hh
index 2682ae3a945ae9aca3139f6803d3913add7ae922..3ee6eb5a91c7a55f14a49d9a19ce27b33d677574 100644
--- a/dumux/multidomain/common/multidomainassembler.hh
+++ b/dumux/multidomain/common/multidomainassembler.hh
@@ -29,22 +29,8 @@
 #include <dune/pdelab/constraints/constraintsparameters.hh>
 #include <dune/pdelab/multidomain/constraints.hh>
 
-
 namespace Dumux {
 
-/*!
- * \brief Prevents the setting of a Dirichlet constraint anywhere
- */
-struct NoDirichletConstraints :
-  public Dune::PDELab::DirichletConstraintsParameters
-{
-  template<typename I>
-  bool isDirichlet(const I& intersection, const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord) const
-  {
-    return false;
-  }
-};
-
 /*!
  * \brief An assembler for the global Jacobian matrix for multidomain models.
  */
@@ -84,13 +70,9 @@ class MultiDomainAssembler
     typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
     typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
 
-    typedef typename GET_PROP_TYPE(SubDomain1TypeTag, Constraints) Constraints1;
-    typedef typename GET_PROP_TYPE(SubDomain2TypeTag, Constraints) Constraints2;
-
     // copying the jacobian assembler is not a good idea
     MultiDomainAssembler(const MultiDomainAssembler &);
 
-
 public:
     //! \brief The constructor
     MultiDomainAssembler()
@@ -111,22 +93,13 @@ public:
         sdProblem1_ = &globalProblem_->sdProblem1();
         sdProblem2_ = &globalProblem_->sdProblem2();
 
-        fem1_ = Dune::make_shared<FEM1>();
-        fem2_ = Dune::make_shared<FEM2>();
-
-        constraints1_ = Dune::make_shared<Constraints1>();
-        constraints2_ = Dune::make_shared<Constraints2>();
+        fem1_ = Dune::make_shared<FEM1>(globalProblem_->sdGridView1());
+        fem2_ = Dune::make_shared<FEM2>(globalProblem_->sdGridView2());
 
         scalarGridFunctionSpace1_ = Dune::make_shared<ScalarGridFunctionSpace1>(globalProblem_->sdGridView1(),
-                                                                                *fem1_, *constraints1_);
+                                                                                *fem1_);
         scalarGridFunctionSpace2_ = Dune::make_shared<ScalarGridFunctionSpace2>(globalProblem_->sdGridView2(),
-                                                                                *fem2_, *constraints2_);
-        // constraints store indices of ghost dofs
-        constraints1_->compute_ghosts(*scalarGridFunctionSpace1_);
-        constraints2_->compute_ghosts(*scalarGridFunctionSpace2_);
-        Valgrind::CheckDefined(*constraints1_);
-        Valgrind::CheckDefined(*constraints2_);
-//        std::cerr  << __FILE__ << ":" << __LINE__ << "\n";
+                                                                                *fem2_);
 
         gridFunctionSpace1_ = Dune::make_shared<GridFunctionSpace1>(*scalarGridFunctionSpace1_);
         gridFunctionSpace2_ = Dune::make_shared<GridFunctionSpace2>(*scalarGridFunctionSpace2_);
@@ -147,17 +120,8 @@ public:
         couplingLocalOperator_ = Dune::make_shared<MultiDomainCouplingLocalOperator>(*globalProblem_);
         mdCoupling_ = Dune::make_shared<MultiDomainCoupling>(*mdSubProblem1_, *mdSubProblem2_, *couplingLocalOperator_);
 
-        // TODO proper constraints stuff
         constraintsTrafo_ = Dune::make_shared<MultiDomainConstraintsTrafo>();
 
-        NoDirichletConstraints dirichletVal;
-        auto constraints = Dune::PDELab::MultiDomain::constraints<Scalar>(*mdGridFunctionSpace_,
-                                                                          Dune::PDELab::MultiDomain::constrainSubProblem(*mdSubProblem1_,
-                                                                                                                         dirichletVal),
-                                                                          Dune::PDELab::MultiDomain::constrainSubProblem(*mdSubProblem2_,
-                                                                                                                         dirichletVal));
-        constraints.assemble(*constraintsTrafo_);
-
         mdGridOperator_ = Dune::make_shared<MultiDomainGridOperator>(*mdGridFunctionSpace_, *mdGridFunctionSpace_,
                                                                      *constraintsTrafo_, *constraintsTrafo_,
                                                                      *mdSubProblem1_, *mdSubProblem2_, *mdCoupling_);
@@ -165,25 +129,19 @@ public:
         matrix_ = Dune::make_shared<JacobianMatrix>(*mdGridOperator_);
         *matrix_ = 0;
 
-        residual_.resize(matrix_->N());
+        residual_ = Dune::make_shared<SolutionVector>(*mdGridFunctionSpace_);
     }
 
     //! \copydoc ImplicitAssembler::assemble()
     void assemble()
     {
-//        std::cerr  << __FILE__ << ":" << __LINE__ << "\n";
-
     	// assemble the matrix
     	*matrix_ = 0;
-
-        residual_ = 0;
     	mdGridOperator_->jacobian(globalProblem_->model().curSol(), *matrix_);
-//    	printmatrix(std::cout, matrix_->base(), "global stiffness matrix", "row", 11, 3);
 
     	// calculate the global residual
-    	residual_ = 0;
-    	mdGridOperator_->residual(globalProblem_->model().curSol(), residual_);
-//    	printvector(std::cout, residual_, "residual", "row", 200, 1, 3);
+    	*residual_ = 0;
+    	mdGridOperator_->residual(globalProblem_->model().curSol(), *residual_);
     }
 
     //! \copydoc ImplicitAssembler::reassembleAll()
@@ -205,9 +163,9 @@ public:
 
     //! \copydoc ImplicitAssembler::residual()
     const SolutionVector &residual() const
-    { return residual_; }
+    { return *residual_; }
     SolutionVector &residual()
-    { return residual_; }
+    { return *residual_; }
 
     /*!
      * \brief Return constant reference to the multidomain gridfunctionspace
@@ -235,9 +193,6 @@ private:
     Dune::shared_ptr<FEM1> fem1_;
     Dune::shared_ptr<FEM2> fem2_;
 
-    Dune::shared_ptr<Constraints1> constraints1_;
-    Dune::shared_ptr<Constraints2> constraints2_;
-
     Dune::shared_ptr<ScalarGridFunctionSpace1> scalarGridFunctionSpace1_;
     Dune::shared_ptr<ScalarGridFunctionSpace2> scalarGridFunctionSpace2_;
 
@@ -262,7 +217,7 @@ private:
 
     Dune::shared_ptr<JacobianMatrix> matrix_;
 
-    SolutionVector residual_;
+    Dune::shared_ptr<SolutionVector> residual_;
 };
 
 } // namespace Dumux
diff --git a/dumux/multidomain/common/multidomainmodel.hh b/dumux/multidomain/common/multidomainmodel.hh
index da236f8b502fda124d715bfe2ad1e6a9b4810798..6001cfda446239e88920b818c93bf684dace3adc 100644
--- a/dumux/multidomain/common/multidomainmodel.hh
+++ b/dumux/multidomain/common/multidomainmodel.hh
@@ -86,18 +86,18 @@ public:
         jacAsm_ = new JacobianAssembler();
         jacAsm_->init(problem);
 
-        uCur_.resize(jacAsm_->residual().size());
-        uPrev_.resize(jacAsm_->residual().size());
+        uCur_ = Dune::make_shared<SolutionVector>(jacAsm_->gridFunctionSpace());
+        uPrev_ = Dune::make_shared<SolutionVector>(jacAsm_->gridFunctionSpace());
 
-        uCur_= 0;
-        uPrev_= 0;
+        *uCur_= 0;
+        *uPrev_= 0;
 
         SplitAndMerge::mergeSolVectors(sdModel1().curSol(),
                                 	   sdModel2().curSol(),
-                                	   uCur_);
+                                	   *uCur_);
         SplitAndMerge::mergeSolVectors(sdModel1().prevSol(),
                                 	   sdModel2().prevSol(),
-                                	   uPrev_);
+                                	   *uPrev_);
     }
 
     /*!
@@ -132,25 +132,25 @@ public:
      * \brief Reference to the current solution as a block vector.
      */
     const SolutionVector &curSol() const
-    { return uCur_; }
+    { return *uCur_; }
 
     /*!
      * \brief Reference to the current solution as a block vector.
      */
     SolutionVector &curSol()
-    { return uCur_; }
+    { return *uCur_; }
 
     /*!
      * \brief Reference to the previous solution as a block vector.
      */
     const SolutionVector &prevSol() const
-    { return uPrev_; }
+    { return *uPrev_; }
 
     /*!
      * \brief Reference to the previous solution as a block vector.
      */
     SolutionVector &prevSol()
-    { return uPrev_; }
+    { return *uPrev_; }
 
     /*!
      * \brief Returns the operator assembler for the global jacobian of
@@ -253,8 +253,8 @@ public:
                 NewtonController &controller)
     {
 #if HAVE_VALGRIND
-        for (size_t i = 0; i < curSol().size(); ++i)
-            Valgrind::CheckDefined(curSol()[i]);
+        for (size_t i = 0; i < curSol().base().size(); ++i)
+            Valgrind::CheckDefined(curSol().base()[i]);
 #endif // HAVE_VALGRIND
 
         asImp_().updateBegin();
@@ -266,8 +266,8 @@ public:
             asImp_().updateSuccessful();
 
 #if HAVE_VALGRIND
-        for (size_t i = 0; i < curSol().size(); ++i) {
-            Valgrind::CheckDefined(curSol()[i]);
+        for (size_t i = 0; i < curSol().base().size(); ++i) {
+            Valgrind::CheckDefined(curSol().base()[i]);
         }
 #endif // HAVE_VALGRIND
 
@@ -285,7 +285,7 @@ public:
         sdModel1().updateBegin();
         sdModel2().updateBegin();
 
-        SplitAndMerge::mergeSolVectors(sdModel1().curSol(), sdModel2().curSol(), uCur_);
+        SplitAndMerge::mergeSolVectors(sdModel1().curSol(), sdModel2().curSol(), *uCur_);
     }
 
     //! \copydoc Dumux::ImplicitModel::updateSuccessful()
@@ -307,8 +307,8 @@ public:
     void advanceTimeLevel()
     {
         // merge the two sub-vectors together
-        SplitAndMerge::mergeSolVectors(sdModel1().curSol(), sdModel2().curSol(), uCur_);
-        SplitAndMerge::mergeSolVectors(sdModel1().prevSol(), sdModel2().prevSol(), uPrev_);
+        SplitAndMerge::mergeSolVectors(sdModel1().curSol(), sdModel2().curSol(), *uCur_);
+        SplitAndMerge::mergeSolVectors(sdModel1().prevSol(), sdModel2().prevSol(), *uPrev_);
     };
 
     //! \copydoc Dumux::ImplicitModel::updateFailed()
@@ -318,7 +318,7 @@ public:
         sdModel2().updateFailed();
 
         // merge the two sub-vectors together
-        SplitAndMerge::mergeSolVectors(sdModel1().curSol(), sdModel2().curSol(), uCur_);
+        SplitAndMerge::mergeSolVectors(sdModel1().curSol(), sdModel2().curSol(), *uCur_);
     };
 
      /*!
@@ -332,7 +332,7 @@ public:
         sdModel2().updateFailedTry();
 
         // merge the two sub-vectors together
-        SplitAndMerge::mergeSolVectors(sdModel1().curSol(), sdModel2().curSol(), uCur_);
+        SplitAndMerge::mergeSolVectors(sdModel1().curSol(), sdModel2().curSol(), *uCur_);
     };
 
     /*!
@@ -395,8 +395,8 @@ protected:
 
     // cur is the current solution, prev the solution of the previous
     // time step
-    SolutionVector uCur_;
-    SolutionVector uPrev_;
+    Dune::shared_ptr<SolutionVector> uCur_;
+    Dune::shared_ptr<SolutionVector> uPrev_;
 
     bool wasRestarted_;
 };
diff --git a/dumux/multidomain/common/multidomainnewtoncontroller.hh b/dumux/multidomain/common/multidomainnewtoncontroller.hh
index fafebcfba636ca002bb7fc2ae1775997deccae79..0b6b94683db04090469b7495ac5e89519d57cc5b 100644
--- a/dumux/multidomain/common/multidomainnewtoncontroller.hh
+++ b/dumux/multidomain/common/multidomainnewtoncontroller.hh
@@ -119,10 +119,10 @@ public:
         SolutionVector uNewI = uLastIter;
         uNewI -= deltaU;
 
-        for (unsigned int i = 0; i < uLastIter.size(); ++i) {
-            for (unsigned int j = 0; j < uLastIter[i].size(); ++j) {
-                Scalar vertexError = std::abs(deltaU[i][j]);
-                vertexError /= std::max<Scalar>(1.0, std::abs(uLastIter[i][j] + uNewI[i][j])/2);
+        for (unsigned int i = 0; i < uLastIter.base().size(); ++i) {
+            for (unsigned int j = 0; j < uLastIter.base()[i].size(); ++j) {
+                Scalar vertexError = std::abs(deltaU.base()[i][j]);
+                vertexError /= std::max<Scalar>(1.0, std::abs(uLastIter.base()[i][j] + uNewI.base()[i][j])/2);
 
                 this->error_ = std::max(this->error_, vertexError);
             }
@@ -151,7 +151,7 @@ public:
         // close to the final value, a reduction of 6 orders of
         // magnitude in the defect should be sufficient...
         try {
-            int converged = linearSolver_.solve(A, x, b);
+            int converged = linearSolver_.solve(A.base(), x.base(), b.base());
 
 #if HAVE_MPI
             // make sure all processes converged
@@ -185,7 +185,7 @@ public:
             Dumux::NumericalProblem p;
             std::string msg;
             std::ostringstream ms(msg);
-            ms << e.what() << "M=" << A[e.r][e.c];
+            ms << e.what() << "M=" << A.base()[e.r][e.c];
             p.message(ms.str());
             throw p;
         }
diff --git a/dumux/multidomain/common/multidomainpropertydefaults.hh b/dumux/multidomain/common/multidomainpropertydefaults.hh
index ae8bdb8bd29bf6d4a1d077e3d50027d1f3222576..54eeb13ae934ab9f51c090572e45da552670d048 100644
--- a/dumux/multidomain/common/multidomainpropertydefaults.hh
+++ b/dumux/multidomain/common/multidomainpropertydefaults.hh
@@ -36,8 +36,6 @@
 #include <dune/pdelab/multidomain/coupling.hh>
 #include <dune/pdelab/multidomain/gridoperator.hh>
 
-#include <dune/pdelab/finiteelementmap/p1fem.hh>
-#include <dune/pdelab/finiteelementmap/q1fem.hh>
 #include <dune/pdelab/gridoperator/gridoperator.hh>
 
 #include "subdomainpropertydefaults.hh"
@@ -71,15 +69,17 @@ SET_PROP(MultiDomain, MultiDomainGridFunctionSpace)
 {
 private:
     typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid;
+    typedef typename Dune::PDELab::LexicographicOrderingTag OrderingTag;
     typedef typename GET_PROP_TYPE(TypeTag, SubDomain1TypeTag) SubTypeTag1;
     typedef typename GET_PROP_TYPE(TypeTag, SubDomain2TypeTag) SubTypeTag2;
     typedef typename GET_PROP_TYPE(SubTypeTag1, GridFunctionSpace) GridFunctionSpace1;
     typedef typename GET_PROP_TYPE(SubTypeTag2, GridFunctionSpace) GridFunctionSpace2;
 public:
     typedef Dune::PDELab::MultiDomain::MultiDomainGridFunctionSpace<MDGrid,
-																	Dune::PDELab::ISTLVectorBackend<1>,
-																	GridFunctionSpace1,
-																	GridFunctionSpace2> type;
+                                                                    Dune::PDELab::ISTLVectorBackend<>,
+                                                                    OrderingTag,
+                                                                    GridFunctionSpace1,
+                                                                    GridFunctionSpace2> type;
 };
 
 // set the subdomain equality condition by default
@@ -140,16 +140,17 @@ public:
 SET_PROP(MultiDomain, MultiDomainConstraintsTrafo)
 {
 private:
-    typedef typename GET_PROP_TYPE(TypeTag, SubDomain1TypeTag) SubTypeTag1;
+    typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGridFunctionSpace) MDGridFunctionSpace;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
 public:
-    typedef typename GET_PROP_TYPE(SubTypeTag1, ConstraintsTrafo) type;
+    typedef typename MDGridFunctionSpace::template ConstraintsContainer<Scalar>::Type type;
 };
 
 SET_PROP(MultiDomain, MultiDomainGridOperator)
 {
 private:
     typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGridFunctionSpace) MDGridFunctionSpace;
-    typedef Dune::PDELab::ISTLBCRSMatrixBackend<1,1> MBE;
+    typedef Dune::PDELab::ISTLMatrixBackend MBE;
     typedef typename GET_PROP_TYPE(TypeTag, MultiDomainSubProblem1) MDSubProblem1;
     typedef typename GET_PROP_TYPE(TypeTag, MultiDomainSubProblem2) MDSubProblem2;
     typedef typename GET_PROP_TYPE(TypeTag, MultiDomainCoupling) MDCoupling;
diff --git a/dumux/multidomain/common/splitandmerge.hh b/dumux/multidomain/common/splitandmerge.hh
index 8e5d32ffb399b8967739a9a743cffa5af7cce83e..91fe462a6f7d4289fd93b1dcaf36d753a18ab1fd 100644
--- a/dumux/multidomain/common/splitandmerge.hh
+++ b/dumux/multidomain/common/splitandmerge.hh
@@ -65,8 +65,8 @@ public:
                                 const SolutionVector2 &vec2,
                                 SolutionVector &dest)
     {
-//        printvector(std::cout, vec1, "vec1", "row", 200, 1, 3);
-//        printvector(std::cout, vec2, "vec2", "row", 200, 1, 3);
+//         printvector(std::cout, vec1, "vec1", "row", 200, 1, 3);
+//         printvector(std::cout, vec2, "vec2", "row", 200, 1, 3);
 
         int nDofs1 = vec1.size();
         int nDofs2 = vec2.size();
@@ -74,19 +74,19 @@ public:
         for (int i = 0; i < nDofs1; ++i)
             for (int j = 0; j < numEq1; j++)
             {
-				dest[i*numEq1 + j][0] = vec1[i][j];
+				dest.base()[i*numEq1 + j][0] = vec1[i][j];
 
-                Valgrind::CheckDefined(dest[i*numEq1 + j][0]);
+                Valgrind::CheckDefined(dest.base()[i*numEq1 + j][0]);
             }
 
         for (int i = 0; i < nDofs2; ++i)
             for (int j = 0; j < numEq2; j++)
             {
-				dest[nDofs1*numEq1 + i*numEq2 + j][0] = vec2[i][j];
+				dest.base()[nDofs1*numEq1 + i*numEq2 + j][0] = vec2[i][j];
 
-                Valgrind::CheckDefined(dest[nDofs1*numEq1 + i*numEq2 + j][0]);
+                Valgrind::CheckDefined(dest.base()[nDofs1*numEq1 + i*numEq2 + j][0]);
             }
-//        printvector(std::cout, dest, "dest", "row", 200, 1, 3);
+//         printvector(std::cout, dest.base(), "dest", "row", 200, 1, 3);
     }
 
     /*!
@@ -102,21 +102,21 @@ public:
                                SolutionVector1 &dest1,
                                SolutionVector2 &dest2)
     {
-//        printvector(std::cout, vec, "vec", "row", 200, 1, 3);
+//         printvector(std::cout, vec.base(), "vec", "row", 200, 1, 3);
 
         int nDofs1 = dest1.size();
         int nDofs2 = dest2.size();
 
         for (int i = 0; i < nDofs1; ++i)
             for (int j = 0; j < numEq1; j++)
-				dest1[i][j] = vec[i*numEq1 + j][0];
+				dest1[i][j] = vec.base()[i*numEq1 + j][0];
 
         for (int i = 0; i < nDofs2; ++i)
             for (int j = 0; j < numEq2; j++)
-				dest2[i][j] = vec[nDofs1*numEq1 + i*numEq2 + j][0];
+				dest2[i][j] = vec.base()[nDofs1*numEq1 + i*numEq2 + j][0];
 
-//        printvector(std::cout, dest1, "dest1", "row", 200, 1, 3);
-//        printvector(std::cout, dest2, "dest2", "row", 200, 1, 3);
+//         printvector(std::cout, dest1, "dest1", "row", 200, 1, 3);
+//         printvector(std::cout, dest2, "dest2", "row", 200, 1, 3);
     }
 
     /*!
@@ -136,19 +136,19 @@ public:
     {
         int nDofs1 = vec1.size();
         int nDofs2 = vec2.size();
-//        printvector(std::cout, vec1, "vec1", "row", 200, 1, 3);
-//        printvector(std::cout, vec2, "vec2", "row", 200, 1, 3);
+//         printvector(std::cout, vec1, "vec1", "row", 200, 1, 3);
+//         printvector(std::cout, vec2, "vec2", "row", 200, 1, 3);
 
         for (int i = 0; i < nDofs1; ++i)
             for (int j = 0; j < numEq1; j++)
-                dest[i*numEq1 + j] = vec1[i][j];
+                dest.base()[i*numEq1 + j] = vec1[i][j];
 
         for (int i = 0; i < nDofs2; ++i)
         {
             for (int j = numEq1; j < numEq2; j++)
-                dest[nDofs1*numEq1 + i*(numEq2-numEq1) + (j - numEq1)] = vec2[i][j];
+                dest.base()[nDofs1*numEq1 + i*(numEq2-numEq1) + (j - numEq1)] = vec2[i][j];
         }
-//        printvector(std::cout, dest, "dest", "row", 200, 1, 3);
+//         printvector(std::cout, dest.base(), "dest", "row", 200, 1, 3);
     }
 
     /*!
@@ -169,24 +169,24 @@ public:
         int nDofs1 = dest1.size();
         int nDofs2 = dest2.size();
 
-//        printvector(std::cout, vec, "vec", "row", 200, 1, 3);
+//         printvector(std::cout, vec.base(), "vec", "row", 200, 1, 3);
         for (int i = 0; i < nDofs1; ++i)
             for (int j = 0; j < numEq1; j++)
-                dest1[i][j] = vec[i*numEq1 + j];
+                dest1[i][j] = vec.base()[i*numEq1 + j];
 
         for (int i = 0; i < nDofs2; ++i)
         {
             int blockIdxCoupled = subDOFToCoupledDOF[i];
             for (int j = 0; j < numEq1; j++)
             {
-                dest2[i][j] = vec[blockIdxCoupled*numEq1 + j];
+                dest2[i][j] = vec.base()[blockIdxCoupled*numEq1 + j];
             }
 
             for (int j = numEq1; j < numEq2; j++)
-                dest2[i][j] = vec[nDofs1*numEq1 + i*(numEq2-numEq1) + (j - numEq1)];
+                dest2[i][j] = vec.base()[nDofs1*numEq1 + i*(numEq2-numEq1) + (j - numEq1)];
         }
-//        printvector(std::cout, dest1, "dest1", "row", 200, 1, 3);
-//        printvector(std::cout, dest2, "dest2", "row", 200, 1, 3);
+//         printvector(std::cout, dest1, "dest1", "row", 200, 1, 3);
+//         printvector(std::cout, dest2, "dest2", "row", 200, 1, 3);
     }
 
 
diff --git a/dumux/multidomain/common/subdomainproperties.hh b/dumux/multidomain/common/subdomainproperties.hh
index 5f59a361173fbc3fcd1e62070cbb523c2ac88362..38b5915f2ca24964726ff4ad258e945c38b49b54 100644
--- a/dumux/multidomain/common/subdomainproperties.hh
+++ b/dumux/multidomain/common/subdomainproperties.hh
@@ -31,7 +31,7 @@ namespace Properties
 // \{
 
 //////////////////////////////////////////////////////////////////
-// Type tags tags
+// Type tags
 //////////////////////////////////////////////////////////////////
 
 //! The type tag from which sub-problems of coupling models inherit
@@ -49,18 +49,9 @@ NEW_PROP_TAG(ScalarGridFunctionSpace);
 //! Specifies the grid function space used for sub-problems
 NEW_PROP_TAG(GridFunctionSpace);
 
-//! Specifies the grid operator used for sub-problems
-NEW_PROP_TAG(GridOperator);
-
-//! Specifies the grid operator space used for sub-problems
-NEW_PROP_TAG(GridOperatorSpace);
-
 //! Specifies the type of the constraints
 NEW_PROP_TAG(Constraints);
 
-//! Specifies the type of the constraints transformation
-NEW_PROP_TAG(ConstraintsTrafo);
-
 //! Specifies the local finite element space
 NEW_PROP_TAG(LocalFEMSpace);
 
diff --git a/dumux/multidomain/common/subdomainpropertydefaults.hh b/dumux/multidomain/common/subdomainpropertydefaults.hh
index 82d55af9381ed412823f54e0523efa9cf7694d18..22f13955585c3b05d8dd83c24a1200b077bdd301 100644
--- a/dumux/multidomain/common/subdomainpropertydefaults.hh
+++ b/dumux/multidomain/common/subdomainpropertydefaults.hh
@@ -25,8 +25,9 @@
 
 #include <dune/grid/multidomaingrid.hh>
 #include <dune/pdelab/backend/istlvectorbackend.hh>
-#include <dune/pdelab/finiteelementmap/q1fem.hh>
-#include <dune/pdelab/finiteelementmap/conformingconstraints.hh>
+#include <dune/pdelab/backend/istlmatrixbackend.hh>
+#include <dune/pdelab/finiteelementmap/qkfem.hh>
+#include <dune/pdelab/constraints/conforming.hh>
 
 #include "subdomainproperties.hh"
 #include "multidomainproperties.hh"
@@ -66,6 +67,9 @@ public:
     typedef typename GET_PROP_TYPE(MultiDomainTypeTag, TimeManager) type;
 };
 
+// set the constraints for the sub-models
+SET_TYPE_PROP(SubDomain, Constraints, Dune::PDELab::NoConstraints);
+
 // set the grid functions space for the sub-models
 SET_PROP(SubDomain, ScalarGridFunctionSpace)
 {
@@ -76,7 +80,7 @@ SET_PROP(SubDomain, ScalarGridFunctionSpace)
     enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
  public:
     typedef Dune::PDELab::GridFunctionSpace<GridView, FEM, Constraints,
-        Dune::PDELab::ISTLVectorBackend<1> > type;
+        Dune::PDELab::ISTLVectorBackend<> > type;
 };
 
 // set the grid functions space for the sub-models
@@ -84,41 +88,10 @@ SET_PROP(SubDomain, GridFunctionSpace)
 {private:
     typedef typename GET_PROP_TYPE(TypeTag, ScalarGridFunctionSpace) ScalarGridFunctionSpace;
     enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
+    typedef typename Dune::PDELab::EntityBlockedOrderingTag OrderingTag;
+    typedef typename Dune::PDELab::ISTLVectorBackend<> VBE;
 public:
-    typedef Dune::PDELab::PowerGridFunctionSpace<ScalarGridFunctionSpace, numEq, Dune::PDELab::GridFunctionSpaceBlockwiseMapper> type;
-};
-
-// set the grid function space for the sub-models
-SET_TYPE_PROP(SubDomain, Constraints, Dune::PDELab::NoConstraints);
-
-// set the grid functions space for the sub-models
-SET_PROP(SubDomain, ConstraintsTrafo)
-{private:
-  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-  typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
-public:
-    typedef typename GridFunctionSpace::template ConstraintsContainer<Scalar>::Type type;
-};
-
-// set the grid operator space used for submodels
-// DEPRECATED: use GridOperator instead
-SET_PROP(SubDomain, GridOperatorSpace)
-{
-private:
-    typedef typename GET_PROP_TYPE(TypeTag, ConstraintsTrafo) ConstraintsTrafo;
-    typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
-    typedef typename GET_PROP_TYPE(TypeTag, LocalOperator) LocalOperator;
-    enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
-
-public:
-    typedef Dune::PDELab::GridOperatorSpace<GridFunctionSpace,
-        GridFunctionSpace,
-        LocalOperator,
-        ConstraintsTrafo,
-        ConstraintsTrafo,
-        Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
-        true
-        > type;
+    typedef Dune::PDELab::PowerGridFunctionSpace<ScalarGridFunctionSpace, numEq, VBE, OrderingTag> type;
 };
 
 // use the local FEM space associated with cubes by default
@@ -129,7 +102,7 @@ SET_PROP(SubDomain, LocalFEMSpace)
     enum{dim = GridView::dimension};
 
 public:
-    typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim>  type;
+    typedef Dune::PDELab::QkLocalFiniteElementMap<GridView,Scalar,Scalar,1>  type;
 };
 
 SET_PROP(SubDomain, ParameterTree)
diff --git a/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh b/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
index ad44c24a743b54aef2f2ab41c54d6496e8c04628..36b985ce18e683ea8636a813d130c0124048d03d 100644
--- a/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
@@ -29,9 +29,6 @@
 #include "config.h"
 #endif
 
-#include <dune/pdelab/finiteelementmap/conformingconstraints.hh>
-#include <dune/pdelab/gridoperator/gridoperator.hh>
-
 #include <dumux/implicit/common/implicitporousmediaproblem.hh>
 #include <dumux/implicit/2p2cni/2p2cnimodel.hh>
 #include <dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh>
@@ -67,35 +64,6 @@ SET_PROP(TwoPTwoCNISubProblem, ReplaceCompEqIdx)
     static const int value = Indices::contiNEqIdx;
 };
 
-// Set the model parameter group for the TypeTag
-//SET_PROP(TwoPTwoCNISubProblem, ModelParameterGroup)
-//{ static const char *value() { return "PM"; }; };
-
-// set the constraints for multidomain / pdelab
-SET_TYPE_PROP(TwoPTwoCNISubProblem, Constraints,
-        Dune::PDELab::NonoverlappingConformingDirichletConstraints);
-
-// set the grid operator
-SET_PROP(TwoPTwoCNISubProblem, GridOperator)
-{
-private:
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, ConstraintsTrafo) ConstraintsTrafo;
-    typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
-    typedef Dumux::PDELab::MultiDomainLocalOperator<TypeTag> LocalOperator;
-
-    enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
-
-public:
-    typedef Dune::PDELab::GridOperator<GridFunctionSpace,
-            GridFunctionSpace, LocalOperator,
-            Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
-            Scalar, Scalar, Scalar,
-            ConstraintsTrafo,
-            ConstraintsTrafo,
-            true> type;
-};
-
 // the fluidsystem is set in the coupled problem
 SET_PROP(TwoPTwoCNISubProblem, FluidSystem)
 {
diff --git a/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh b/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
index fb6de08f01c2c95dad93130526ee62a5303e7bbe..9dfc2c0e48973e78925dbd406b8f69a2a1b55ee6 100644
--- a/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
@@ -25,9 +25,6 @@
 #ifndef DUMUX_STOKES2CNI_SUBPROBLEM_HH
 #define DUMUX_STOKES2CNI_SUBPROBLEM_HH
 
-#include <dune/pdelab/finiteelementmap/conformingconstraints.hh>
-#include <dune/pdelab/gridoperator/gridoperator.hh>
-
 #include <dumux/freeflow/stokesncni/stokesncnimodel.hh>
 #include <dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh>
 #include <dumux/multidomain/common/subdomainpropertydefaults.hh>
@@ -66,28 +63,6 @@ public:
     typedef typename MaterialLaw::Params type;
 };
 
-SET_TYPE_PROP(Stokes2cniSubProblem, Constraints,
-		Dune::PDELab::NonoverlappingConformingDirichletConstraints);
-
-// set the grid operator
-SET_PROP(Stokes2cniSubProblem, GridOperator)
-{
-private:
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, ConstraintsTrafo) ConstraintsTrafo;
-    typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
-    typedef Dumux::PDELab::MultiDomainLocalOperator<TypeTag> LocalOperator;
-
-    enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
-public:
-    typedef Dune::PDELab::GridOperator<GridFunctionSpace,
-            GridFunctionSpace, LocalOperator,
-            Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
-            Scalar, Scalar, Scalar,
-            ConstraintsTrafo, ConstraintsTrafo,
-            true> type;
-};
-
 SET_PROP(Stokes2cniSubProblem, FluidSystem)
 {
 private:
diff --git a/test/multidomain/2cstokes2p2c/2p2csubproblem.hh b/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
index 0a4d7a4f1e034e29ffadbf7d9d9b82e2d4999010..544edd41e37c8f66e00937ec4ab9fd5baf5ef6c4 100644
--- a/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
+++ b/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
@@ -25,8 +25,6 @@
 #ifndef DUMUX_2P2C_SUBPROBLEM_HH
 #define DUMUX_2P2C_SUBPROBLEM_HH
 
-#include <dune/pdelab/gridoperator/gridoperator.hh>
-
 #include <dumux/implicit/2p2c/2p2cindices.hh>
 #include <dumux/implicit/common/implicitporousmediaproblem.hh>
 #include <dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh>
@@ -62,31 +60,6 @@ SET_PROP(TwoPTwoCSubProblem, ReplaceCompEqIdx)
     static const int value = Indices::contiNEqIdx;
 };
 
-// set the constraints for multidomain / pdelab // not required here
-SET_TYPE_PROP(TwoPTwoCSubProblem, Constraints,
-        Dune::PDELab::NonoverlappingConformingDirichletConstraints);
-
-// set the grid operator
-SET_PROP(TwoPTwoCSubProblem, GridOperator)
-{
- private:
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, ConstraintsTrafo) ConstraintsTrafo;
-    typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
-    typedef Dumux::PDELab::MultiDomainLocalOperator<TypeTag> LocalOperator;
-
-    enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
-
- public:
-    typedef Dune::PDELab::GridOperator<GridFunctionSpace,
-        GridFunctionSpace, LocalOperator,
-        Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
-        Scalar, Scalar, Scalar,
-        ConstraintsTrafo,
-        ConstraintsTrafo,
-        true> type;
-};
-
 // the fluidsystem is set in the coupled problem
 SET_PROP(TwoPTwoCSubProblem, FluidSystem)
 {
diff --git a/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh b/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
index f5bfe913c96eef04398ec93ae6796c3657d33911..73713a2d71a952a784c5372b8874d06915579fa3 100644
--- a/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
+++ b/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
@@ -25,8 +25,6 @@
 #ifndef DUMUX_STOKES2C_SUBPROBLEM_HH
 #define DUMUX_STOKES2C_SUBPROBLEM_HH
 
-#include <dune/pdelab/gridoperator/gridoperator.hh>
-
 #include <dumux/freeflow/stokesnc/stokesncmodel.hh>
 #include <dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh>
 #include <dumux/multidomain/common/subdomainpropertydefaults.hh>
@@ -64,27 +62,6 @@ SET_TYPE_PROP(Stokes2cSubProblem,
               LocalResidual,
               StokesncCouplingLocalResidual<TypeTag>);
 
-SET_TYPE_PROP(Stokes2cSubProblem, Constraints,
-        Dune::PDELab::NonoverlappingConformingDirichletConstraints);
-
-SET_PROP(Stokes2cSubProblem, GridOperator)
-{
- private:
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, ConstraintsTrafo) ConstraintsTrafo;
-    typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
-    typedef Dumux::PDELab::MultiDomainLocalOperator<TypeTag> LocalOperator;
-
-    enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
- public:
-    typedef Dune::PDELab::GridOperator<GridFunctionSpace,
-        GridFunctionSpace, LocalOperator,
-        Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
-        Scalar, Scalar, Scalar,
-        ConstraintsTrafo, ConstraintsTrafo,
-        true> type;
-};
-
 SET_PROP(Stokes2cSubProblem, FluidSystem)
 {
 private: