diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9dfe1b77652d438dac3f34f17894fdc67c224184..8a36b953e73d84ee57a310777d4b77c480b3e99a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -17,6 +17,47 @@ Differences Between DuMuX 3.2 and DuMuX 3.1
 - __Van Genuchten__: Corrected VanGenuchten-Mualem implementation of `dkrn/dSw`
 - __AMGBackend__: The internal structure of the AMGBackend and the ParallelISTLHelper has been overhauled, as only used by the AMG, we did not make the changes backwards-compatible
 
+- __Change matrix block arrangement for staggered models__: The matrix block structure has been adapted such that it complies with the literature standard, i.e., having the velocity block (C) on `M[0][0]`
+rather than on `M[1][1]`. This also requires re-arranging the submodels and properties in dumux-multidomain such that the face-related classes and vector entries now appear before the cell-centered ones.
+
+```math
+M = \begin{pmatrix}
+  A & B^T\\
+  B & C
+\end{pmatrix}
+
+\qquad => \qquad
+
+M = \begin{pmatrix}
+    C & B^T\\
+    B & A
+    \end{pmatrix}
+```
+
+Backwards-compatibility can only be provided to a certain extent. The following changes need to made in the main file:
+
+1.) change the order of the arguments for the `assembler` such that it reads:
+```c++
+auto assembler = std::make_shared<Assembler>(std::make_tuple(ffProblem, ffProblem, otherProblem, ...),
+                                             std::make_tuple(ffGridGeometry->faceFVGridGeometryPtr(),
+                                                             ffFvGridGeometry->cellCenterFVGridGeometryPtr(),
+                                                             otherFvGridGeometry, ...),
+                                             std::make_tuple(ffGridVariables->faceGridVariablesPtr(),
+                                                             ffGridVariables->cellCenterGridVariablesPtr(),
+                                                             otherGridVariables, ...),
+                                             couplingManager,
+                                             timeLoop, solOld);
+
+// Not changing the arguments will yield a deprecation warning stating this hint but the code still compiles and runs.
+```
+
+2.) change the order of arguments in the `partial` function:
+```c++
+ffSol = partial(sol, ffFaceIdx, ffCellCenterIdx);
+
+// Not changing the argument will rise a compiler error which makes the MR not fully backwards-compatible.
+```
+
 ### Deprecated properties, to be removed after 3.2:
 
 ### Deprecated classes/files, to be removed after 3.2:
diff --git a/dumux/assembly/staggeredfvassembler.hh b/dumux/assembly/staggeredfvassembler.hh
index 786b99a533fed2679c420768e77e9093d2c19147..21f7f8d52d713771424fb71cdc908d3a1ecbacb5 100644
--- a/dumux/assembly/staggeredfvassembler.hh
+++ b/dumux/assembly/staggeredfvassembler.hh
@@ -25,9 +25,10 @@
 #ifndef DUMUX_STAGGERED_FV_ASSEMBLER_HH
 #define DUMUX_STAGGERED_FV_ASSEMBLER_HH
 
-#include <type_traits>
+#include <tuple>
+#include <memory>
 
-#include <dune/istl/matrixindexset.hh>
+#include <dune/common/indices.hh>
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/timeloop.hh>
@@ -66,7 +67,6 @@ class StaggeredFVAssembler: public MultiDomainFVAssembler<StaggeredMultiDomainTr
 
 public:
     using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
-    using FVGridGeometry [[deprecated("Use GridGeometry instead. FVGridGeometry will be removed after 3.1!")]] = GridGeometry;
     using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
     using CouplingManager = typename ParentType::CouplingManager;
 
@@ -75,8 +75,8 @@ public:
                          std::shared_ptr<const GridGeometry> gridGeometry,
                          std::shared_ptr<GridVariables> gridVariables)
     : ParentType(std::make_tuple(problem, problem),
-                 std::make_tuple(gridGeometry->cellCenterFVGridGeometryPtr(), gridGeometry->faceFVGridGeometryPtr()),
-                 std::make_tuple(gridVariables->cellCenterGridVariablesPtr(), gridVariables->faceGridVariablesPtr()),
+                 std::make_tuple(gridGeometry->faceFVGridGeometryPtr(), gridGeometry->cellCenterFVGridGeometryPtr()),
+                 std::make_tuple(gridVariables->faceGridVariablesPtr(), gridVariables->cellCenterGridVariablesPtr()),
                  std::make_shared<CouplingManager>())
     {
         static_assert(isImplicit, "Explicit assembler for stationary problem doesn't make sense!");
@@ -90,8 +90,8 @@ public:
                          std::shared_ptr<GridVariables> gridVariables,
                          std::shared_ptr<const TimeLoop> timeLoop)
     : ParentType(std::make_tuple(problem, problem),
-                 std::make_tuple(gridGeometry->cellCenterFVGridGeometryPtr(), gridGeometry->faceFVGridGeometryPtr()),
-                 std::make_tuple(gridVariables->cellCenterGridVariablesPtr(), gridVariables->faceGridVariablesPtr()),
+                 std::make_tuple(gridGeometry->faceFVGridGeometryPtr(), gridGeometry->cellCenterFVGridGeometryPtr()),
+                 std::make_tuple(gridVariables->faceGridVariablesPtr(), gridVariables->cellCenterGridVariablesPtr()),
                  std::make_shared<CouplingManager>(),
                  timeLoop)
     {
@@ -106,8 +106,8 @@ public:
                          std::shared_ptr<const TimeLoop> timeLoop,
                          const SolutionVector& prevSol)
     : ParentType(std::make_tuple(problem, problem),
-                 std::make_tuple(gridGeometry->cellCenterFVGridGeometryPtr(), gridGeometry->faceFVGridGeometryPtr()),
-                 std::make_tuple(gridVariables->cellCenterGridVariablesPtr(), gridVariables->faceGridVariablesPtr()),
+                 std::make_tuple(gridGeometry->faceFVGridGeometryPtr(), gridGeometry->cellCenterFVGridGeometryPtr()),
+                 std::make_tuple(gridVariables->faceGridVariablesPtr(), gridVariables->cellCenterGridVariablesPtr()),
                  std::make_shared<CouplingManager>(),
                  timeLoop,
                  prevSol)
@@ -116,18 +116,12 @@ public:
         this->couplingManager_->setSubProblems(std::make_tuple(problem, problem));
     }
 
-
     auto& gridVariables()
     { return ParentType::gridVariables(Dune::index_constant<0>()); }
 
     const auto& gridVariables() const
     { return ParentType::gridVariables(Dune::index_constant<0>()); }
 
-    //! The global finite volume geometry
-    [[deprecated("Use gridGeometry() instead. fvGridGeometry() will be removed after 3.1!")]]
-    const GridGeometry& fvGridGeometry() const
-    { return gridGeometry(); }
-
     const GridGeometry& gridGeometry() const
     { return ParentType::gridGeometry(Dune::index_constant<0>()).actualGridGeometry(); }
 
diff --git a/dumux/common/typetraits/utility.hh b/dumux/common/typetraits/utility.hh
index 5c7097214310adfc7fda5d7eb5a47cf711c69153..20d69bb10291a164d0dd1c7817da71c299d5ac34 100644
--- a/dumux/common/typetraits/utility.hh
+++ b/dumux/common/typetraits/utility.hh
@@ -70,6 +70,34 @@ template <std::size_t n, std::size_t e>
 using makeIncompleteIntegerSequence =
         typename Detail::ConcatSeq<decltype(std::make_index_sequence<e>{}), e + 1, decltype(std::make_index_sequence<(n > e) ? (n - e - 1) : 0>{})>::type;
 
+/*
+ * \ingroup TypeTraits
+ * \brief add an offset to an index sequence
+ * \tparam offset the offset
+ * \tparam is the index sequence
+ *
+ * see https://stackoverflow.com/a/35625414
+ */
+template <std::size_t offset, std::size_t ... is>
+constexpr std::index_sequence<(offset + is)...> addOffsetToIndexSequence(std::index_sequence<is...>)
+{ return {}; }
+
+/*
+ * \ingroup TypeTraits
+ * \brief create an index sequence starting from an offset
+ * \tparam offset the offset
+ * \tparam n the length of the sequence
+ *
+ * example: makeIndexSequenceWithOffset<2, 3> = [2,3,4]
+ *
+ * see https://stackoverflow.com/a/35625414
+ */
+template <std::size_t offset, std::size_t n>
+constexpr auto makeIndexSequenceWithOffset()
+{
+    return addOffsetToIndexSequence<offset>(std::make_index_sequence<n>{});
+}
+
 } // end namespace Dumux
 
 #endif
diff --git a/dumux/discretization/staggered.hh b/dumux/discretization/staggered.hh
index 585bab9d8f24ec5e23273202ea18b91427bc5160..ca73668c30765f4ce4cce780e5f8d846716db1ba 100644
--- a/dumux/discretization/staggered.hh
+++ b/dumux/discretization/staggered.hh
@@ -159,10 +159,10 @@ template<class TypeTag>
 struct SolutionVector<TypeTag, TTag::StaggeredModel>
 {
 private:
-    using CellCenterSolutionVector = GetPropType<TypeTag, Properties::CellCenterSolutionVector>;
     using FaceSolutionVector = GetPropType<TypeTag, Properties::FaceSolutionVector>;
+    using CellCenterSolutionVector = GetPropType<TypeTag, Properties::CellCenterSolutionVector>;
 public:
-    using type = Dune::MultiTypeBlockVector<CellCenterSolutionVector, FaceSolutionVector>;
+    using type = Dune::MultiTypeBlockVector<FaceSolutionVector, CellCenterSolutionVector>;
 };
 
 //! Set the type of a global jacobian matrix from the solution types TODO: move to LinearAlgebra traits
@@ -191,11 +191,11 @@ public:
     using MatrixBlockFaceToCC = typename Dune::BCRSMatrix<MatrixLittleBlockFaceToCC>;
 
     // the row types
-    using RowCellCenter = typename Dune::MultiTypeBlockVector<MatrixBlockCCToCC, MatrixBlockCCToFace>;
-    using RowFace = typename Dune::MultiTypeBlockVector<MatrixBlockFaceToCC, MatrixBlockFaceToFace>;
+    using RowFace = typename Dune::MultiTypeBlockVector<MatrixBlockFaceToFace, MatrixBlockFaceToCC>;
+    using RowCellCenter = typename Dune::MultiTypeBlockVector<MatrixBlockCCToFace, MatrixBlockCCToCC>;
 
     // the jacobian matrix
-    using type = typename Dune::MultiTypeBlockMatrix<RowCellCenter, RowFace>;
+    using type = typename Dune::MultiTypeBlockMatrix<RowFace, RowCellCenter>;
 };
 
 } // namespace Properties
diff --git a/dumux/discretization/staggered/freeflow/facevariables.hh b/dumux/discretization/staggered/freeflow/facevariables.hh
index 89cdfa2e6fc00462b263d0d9cdba5d417c13a449..6eebd6e3c2441838a40ed4aa0b718bed0784d50c 100644
--- a/dumux/discretization/staggered/freeflow/facevariables.hh
+++ b/dumux/discretization/staggered/freeflow/facevariables.hh
@@ -26,6 +26,7 @@
 
 #include <array>
 #include <vector>
+#include <type_traits>
 
 namespace Dumux {
 
@@ -88,14 +89,17 @@ public:
     * \param fvGeometry The finite-volume geometry
     * \param scvf The sub-control volume face of interest
     */
-    template<class SolVector, class Problem, class Element,
+    template<class FaceSolution, class Problem, class Element,
              class FVElementGeometry, class SubControlVolumeFace>
-    void update(const SolVector& faceSol,
+    void update(const FaceSolution& faceSol,
                 const Problem& problem,
                 const Element& element,
                 const FVElementGeometry& fvGeometry,
                 const SubControlVolumeFace& scvf)
     {
+        static_assert(std::decay_t<decltype(faceSol[0])>::dimension == 1,
+                      "\n\n\nVelocity primary variable must be a scalar value. \n\n Make sure to use\n\n ffSol = partial(sol, ffFaceIdx, ffCellCenterIdx);\n\n");
+
         inAxisVelocities_.self = faceSol[scvf.dofIndex()];
         inAxisVelocities_.opposite = faceSol[scvf.dofIndexOpposingFace()];
 
diff --git a/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh b/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
index f321e56e2cdf2d90062b24ba68d5d80461525c2d..051ca1c4133e54b4bc87429cd5e74020a176de27 100644
--- a/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
+++ b/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
@@ -51,8 +51,8 @@ struct StaggeredFreeFlowDefaultFVGridGeometryTraits
 
     struct DofTypeIndices
     {
-        using CellCenterIdx = Dune::index_constant<0>;
-        using FaceIdx = Dune::index_constant<1>;
+        using FaceIdx = Dune::index_constant<0>;
+        using CellCenterIdx = Dune::index_constant<1>;
     };
 
     template<class GridGeometry>
diff --git a/dumux/discretization/staggered/freeflow/gridvolumevariables.hh b/dumux/discretization/staggered/freeflow/gridvolumevariables.hh
index 7170c40cd65878be68dcc7583b57c5a8a0f6c18b..badedcf28729620dffd1ae17fe2df2e7ff883aa7 100644
--- a/dumux/discretization/staggered/freeflow/gridvolumevariables.hh
+++ b/dumux/discretization/staggered/freeflow/gridvolumevariables.hh
@@ -162,8 +162,13 @@ public:
     template<class GridGeometry, class SolutionVector>
     void update(const GridGeometry& gridGeometry, const SolutionVector& sol)
     {
-        auto numScv = gridGeometry.numScv();
-        volumeVariables_.resize(numScv);
+        if (sol.size() != gridGeometry.numScv())
+            DUNE_THROW(Dune::InvalidStateException, "The solution vector passed to the GridVolumeVariables has the wrong size.\n"
+                                                     << "Make sure to initialize the gridVariables correctly: \n\n"
+                                                     << "auto ffSol = partial(sol, ffFaceIdx, ffCellCenterIdx); \n"
+                                                     << "ffGridVariables->init(ffSol);\n\n");
+
+        volumeVariables_.resize(gridGeometry.numScv());
 
         for (const auto& element : elements(gridGeometry.gridView()))
         {
@@ -240,7 +245,14 @@ public:
     StaggeredGridVolumeVariables(const Problem& problem) : problemPtr_(&problem) {}
 
     template<class GridGeometry, class SolutionVector>
-    void update(const GridGeometry& gridGeometry, const SolutionVector& sol) {}
+    void update(const GridGeometry& gridGeometry, const SolutionVector& sol)
+    {
+        if (sol.size() != gridGeometry.numScv())
+            DUNE_THROW(Dune::InvalidStateException, "The solution vector passed to the GridVolumeVariables has the wrong size.\n"
+                                                     << "Make sure to initialize the gridVariables correctly: \n\n"
+                                                     << "auto ffSol = partial(sol, ffFaceIdx, ffCellCenterIdx); \n"
+                                                     << "ffGridVariables->init(ffSol);\n\n");
+    }
 
     const Problem& problem() const
     { return *problemPtr_;}
diff --git a/dumux/discretization/staggered/gridfacevariables.hh b/dumux/discretization/staggered/gridfacevariables.hh
index f6c5eb63155dd862f2bdd9ff8e8390011917f4b8..ccc8859b48176b79a1e84de90a0547fec5257763 100644
--- a/dumux/discretization/staggered/gridfacevariables.hh
+++ b/dumux/discretization/staggered/gridfacevariables.hh
@@ -27,6 +27,7 @@
 //! make the local view function available whenever we use this class
 #include <dumux/discretization/localview.hh>
 #include <dumux/discretization/staggered/elementfacevariables.hh>
+#include <dumux/discretization/staggered/facesolution.hh>
 
 namespace Dumux {
 
@@ -82,18 +83,18 @@ public:
 
     //! Update all face variables
     template<class GridGeometry, class SolutionVector>
-    void update(const GridGeometry& gridGeometry, const SolutionVector& faceSol)
+    void update(const GridGeometry& gridGeometry, const SolutionVector& sol)
     {
-
         faceVariables_.resize(gridGeometry.numScvf());
 
-        for(auto&& element : elements(gridGeometry.gridView()))
+        for (auto&& element : elements(gridGeometry.gridView()))
         {
             auto fvGeometry = localView(gridGeometry);
             fvGeometry.bindElement(element);
 
-            for(auto&& scvf : scvfs(fvGeometry))
+            for (auto&& scvf : scvfs(fvGeometry))
             {
+                const auto faceSol = StaggeredFaceSolution(scvf, sol, gridGeometry);
                 faceVariables_[scvf.index()].update(faceSol, problem(), element, fvGeometry, scvf);
             }
         }
diff --git a/dumux/discretization/staggered/gridvariables.hh b/dumux/discretization/staggered/gridvariables.hh
index 4475e690cd3427b6cfa06ad099c349500aea8363..be35f83ac582c677e4e144c0f35849e74dc48bb3 100644
--- a/dumux/discretization/staggered/gridvariables.hh
+++ b/dumux/discretization/staggered/gridvariables.hh
@@ -47,7 +47,6 @@ public:
     //! export primary variable type
     using PrimaryVariables = typename VolumeVariables::PrimaryVariables;
     using GridGeometry = typename ActualGridVariables::GridGeometry;
-    using FVGridGeometry [[deprecated ("Use GridGeometry instead. Will be removed after 3.1!")]]= GridGeometry;
 
     explicit StaggeredGridVariablesView(ActualGridVariables* gridVariables)
     : gridVariables_(gridVariables) {}
@@ -92,10 +91,6 @@ public:
     GridFaceVariables& prevGridFaceVars()
     { return gridVariables_->prevGridFaceVars(); }
 
-    //! return the fv grid geometry
-    [[deprecated("Use gridGeometry() instead. fvGridGeometry() will be removed after 3.1!")]]
-    const GridGeometry& fvGridGeometry() const
-    { return (*gridVariables_->gridGeometry_);    }
     //! return the fv grid geometry
     const GridGeometry& gridGeometry() const
     { return (*gridVariables_->gridGeometry_);    }
@@ -133,14 +128,6 @@ public:
         this->prevGridVolVars().update(this->gridGeometry(), curSol);
     }
 
-    //! initialize all variables (instationary case)
-    template<class SolVector>
-    [[deprecated("Use init with one argument")]]
-    void init(const SolVector& curSol, const SolVector& initSol)
-    {
-        init(initSol);
-    }
-
     //! update the volume variables and the flux variables cache
     template<class SolVector>
     void update(const SolVector& curSol)
@@ -178,14 +165,6 @@ public:
         this->prevGridFaceVars().update(this->gridGeometry(), curSol);
     }
 
-    //! initialize all variables (instationary case)
-    template<class SolVector>
-    [[deprecated("Use init with one argument")]]
-    void init(const SolVector& curSol, const SolVector& initSol)
-    {
-        init(initSol);
-    }
-
     //! update the face variables
     template<class SolVector>
     void update(const SolVector& curSol)
@@ -260,14 +239,6 @@ public:
         prevGridFaceVariables_.update(*this->gridGeometry_, curSol[faceIdx]);
     }
 
-    //! initialize all variables (instationary case)
-    template<class SolutionVector>
-    [[deprecated("Use init with one argument")]]
-    void init(const SolutionVector& curSol, const SolutionVector& initSol)
-    {
-        init(initSol);
-    }
-
     //! Sets the current state as the previous for next time step
     //! this has to be called at the end of each time step
     void advanceTimeStep()
diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
index 52ffacc88fbc32b9c4d066979c2e4c908272e58d..ed6f316e384536445503fbc5132fdc3fffdf07c1 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
@@ -537,10 +537,10 @@ class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEne
 {
     using ParentType = StokesDarcyCouplingDataImplementationBase<MDTraits, CouplingManager>;
     using Scalar = typename MDTraits::Scalar;
-    static constexpr auto stokesIdx = typename MDTraits::template SubDomain<0>::Index();
-    static constexpr auto darcyIdx = typename MDTraits::template SubDomain<2>::Index();
-    static constexpr auto stokesCellCenterIdx = stokesIdx;
-    static constexpr auto stokesFaceIdx = typename MDTraits::template SubDomain<1>::Index();
+    static constexpr auto stokesIdx = CouplingManager::stokesIdx;
+    static constexpr auto darcyIdx = CouplingManager::darcyIdx;
+    static constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+    static constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
 
     // the sub domain type tags
     template<std::size_t id>
@@ -699,10 +699,10 @@ class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEne
 {
     using ParentType = StokesDarcyCouplingDataImplementationBase<MDTraits, CouplingManager>;
     using Scalar = typename MDTraits::Scalar;
-    static constexpr auto stokesIdx = typename MDTraits::template SubDomain<0>::Index();
-    static constexpr auto darcyIdx = typename MDTraits::template SubDomain<2>::Index();
-    static constexpr auto stokesCellCenterIdx = stokesIdx;
-    static constexpr auto stokesFaceIdx = typename MDTraits::template SubDomain<1>::Index();
+    static constexpr auto stokesIdx = CouplingManager::stokesIdx;
+    static constexpr auto darcyIdx = CouplingManager::darcyIdx;
+    static constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+    static constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
 
     // the sub domain type tags
     template<std::size_t id>
diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh b/dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh
index 47c20dd3b5882ded199337d32aa3a158ed850b46..57d3a862d851d14dc6e8309e4245b2fbd0914053 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh
@@ -51,11 +51,9 @@ class StokesDarcyCouplingManager
     using ParentType = StaggeredCouplingManager<MDTraits>;
 
 public:
-    static constexpr auto stokesCellCenterIdx = typename MDTraits::template SubDomain<0>::Index();
-    static constexpr auto stokesFaceIdx = typename MDTraits::template SubDomain<1>::Index();
-    static constexpr auto cellCenterIdx = typename MDTraits::template SubDomain<0>::Index();
-    static constexpr auto faceIdx = typename MDTraits::template SubDomain<1>::Index();
-    static constexpr auto stokesIdx = stokesCellCenterIdx;
+    static constexpr auto stokesFaceIdx = typename MDTraits::template SubDomain<0>::Index();
+    static constexpr auto stokesCellCenterIdx = typename MDTraits::template SubDomain<1>::Index();
+    static constexpr auto stokesIdx = stokesFaceIdx;
     static constexpr auto darcyIdx = typename MDTraits::template SubDomain<2>::Index();
 
 private:
@@ -94,7 +92,7 @@ private:
 
     using VelocityVector = typename Element<stokesIdx>::Geometry::GlobalCoordinate;
 
-    using CouplingMapper = StokesDarcyCouplingMapper<MDTraits>;
+    using CouplingMapper = StokesDarcyCouplingMapper;
 
     struct StationaryStokesCouplingContext
     {
@@ -322,7 +320,7 @@ public:
                                const LocalAssemblerI& localAssemblerI,
                                Dune::index_constant<stokesFaceIdx> domainJ,
                                const std::size_t dofIdxGlobalJ,
-                               const PrimaryVariables<1>& priVars,
+                               const PrimaryVariables<stokesFaceIdx>& priVars,
                                int pvIdxJ)
     {
         this->curSol()[domainJ][dofIdxGlobalJ] = priVars;
diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingmapper.hh b/dumux/multidomain/boundary/stokesdarcy/couplingmapper.hh
index 748060bf98390ec3ee84f2b5e451d0eda083b4eb..19568f399eb5ac45602a9ac0a4699ae780020799 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingmapper.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingmapper.hh
@@ -29,7 +29,6 @@
 #include <unordered_map>
 #include <vector>
 
-#include <dune/common/exceptions.hh>
 #include <dumux/discretization/method.hh>
 
 namespace Dumux {
@@ -38,20 +37,8 @@ namespace Dumux {
  * \ingroup StokesDarcyCoupling
  * \brief Coupling mapper for Stokes and Darcy domains with equal dimension.
  */
-template<class MDTraits>
 class StokesDarcyCouplingMapper
 {
-    using Scalar = typename MDTraits::Scalar;
-
-public:
-    static constexpr auto stokesCellCenterIdx = typename MDTraits::template SubDomain<0>::Index();
-    static constexpr auto stokesFaceIdx = typename MDTraits::template SubDomain<1>::Index();
-    static constexpr auto cellCenterIdx = typename MDTraits::template SubDomain<0>::Index();
-    static constexpr auto faceIdx = typename MDTraits::template SubDomain<1>::Index();
-    static constexpr auto stokesIdx = stokesCellCenterIdx;
-    static constexpr auto darcyIdx = typename MDTraits::template SubDomain<2>::Index();
-
-private:
     struct ElementMapInfo
     {
         std::size_t eIdx;
@@ -71,8 +58,8 @@ public:
                                         Stencils& stokesCellCenterToDarcyStencils,
                                         Stencils& stokesFaceToDarcyStencils)
     {
-        const auto& stokesFvGridGeometry = couplingManager.problem(stokesIdx).gridGeometry();
-        const auto& darcyFvGridGeometry = couplingManager.problem(darcyIdx).gridGeometry();
+        const auto& stokesFvGridGeometry = couplingManager.problem(CouplingManager::stokesIdx).gridGeometry();
+        const auto& darcyFvGridGeometry = couplingManager.problem(CouplingManager::darcyIdx).gridGeometry();
 
         static_assert(std::decay_t<decltype(stokesFvGridGeometry)>::discMethod == DiscretizationMethod::staggered,
                       "The free flow domain must use the staggered discretization");
@@ -126,7 +113,7 @@ public:
                 // find the corresponding Darcy sub control volume face
                 for (const auto& darcyScvf : scvfs(darcyFvGeometry))
                 {
-                    const Scalar distance = (darcyScvf.center() - scvf.center()).two_norm();
+                    const auto distance = (darcyScvf.center() - scvf.center()).two_norm();
 
                     if (distance < eps)
                     {
diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh
index 13cf56370a65105b1aa3a72b46f543b1ac29d6be..db59dce18ffe842fd8faeca6abca96059088b054 100644
--- a/dumux/multidomain/fvassembler.hh
+++ b/dumux/multidomain/fvassembler.hh
@@ -34,6 +34,7 @@
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/timeloop.hh>
+#include <dumux/common/typetraits/utility.hh>
 #include <dumux/discretization/method.hh>
 #include <dumux/assembly/diffmethod.hh>
 #include <dumux/assembly/jacobianpattern.hh>
@@ -43,6 +44,8 @@
 #include "subdomainboxlocalassembler.hh"
 #include "subdomainstaggeredlocalassembler.hh"
 
+#include <dumux/discretization/method.hh>
+
 namespace Dumux {
 
 /*!
@@ -75,9 +78,6 @@ public:
     template<std::size_t id>
     using GridGeometry = typename MDTraits::template SubDomain<id>::GridGeometry;
 
-    template<std::size_t id>
-    using FVGridGeometry [[deprecated("Use GridGeometry instead. FVGridGeometry will be removed after 3.1!")]] = GridGeometry<id>;
-
     template<std::size_t id>
     using Problem = typename MDTraits::template SubDomain<id>::Problem;
 
@@ -197,6 +197,48 @@ public:
         std::cout << "Instantiated assembler for an instationary problem." << std::endl;
     }
 
+    template<class DeprecatedGridGeometryTuple, class DeprecatedGridVariablesTuple,
+             bool isStaggered = std::tuple_element_t<0, GridGeometryTuple>::element_type::discMethod == DiscretizationMethod::staggered,
+             bool isDeprecated = !std::is_convertible<DeprecatedGridVariablesTuple, GridVariablesTuple>::value,
+             std::enable_if_t<isStaggered && isDeprecated, int> = 0>
+    [[deprecated("Please change the order within your tuples for gridGeometry and gridVariables. Will be removed after release 3.2!")]]
+    MultiDomainFVAssembler(ProblemTuple&& problem,
+                           DeprecatedGridGeometryTuple&& gridGeometry,
+                           DeprecatedGridVariablesTuple&& gridVariables,
+                           std::shared_ptr<CouplingManager> couplingManager)
+   : MultiDomainFVAssembler(std::forward<ProblemTuple>(problem), gridGeometry, gridVariables, couplingManager,
+                            makeIndexSequenceWithOffset<2,std::tuple_size<DeprecatedGridGeometryTuple>::value - 2>())
+   {}
+
+    template<class DeprecatedGridGeometryTuple, class DeprecatedGridVariablesTuple,
+             bool isStaggered = std::tuple_element_t<0, GridGeometryTuple>::element_type::discMethod == DiscretizationMethod::staggered,
+             bool isDeprecated = !std::is_convertible<DeprecatedGridVariablesTuple, GridVariablesTuple>::value,
+             std::enable_if_t<isStaggered && isDeprecated, int> = 0>
+    [[deprecated("Please change the order within your tuples for gridGeometry and gridVariables. Will be removed after release 3.2!")]]
+    MultiDomainFVAssembler(ProblemTuple&& problem,
+                           DeprecatedGridGeometryTuple&& gridGeometry,
+                           DeprecatedGridVariablesTuple&& gridVariables,
+                           std::shared_ptr<CouplingManager> couplingManager,
+                           std::shared_ptr<const TimeLoop> timeLoop,
+                           const SolutionVector& prevSol)
+   : MultiDomainFVAssembler(std::forward<ProblemTuple>(problem), gridGeometry, gridVariables, couplingManager, timeLoop, prevSol,
+                            makeIndexSequenceWithOffset<2,std::tuple_size<DeprecatedGridGeometryTuple>::value - 2>())
+   {}
+
+    template<class DeprecatedGridGeometryTuple, class DeprecatedGridVariablesTuple,
+             bool isStaggered = std::tuple_element_t<0, GridGeometryTuple>::element_type::discMethod == DiscretizationMethod::staggered,
+             bool isDeprecated = !std::is_convertible<DeprecatedGridVariablesTuple, GridVariablesTuple>::value,
+             std::enable_if_t<isStaggered && isDeprecated, int> = 0>
+    [[deprecated("Please change the order within your tuples for gridGeometry and gridVariables. Will be removed after release 3.2!")]]
+    MultiDomainFVAssembler(ProblemTuple&& problem,
+                           DeprecatedGridGeometryTuple&& gridGeometry,
+                           DeprecatedGridVariablesTuple&& gridVariables,
+                           std::shared_ptr<CouplingManager> couplingManager,
+                           std::shared_ptr<const TimeLoop> timeLoop)
+   : MultiDomainFVAssembler(std::forward<ProblemTuple>(problem), gridGeometry, gridVariables, couplingManager, timeLoop, nullptr,
+                            makeIndexSequenceWithOffset<2,std::tuple_size<DeprecatedGridGeometryTuple>::value - 2>())
+   {}
+
     /*!
      * \brief Assembles the global Jacobian of the residual
      *        and the residual for the current solution.
@@ -358,12 +400,6 @@ public:
     const auto& problem(Dune::index_constant<i> domainId) const
     { return *std::get<domainId>(problemTuple_); }
 
-    //! the finite volume grid geometry of domain i
-    template<std::size_t i>
-    [[deprecated("Use gridGeometry() instead. fvGridGeometry() will be removed after 3.1!")]]
-    const auto& fvGridGeometry(Dune::index_constant<i> domainId) const
-    { return gridGeometry(domainId); }
-
     //! the finite volume grid geometry of domain i
     template<std::size_t i>
     const auto& gridGeometry(Dune::index_constant<i> domainId) const
@@ -525,6 +561,71 @@ private:
                                                         domainJ, gridGeometry(domainJ));
     }
 
+    template<class DeprecatedGridGeometryTuple, class DeprecatedGridVariablesTuple, std::size_t... I>
+    [[deprecated("\n\n Please change the order within your tuples for gridGeometry and gridVariables. Usage:\n\n"
+                 "auto assembler = std::make_shared<Assembler>(std::make_tuple(ffProblem, ffProblem, otherProblem, ...),\n"
+                 "                                             std::make_tuple(ffGridGeometry->faceFVGridGeometryPtr(),\n"
+                 "                                                             ffFvGridGeometry->cellCenterFVGridGeometryPtr(),\n"
+                 "                                                             otherFvGridGeometry, ...),\n"
+                 "                                             std::make_tuple(ffGridVariables->faceGridVariablesPtr(),\n"
+                 "                                                             ffGridVariables->cellCenterGridVariablesPtr(),\n"
+                 "                                                             otherGridVariables, ...),\n"
+                 "                                             couplingManager);\n\n"
+                 "Will be removed after release 3.2!")]]
+    MultiDomainFVAssembler(ProblemTuple&& problem,
+                           DeprecatedGridGeometryTuple&& gridGeometry,
+                           DeprecatedGridVariablesTuple&& gridVariables,
+                           std::shared_ptr<CouplingManager> couplingManager,
+                           std::index_sequence<I...>)
+   : couplingManager_(couplingManager)
+   , problemTuple_(problem)
+   , gridGeometryTuple_(std::make_tuple(std::move(std::get<1>(gridGeometry)),
+                                        std::move(std::get<0>(gridGeometry)),
+                                        std::move(std::get<I>(gridGeometry))...))
+   , gridVariablesTuple_(std::make_tuple(std::move(std::get<1>(gridVariables)),
+                                         std::move(std::get<0>(gridVariables)),
+                                         std::move(std::get<I>(gridVariables))...))
+   , timeLoop_()
+   , isStationaryProblem_(true)
+   {
+       static_assert(isImplicit(), "Explicit assembler for stationary problem doesn't make sense!");
+       std::cout << "Instantiated assembler for a stationary problem." << std::endl;
+   }
+
+    template<class DeprecatedGridGeometryTuple, class DeprecatedGridVariablesTuple, std::size_t... I>
+    [[deprecated("\n\n Please change the order within your tuples for gridGeometry and gridVariables. Usage:\n\n"
+                 "auto assembler = std::make_shared<Assembler>(std::make_tuple(ffProblem, ffProblem, otherProblem, ...),\n"
+                 "                                             std::make_tuple(ffGridGeometry->faceFVGridGeometryPtr(),\n"
+                 "                                                             ffFvGridGeometry->cellCenterFVGridGeometryPtr(),\n"
+                 "                                                             otherFvGridGeometry, ...),\n"
+                 "                                             std::make_tuple(ffGridVariables->faceGridVariablesPtr(),\n"
+                 "                                                             ffGridVariables->cellCenterGridVariablesPtr(),\n"
+                 "                                                             otherGridVariables, ...),\n"
+                 "                                             couplingManager,\n"
+                 "                                             timeLoop, solOld);\n\n"
+                 "Will be removed after release 3.2!")]]
+    MultiDomainFVAssembler(ProblemTuple&& problem,
+                           DeprecatedGridGeometryTuple&& gridGeometry,
+                           DeprecatedGridVariablesTuple&& gridVariables,
+                           std::shared_ptr<CouplingManager> couplingManager,
+                           std::shared_ptr<const TimeLoop> timeLoop,
+                           const SolutionVector& prevSol,
+                           std::index_sequence<I...>)
+   : couplingManager_(couplingManager)
+   , problemTuple_(problem)
+   , gridGeometryTuple_(std::make_tuple(std::move(std::get<1>(gridGeometry)),
+                                        std::move(std::get<0>(gridGeometry)),
+                                        std::move(std::get<I>(gridGeometry))...))
+   , gridVariablesTuple_(std::make_tuple(std::move(std::get<1>(gridVariables)),
+                                         std::move(std::get<0>(gridVariables)),
+                                         std::move(std::get<I>(gridVariables))...))
+   , timeLoop_(timeLoop)
+   , prevSol_(&prevSol)
+   , isStationaryProblem_(false)
+   {
+       std::cout << "Instantiated assembler for an instationary problem." << std::endl;
+   }
+
     //! pointer to the problem to be solved
     ProblemTuple problemTuple_;
 
diff --git a/dumux/multidomain/staggeredtraits.hh b/dumux/multidomain/staggeredtraits.hh
index 12172fefa3027ac3c0f9ba15f510c7b317c746dd..0d9edcc8554214f6f603debea8c118833a90fb2e 100644
--- a/dumux/multidomain/staggeredtraits.hh
+++ b/dumux/multidomain/staggeredtraits.hh
@@ -53,11 +53,11 @@ struct SubDomainFVGridGeometryImpl
 
 template<template<std::size_t> class SubDomainTypeTag>
 struct SubDomainFVGridGeometryImpl<SubDomainTypeTag, 0>
-{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::GridGeometry>::CellCenterFVGridGeometryType; };
+{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::GridGeometry>::FaceFVGridGeometryType; };
 
 template<template<std::size_t> class SubDomainTypeTag>
 struct SubDomainFVGridGeometryImpl<SubDomainTypeTag, 1>
-{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::GridGeometry>::FaceFVGridGeometryType; };
+{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::GridGeometry>::CellCenterFVGridGeometryType; };
 //////////////////////////////////////////////////////////
 
 //////////////////////////////////////////////////////////
@@ -67,11 +67,11 @@ struct SubDomainGridVariablesImpl
 
 template<template<std::size_t> class SubDomainTypeTag>
 struct SubDomainGridVariablesImpl<SubDomainTypeTag, 0>
-{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::GridVariables>::CellCenterGridVariablesType; };
+{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::GridVariables>::FaceGridVariablesType; };
 
 template<template<std::size_t> class SubDomainTypeTag>
 struct SubDomainGridVariablesImpl<SubDomainTypeTag, 1>
-{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::GridVariables>::FaceGridVariablesType; };
+{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::GridVariables>::CellCenterGridVariablesType; };
 //////////////////////////////////////////////////////////
 
 //////////////////////////////////////////////////////////
@@ -81,11 +81,11 @@ struct SubDomainPrimaryVariablesImpl
 
 template<template<std::size_t> class SubDomainTypeTag>
 struct SubDomainPrimaryVariablesImpl<SubDomainTypeTag, 0>
-{ using type = GetPropType<SubDomainTypeTag<0>, Properties::CellCenterPrimaryVariables>; };
+{ using type = GetPropType<SubDomainTypeTag<0>, Properties::FacePrimaryVariables>; };
 
 template<template<std::size_t> class SubDomainTypeTag>
 struct SubDomainPrimaryVariablesImpl<SubDomainTypeTag, 1>
-{ using type = GetPropType<SubDomainTypeTag<0>, Properties::FacePrimaryVariables>; };
+{ using type = GetPropType<SubDomainTypeTag<0>, Properties::CellCenterPrimaryVariables>; };
 //////////////////////////////////////////////////////////
 
 template<class Scalar, int numEq>
@@ -105,12 +105,12 @@ struct SubDomainJacobianMatrixImpl
 template<template<std::size_t> class SubDomainTypeTag>
 struct SubDomainJacobianMatrixImpl<SubDomainTypeTag, 0>
 { using type = typename JacobianTypeImpl<GetPropType<SubDomainTypeTag<0>, Properties::Scalar>,
-                                         getPropValue<SubDomainTypeTag<0>, Properties::NumEqCellCenter>()>::type; };
+                                         getPropValue<SubDomainTypeTag<0>, Properties::NumEqFace>()>::type; };
 
 template<template<std::size_t> class SubDomainTypeTag>
 struct SubDomainJacobianMatrixImpl<SubDomainTypeTag, 1>
 { using type = typename JacobianTypeImpl<GetPropType<SubDomainTypeTag<1>, Properties::Scalar>,
-                                         getPropValue<SubDomainTypeTag<0>, Properties::NumEqFace>()>::type; };
+                                         getPropValue<SubDomainTypeTag<0>, Properties::NumEqCellCenter>()>::type; };
 //////////////////////////////////////////////////////////
 
 //////////////////////////////////////////////////////////
@@ -120,11 +120,11 @@ struct SubDomainSolutionVectorImpl
 
 template<template<std::size_t> class SubDomainTypeTag>
 struct SubDomainSolutionVectorImpl<SubDomainTypeTag, 0>
-{ using type = GetPropType<SubDomainTypeTag<0>, Properties::CellCenterSolutionVector>; };
+{ using type = GetPropType<SubDomainTypeTag<0>, Properties::FaceSolutionVector>; };
 
 template<template<std::size_t> class SubDomainTypeTag>
 struct SubDomainSolutionVectorImpl<SubDomainTypeTag, 1>
-{ using type = GetPropType<SubDomainTypeTag<0>, Properties::FaceSolutionVector>; };
+{ using type = GetPropType<SubDomainTypeTag<0>, Properties::CellCenterSolutionVector>; };
 //////////////////////////////////////////////////////////
 
 } // end namespace Staggered
diff --git a/dumux/multidomain/subdomainstaggeredlocalassembler.hh b/dumux/multidomain/subdomainstaggeredlocalassembler.hh
index c14390b471022b42e6bec91a1159d02151c4cf1a..c02a12ed781be6d22721b14199ce99aa91261db3 100644
--- a/dumux/multidomain/subdomainstaggeredlocalassembler.hh
+++ b/dumux/multidomain/subdomainstaggeredlocalassembler.hh
@@ -83,8 +83,8 @@ class SubDomainStaggeredLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag
 
 public:
     static constexpr auto domainId = typename Dune::index_constant<id>();
-    static constexpr auto cellCenterId = typename Dune::index_constant<0>();
-    static constexpr auto faceId = typename Dune::index_constant<1>();
+    static constexpr auto cellCenterId = GridGeometry::cellCenterIdx();
+    static constexpr auto faceId = GridGeometry::faceIdx();
 
     static constexpr auto numEqCellCenter = CellCenterResidualValue::dimension;
     static constexpr auto faceOffset = numEqCellCenter;
@@ -131,7 +131,16 @@ public:
         this->asImp_().bindLocalViews();
         this->elemBcTypes().update(problem(), this->element(), this->fvGeometry());
 
-        assembleResidualImpl_(domainId, res);
+        if constexpr (domainId == cellCenterId)
+        {
+            const auto cellCenterGlobalI = problem().gridGeometry().elementMapper().index(this->element());
+            res[cellCenterGlobalI] = this->asImp_().assembleCellCenterResidualImpl();
+        }
+        else
+        {
+            for (auto&& scvf : scvfs(this->fvGeometry()))
+                res[scvf.dofIndex()] +=  this->asImp_().assembleFaceResidualImpl(scvf);
+        }
     }
 
     /*!
@@ -315,25 +324,9 @@ public:
 
 private:
 
-    //! Assembles the residuals for the cell center dofs.
-    template<class SubSol>
-    void assembleResidualImpl_(Dune::index_constant<0>, SubSol& res)
-    {
-        const auto cellCenterGlobalI = problem().gridGeometry().elementMapper().index(this->element());
-        res[cellCenterGlobalI] = this->asImp_().assembleCellCenterResidualImpl();
-    }
-
-    //! Assembles the residuals for the face dofs.
-    template<class SubSol>
-    void assembleResidualImpl_(Dune::index_constant<1>, SubSol& res)
-    {
-        for (auto&& scvf : scvfs(this->fvGeometry()))
-            res[scvf.dofIndex()] +=  this->asImp_().assembleFaceResidualImpl(scvf);
-    }
-
     //! Assembles the residuals and derivatives for the cell center dofs.
     template<class JacobianMatrixRow, class SubSol, class GridVariablesTuple>
-    auto assembleJacobianAndResidualImpl_(Dune::index_constant<0>, JacobianMatrixRow& jacRow, SubSol& res, GridVariablesTuple& gridVariables)
+    auto assembleJacobianAndResidualImpl_(Dune::index_constant<cellCenterId>, JacobianMatrixRow& jacRow, SubSol& res, GridVariablesTuple& gridVariables)
     {
         auto& gridVariablesI = *std::get<domainId>(gridVariables);
         const auto cellCenterGlobalI = problem().gridGeometry().elementMapper().index(this->element());
@@ -355,7 +348,7 @@ private:
 
     //! Assembles the residuals and derivatives for the face dofs.
     template<class JacobianMatrixRow, class SubSol, class GridVariablesTuple>
-    void assembleJacobianAndResidualImpl_(Dune::index_constant<1>, JacobianMatrixRow& jacRow, SubSol& res, GridVariablesTuple& gridVariables)
+    void assembleJacobianAndResidualImpl_(Dune::index_constant<faceId>, JacobianMatrixRow& jacRow, SubSol& res, GridVariablesTuple& gridVariables)
     {
         auto& gridVariablesI = *std::get<domainId>(gridVariables);
         const auto residual = this->asImp_().assembleFaceJacobianAndResidualImpl(jacRow[domainId], gridVariablesI);
@@ -508,8 +501,8 @@ class SubDomainStaggeredLocalAssembler<id, TypeTag, Assembler, DiffMethod::numer
     static constexpr bool enableGridFluxVarsCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
     static constexpr int maxNeighbors = 4*(2*ModelTraits::dim());
     static constexpr auto domainI = Dune::index_constant<id>();
-    static constexpr auto cellCenterId = typename Dune::index_constant<0>();
-    static constexpr auto faceId = typename Dune::index_constant<1>();
+    static constexpr auto cellCenterId = GridGeometry::cellCenterIdx();
+    static constexpr auto faceId = GridGeometry::faceIdx();
 
     static constexpr auto numEq = ModelTraits::numEq();
     static constexpr auto numEqCellCenter = CellCenterPrimaryVariables::dimension;
@@ -665,7 +658,7 @@ public:
             auto evaluateFaceDerivatives = [&](const std::size_t globalJ)
             {
                 // get the faceVars of the face with respect to which we are going to build the derivative
-                auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), this->curElemFaceVars(), scvf);
+                auto& faceVars = getFaceVarAccess_(gridVariables.curGridFaceVars(), this->curElemFaceVars(), scvf);
                 const auto origFaceVars = faceVars;
 
                 for (int pvIdx = 0; pvIdx < numEqFace; ++pvIdx)
@@ -725,7 +718,7 @@ public:
      * \return The element residual at the current solution.
      */
     template<class JacobianBlock, class GridVariables>
-    void assembleJacobianCellCenterCoupling(Dune::index_constant<1> domainJ, JacobianBlock& A,
+    void assembleJacobianCellCenterCoupling(Dune::index_constant<faceId> domainJ, JacobianBlock& A,
                                             const CellCenterResidualValue& origResidual, GridVariables& gridVariables)
     {
         /////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -745,7 +738,7 @@ public:
             const auto globalJ = scvfJ.dofIndex();
 
             // get the faceVars of the face with respect to which we are going to build the derivative
-            auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), this->curElemFaceVars(), scvfJ);
+            auto& faceVars = getFaceVarAccess_(gridVariables.curGridFaceVars(), this->curElemFaceVars(), scvfJ);
             const auto origFaceVars(faceVars);
 
             for (int pvIdx = 0; pvIdx < numEqFace; ++pvIdx)
@@ -847,7 +840,7 @@ public:
      * \return The element residual at the current solution.
      */
     template<class JacobianBlock, class ElementResidualVector, class GridVariables>
-    void assembleJacobianFaceCoupling(Dune::index_constant<0> domainJ, JacobianBlock& A,
+    void assembleJacobianFaceCoupling(Dune::index_constant<cellCenterId> domainJ, JacobianBlock& A,
                                       const ElementResidualVector& origResiduals, GridVariables& gridVariables)
     {
         /////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -1012,15 +1005,13 @@ public:
 
 private:
 
-    template<class T = TypeTag>
-    static typename std::enable_if<!getPropValue<T, Properties::EnableGridFaceVariablesCache>(), FaceVariables&>::type
-    getFaceVarAccess(GridFaceVariables& gridFaceVariables, ElementFaceVariables& elemFaceVars, const SubControlVolumeFace& scvf)
-    { return elemFaceVars[scvf]; }
-
-    template<class T = TypeTag>
-    static typename std::enable_if<getPropValue<T, Properties::EnableGridFaceVariablesCache>(), FaceVariables&>::type
-    getFaceVarAccess(GridFaceVariables& gridFaceVariables, ElementFaceVariables& elemFaceVars, const SubControlVolumeFace& scvf)
-    { return gridFaceVariables.faceVars(scvf.index()); }
+    FaceVariables& getFaceVarAccess_(GridFaceVariables& gridFaceVariables, ElementFaceVariables& elemFaceVars, const SubControlVolumeFace& scvf)
+    {
+        if constexpr (getPropValue<TypeTag, Properties::EnableGridFaceVariablesCache>())
+            return gridFaceVariables.faceVars(scvf.index());
+        else
+            return elemFaceVars[scvf];
+    }
 };
 
 } // end namespace Dumux
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/diffusionlawcomparison/main.cc b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/diffusionlawcomparison/main.cc
index 710f2b4dad028cab82ec91699c64964513720631..2580af165d51c8d745d78c1c4dd9ef2b7ffae91d 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/diffusionlawcomparison/main.cc
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/diffusionlawcomparison/main.cc
@@ -166,7 +166,7 @@ int main(int argc, char** argv) try
     sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
 
     // get a solution vector storing references to the two Stokes solution vectors
-    auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
 
     // apply initial solution for instationary problems
     stokesProblem->applyInitialSolution(stokesSol);
@@ -199,11 +199,11 @@ int main(int argc, char** argv) try
     // the assembler with time loop for instationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
                                                                  darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 stokesGridVariables->faceGridVariablesPtr(),
+                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
+                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
                                                                  darcyGridVariables),
                                                  couplingManager,
                                                  timeLoop,
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/main.cc b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/main.cc
index 4dba90b79b72e5439e2f4bb90e5cabb041a203da..56c348c5a53f330167452c25a51496c20f4d293d 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/main.cc
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/main.cc
@@ -154,7 +154,7 @@ int main(int argc, char** argv) try
     sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
 
     // get a solution vector storing references to the two Stokes solution vectors
-    auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
 
     // apply initial solution for instationary problems
     stokesProblem->applyInitialSolution(stokesSol);
@@ -186,11 +186,11 @@ int main(int argc, char** argv) try
     // the assembler with time loop for instationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
                                                                  darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 stokesGridVariables->faceGridVariablesPtr(),
+                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
+                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
                                                                  darcyGridVariables),
                                                  couplingManager,
                                                  timeLoop,
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc
index 9415666af71da7540e520d07ac9d7d2862ec8091..6bba58042f538a998d9ec22bf578718b9122da9b 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc
@@ -154,7 +154,7 @@ int main(int argc, char** argv) try
     sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
 
     // get a solution vector storing references to the two Stokes solution vectors
-    auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
 
     // apply initial solution for instationary problems
     stokesProblem->applyInitialSolution(stokesSol);
@@ -184,11 +184,11 @@ int main(int argc, char** argv) try
     // the assembler with time loop for instationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
                                                                  darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 stokesGridVariables->faceGridVariablesPtr(),
+                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
+                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
                                                                  darcyGridVariables),
                                                  couplingManager,
                                                  timeLoop,
diff --git a/test/multidomain/boundary/stokesdarcy/1p3c_1p3c/main.cc b/test/multidomain/boundary/stokesdarcy/1p3c_1p3c/main.cc
index 081b610b30ba9c6aa319bc6bcbf5ebb9bddce4ee..8c45d59943ce3bb2ac7e6089cd066dd7093d5653 100644
--- a/test/multidomain/boundary/stokesdarcy/1p3c_1p3c/main.cc
+++ b/test/multidomain/boundary/stokesdarcy/1p3c_1p3c/main.cc
@@ -153,7 +153,7 @@ int main(int argc, char** argv) try
     sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
 
     // get a solution vector storing references to the two Stokes solution vectors
-    auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
 
     // apply initial solution for instationary problems
     stokesProblem->applyInitialSolution(stokesSol);
@@ -186,11 +186,11 @@ int main(int argc, char** argv) try
     // the assembler with time loop for instationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
                                                                  darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 stokesGridVariables->faceGridVariablesPtr(),
+                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
+                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
                                                                  darcyGridVariables),
                                                  couplingManager,
                                                  timeLoop,
diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/main.cc b/test/multidomain/boundary/stokesdarcy/1p_1p/main.cc
index 779af6ea4fddfda7c894c14badc2523c9175eb86..28519da129b5f1633e9b7ca57c35f5c9ab1d385b 100644
--- a/test/multidomain/boundary/stokesdarcy/1p_1p/main.cc
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/main.cc
@@ -136,7 +136,7 @@ int main(int argc, char** argv) try
     sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
 
     // get a solution vector storing references to the two Stokes solution vectors
-    auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
 
     couplingManager->init(stokesProblem, darcyProblem, sol);
 
@@ -162,11 +162,11 @@ int main(int argc, char** argv) try
     // the assembler for a stationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
                                                                  darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 stokesGridVariables->faceGridVariablesPtr(),
+                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
+                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
                                                                  darcyGridVariables),
                                                  couplingManager);
 
diff --git a/test/multidomain/boundary/stokesdarcy/1p_2p/main.cc b/test/multidomain/boundary/stokesdarcy/1p_2p/main.cc
index 400a22ae2696124345a4974ee32207b430c2d9fa..1e0c171ba8af99fd529617ecaac48e0b07749c03 100644
--- a/test/multidomain/boundary/stokesdarcy/1p_2p/main.cc
+++ b/test/multidomain/boundary/stokesdarcy/1p_2p/main.cc
@@ -218,11 +218,11 @@ int main(int argc, char** argv) try
     // the assembler with time loop for instationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
                                                                  darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 stokesGridVariables->faceGridVariablesPtr(),
+                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
+                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
                                                                  darcyGridVariables),
                                                  couplingManager,
                                                  timeLoop, solOld);