diff --git a/dumux/discretization/fvgridvariables.hh b/dumux/discretization/fvgridvariables.hh
index 039fe29d42a7de506562f6a1681922e0b4e3f5e4..ab49a56ee83b866e3812e3bc289b450ae3e9197b 100644
--- a/dumux/discretization/fvgridvariables.hh
+++ b/dumux/discretization/fvgridvariables.hh
@@ -63,7 +63,7 @@ public:
 
     template<class Problem>
     FVGridVariables(std::shared_ptr<Problem> problem,
-                    std::shared_ptr<GridGeometry> fvGridGeometry)
+                    std::shared_ptr<const GridGeometry> fvGridGeometry)
     : fvGridGeometry_(fvGridGeometry)
     , curGridVolVars_(*problem)
     , prevGridVolVars_(*problem)
diff --git a/dumux/geomechanics/poroelastic/couplingmanager.hh b/dumux/geomechanics/poroelastic/couplingmanager.hh
index 0c42b5265eef6d99a300ad4003b1c05ca97d3ae8..de49960670cc1251a2e41ce5beeb715b62711e86 100644
--- a/dumux/geomechanics/poroelastic/couplingmanager.hh
+++ b/dumux/geomechanics/poroelastic/couplingmanager.hh
@@ -58,7 +58,7 @@ class PoroMechanicsCouplingManager : public virtual CouplingManager< MDTraits >
     using ParentType = CouplingManager< MDTraits >;
 
     // the sub-domain type tags
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
     // further types specific to the sub-problems
     template<std::size_t id> using Scalar = GetPropType<SubDomainTypeTag<id>, Properties::Scalar>;
diff --git a/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh
index d2afc37ea074014e4fe6498d83d65d25b857fad8..6cc54be3ae4050506b8b46ed98265bc41fa34886 100644
--- a/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh
+++ b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh
@@ -58,13 +58,13 @@ class DarcyDarcyBoundaryCouplingManager
     using Scalar = typename MDTraits::Scalar;
     using SolutionVector = typename MDTraits::SolutionVector;
 
-    template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<i>;
+    template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomain<i>::TypeTag;
     template<std::size_t i> using Problem = GetPropType<SubDomainTypeTag<i>, Properties::Problem>;
     template<std::size_t i> using PrimaryVariables = GetPropType<SubDomainTypeTag<i>, Properties::PrimaryVariables>;
     template<std::size_t i> using NumEqVector = GetPropType<SubDomainTypeTag<i>, Properties::NumEqVector>;
     template<std::size_t i> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<i>, Properties::GridVolumeVariables>::LocalView;
     template<std::size_t i> using VolumeVariables = typename GetPropType<SubDomainTypeTag<i>, Properties::GridVolumeVariables>::VolumeVariables;
-    template<std::size_t i> using FVGridGeometry = typename MDTraits::template SubDomainFVGridGeometry<i>;
+    template<std::size_t i> using FVGridGeometry = typename MDTraits::template SubDomain<i>::FVGridGeometry;
     template<std::size_t i> using FVElementGeometry = typename FVGridGeometry<i>::LocalView;
     template<std::size_t i> using SubControlVolumeFace = typename FVGridGeometry<i>::SubControlVolumeFace;
     template<std::size_t i> using SubControlVolume = typename FVGridGeometry<i>::SubControlVolume;
@@ -73,7 +73,7 @@ class DarcyDarcyBoundaryCouplingManager
 
     template<std::size_t i>
     static constexpr auto domainIdx()
-    { return typename MDTraits::template DomainIdx<i>{}; }
+    { return typename MDTraits::template SubDomain<i>::Index{}; }
 
     template<std::size_t i>
     static constexpr bool isCCTpfa()
diff --git a/dumux/multidomain/boundary/darcydarcy/couplingmapper.hh b/dumux/multidomain/boundary/darcydarcy/couplingmapper.hh
index 6ddda7ed9d9800dcf840f03c8510386517cd5680..f22e8a1f1ca118c8d1449d52ab6bbaa4849ecc13 100644
--- a/dumux/multidomain/boundary/darcydarcy/couplingmapper.hh
+++ b/dumux/multidomain/boundary/darcydarcy/couplingmapper.hh
@@ -53,14 +53,14 @@ class DarcyDarcyBoundaryCouplingMapper
 {
     using Scalar = typename MDTraits::Scalar;
 
-    template<std::size_t i> using FVGridGeometry = typename MDTraits::template SubDomainFVGridGeometry<i>;
+    template<std::size_t i> using FVGridGeometry = typename MDTraits::template SubDomain<i>::FVGridGeometry;
     template<std::size_t i> using SubControlVolumeFace = typename FVGridGeometry<i>::SubControlVolumeFace;
     template<std::size_t i> using GridView = typename FVGridGeometry<i>::GridView;
     template<std::size_t i> using Element = typename GridView<i>::template Codim<0>::Entity;
 
     template<std::size_t i>
     static constexpr auto domainIdx()
-    { return typename MDTraits::template DomainIdx<i>{}; }
+    { return typename MDTraits::template SubDomain<i>::Index{}; }
 
     template<std::size_t i>
     static constexpr bool isCCTpfa()
diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
index b7f12b80202bdaab5de5b31e256337881ce4cd64..dd03b033768d25990ae3fd84ab6b233ccfc004a1 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
@@ -217,8 +217,8 @@ class StokesDarcyCouplingDataImplementation;
 */
 template<class MDTraits, class CouplingManager>
 using StokesDarcyCouplingData = StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager,
-                                                                      GetPropType<typename MDTraits::template SubDomainTypeTag<0>, Properties::ModelTraits>::enableEnergyBalance(),
-                                                                      (GetPropType<typename MDTraits::template SubDomainTypeTag<0>, Properties::ModelTraits>::numComponents() > 1)>;
+                                                                      GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::enableEnergyBalance(),
+                                                                      (GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::numComponents() > 1)>;
 
 /*!
  * \ingroup MultiDomain
@@ -231,7 +231,7 @@ class StokesDarcyCouplingDataImplementationBase
 {
     using Scalar = typename MDTraits::Scalar;
 
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
     template<std::size_t id> using FVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
     template<std::size_t id> using Element = typename FVGridGeometry<id>::GridView::template Codim<0>::Entity;
     template<std::size_t id> using FVElementGeometry = typename FVGridGeometry<id>::LocalView;
@@ -469,14 +469,14 @@ class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEne
 {
     using ParentType = StokesDarcyCouplingDataImplementationBase<MDTraits, CouplingManager>;
     using Scalar = typename MDTraits::Scalar;
-    static constexpr auto stokesIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto darcyIdx = typename MDTraits::template DomainIdx<2>();
+    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 DomainIdx<1>();
+    static constexpr auto stokesFaceIdx = typename MDTraits::template SubDomain<1>::Index();
 
     // the sub domain type tags
     template<std::size_t id>
-    using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
     template<std::size_t id> using FVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
     template<std::size_t id> using Element = typename FVGridGeometry<id>::GridView::template Codim<0>::Entity;
@@ -633,14 +633,14 @@ class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEne
 {
     using ParentType = StokesDarcyCouplingDataImplementationBase<MDTraits, CouplingManager>;
     using Scalar = typename MDTraits::Scalar;
-    static constexpr auto stokesIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto darcyIdx = typename MDTraits::template DomainIdx<2>();
+    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 DomainIdx<1>();
+    static constexpr auto stokesFaceIdx = typename MDTraits::template SubDomain<1>::Index();
 
     // the sub domain type tags
     template<std::size_t id>
-    using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
     template<std::size_t id> using FVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
     template<std::size_t id> using Element = typename FVGridGeometry<id>::GridView::template Codim<0>::Entity;
diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh b/dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh
index df2e4d1e4ec218ef5f78c2a933d0f571b51f7f1d..16ed0624c1b041d26cdceb27ceb4a7069e111de1 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh
@@ -55,27 +55,27 @@ class StokesDarcyCouplingManager
     using ParentType = StaggeredCouplingManagerBase<MDTraits, StokesDarcyCouplingManager<MDTraits>>;
 
 public:
-    static constexpr auto stokesCellCenterIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto stokesFaceIdx = typename MDTraits::template DomainIdx<1>();
-    static constexpr auto cellCenterIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto faceIdx = typename MDTraits::template DomainIdx<1>();
+    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 DomainIdx<2>();
+    static constexpr auto darcyIdx = typename MDTraits::template SubDomain<2>::Index();
 
 private:
 
     using SolutionVector = typename MDTraits::SolutionVector;
 
     // obtain the type tags of the sub problems
-    using StokesTypeTag = typename MDTraits::template SubDomainTypeTag<0>;
-    using DarcyTypeTag = typename MDTraits::template SubDomainTypeTag<2>;
+    using StokesTypeTag = typename MDTraits::template SubDomain<0>::TypeTag;
+    using DarcyTypeTag = typename MDTraits::template SubDomain<2>::TypeTag;
 
     using CouplingStencils = std::unordered_map<std::size_t, std::vector<std::size_t> >;
     using CouplingStencil = CouplingStencils::mapped_type;
 
     // the sub domain type tags
     template<std::size_t id>
-    using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
     static constexpr bool isCompositional = GetPropType<SubDomainTypeTag<0>, Properties::ModelTraits>::numComponents() > 1;
 
@@ -91,7 +91,7 @@ private:
     template<std::size_t id> using ElementFluxVariablesCache = typename GetPropType<SubDomainTypeTag<id>, Properties::GridFluxVariablesCache>::LocalView;
     template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
     template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
-    template<std::size_t id> using PrimaryVariables = typename MDTraits::template PrimaryVariables<id>;
+    template<std::size_t id> using PrimaryVariables = typename MDTraits::template SubDomain<id>::PrimaryVariables;
     template<std::size_t id> using SubControlVolumeFace  = typename FVElementGeometry<id>::SubControlVolumeFace;
 
     using CellCenterSolutionVector = GetPropType<StokesTypeTag, Properties::CellCenterSolutionVector>;
diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingmapper.hh b/dumux/multidomain/boundary/stokesdarcy/couplingmapper.hh
index ea9ae18f42d0e7c5b215fe78b4754667471c3ac0..b6d68ef6065c1f52c2289790622fb7656a548c6c 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingmapper.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingmapper.hh
@@ -54,18 +54,18 @@ class StokesDarcyCouplingMapper
     using Scalar = typename MDTraits::Scalar;
 
 public:
-    static constexpr auto stokesCellCenterIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto stokesFaceIdx = typename MDTraits::template DomainIdx<1>();
-    static constexpr auto cellCenterIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto faceIdx = typename MDTraits::template DomainIdx<1>();
+    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 DomainIdx<2>();
+    static constexpr auto darcyIdx = typename MDTraits::template SubDomain<2>::Index();
 
 
 private:
     // obtain the type tags of the sub problems
-    using StokesTypeTag = typename MDTraits::template SubDomainTypeTag<0>;
-    using DarcyTypeTag = typename MDTraits::template SubDomainTypeTag<2>;
+    using StokesTypeTag = typename MDTraits::template SubDomain<0>::TypeTag;
+    using DarcyTypeTag = typename MDTraits::template SubDomain<2>::TypeTag;
 
     struct ElementMapInfo
     {
@@ -76,7 +76,7 @@ private:
 
     // the sub domain type tags
     template<std::size_t id>
-    using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
     using CouplingManager = GetPropType<StokesTypeTag, Properties::CouplingManager>;
 
     static_assert(GetPropType<SubDomainTypeTag<stokesIdx>, Properties::FVGridGeometry>::discMethod == DiscretizationMethod::staggered,
diff --git a/dumux/multidomain/couplingmanager.hh b/dumux/multidomain/couplingmanager.hh
index 9e220aa31a359133729e0021046e74afaf57ff13..37f2b20f5f8c2a4947596547d7d1943b3f97f803 100644
--- a/dumux/multidomain/couplingmanager.hh
+++ b/dumux/multidomain/couplingmanager.hh
@@ -44,13 +44,13 @@ namespace Dumux {
 template<class Traits>
 class CouplingManager
 {
-    template<std::size_t id> using SubDomainTypeTag = typename Traits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename Traits::template SubDomain<id>::TypeTag;
     template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
     template<std::size_t id> using GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>::GridView;
     template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
     template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
     template<std::size_t id> using ProblemWeakPtr = std::weak_ptr<const Problem<id>>;
-    using Problems = typename Traits::template MakeTuple<ProblemWeakPtr>;
+    using Problems = typename Traits::template Tuple<ProblemWeakPtr>;
 
 public:
     //! default type used for coupling element stencils
diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh
index 41b0ae8eeaab3ba739399658fe40d51737d6461b..c7a4b65e65c9d2da5ab4a6a73cd9afdd0f6043ad 100644
--- a/dumux/multidomain/embedded/couplingmanager1d3d.hh
+++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh
@@ -54,7 +54,7 @@ template<class MDTraits>
 struct CircleAveragePointSourceTraits
 {
 private:
-    template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<i>;
+    template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomain<i>::TypeTag;
     template<std::size_t i> using FVGridGeometry = GetPropType<SubDomainTypeTag<i>, Properties::FVGridGeometry>;
     template<std::size_t i> using NumEqVector = GetPropType<SubDomainTypeTag<i>, Properties::NumEqVector>;
 public:
@@ -98,11 +98,11 @@ class EmbeddedCouplingManager1d3d<MDTraits, EmbeddedCouplingMode::line>
     using Scalar = typename MDTraits::Scalar;
     using SolutionVector = typename MDTraits::SolutionVector;
 
-    static constexpr auto bulkIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto lowDimIdx = typename MDTraits::template DomainIdx<1>();
+    static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
+    static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
 
     // the sub domain type aliases
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
     template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
     template<std::size_t id> using FVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
     template<std::size_t id> using GridView = typename FVGridGeometry<id>::GridView;
@@ -198,11 +198,11 @@ class EmbeddedCouplingManager1d3d<MDTraits, EmbeddedCouplingMode::average>
     using SolutionVector = typename MDTraits::SolutionVector;
     using PointSourceData = typename ParentType::PointSourceTraits::PointSourceData;
 
-    static constexpr auto bulkIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto lowDimIdx = typename MDTraits::template DomainIdx<1>();
+    static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
+    static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
 
     // the sub domain type aliases
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
     template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
     template<std::size_t id> using FVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
     template<std::size_t id> using GridView = typename FVGridGeometry<id>::GridView;
@@ -576,11 +576,11 @@ class EmbeddedCouplingManager1d3d<MDTraits, EmbeddedCouplingMode::cylindersource
     using SolutionVector = typename MDTraits::SolutionVector;
     using PointSourceData = typename ParentType::PointSourceTraits::PointSourceData;
 
-    static constexpr auto bulkIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto lowDimIdx = typename MDTraits::template DomainIdx<1>();
+    static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
+    static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
 
     // the sub domain type aliases
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
     template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
     template<std::size_t id> using FVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
     template<std::size_t id> using GridView = typename FVGridGeometry<id>::GridView;
@@ -853,11 +853,11 @@ class EmbeddedCouplingManager1d3d<MDTraits, EmbeddedCouplingMode::kernel>
     using SolutionVector = typename MDTraits::SolutionVector;
     using PointSourceData = typename ParentType::PointSourceTraits::PointSourceData;
 
-    static constexpr auto bulkIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto lowDimIdx = typename MDTraits::template DomainIdx<1>();
+    static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
+    static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
 
     // the sub domain type aliases
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
     template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
     template<std::size_t id> using FVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
     template<std::size_t id> using GridView = typename FVGridGeometry<id>::GridView;
diff --git a/dumux/multidomain/embedded/couplingmanagerbase.hh b/dumux/multidomain/embedded/couplingmanagerbase.hh
index fa5fec57f1a3b56a12cd9316dfdde495d25a93ed..c56ad8e041a1e10e95785e6654805a1e0166c395 100644
--- a/dumux/multidomain/embedded/couplingmanagerbase.hh
+++ b/dumux/multidomain/embedded/couplingmanagerbase.hh
@@ -51,7 +51,7 @@ template<class MDTraits>
 struct DefaultPointSourceTraits
 {
 private:
-    template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<i>;
+    template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomain<i>::TypeTag;
     template<std::size_t i> using FVGridGeometry = GetPropType<SubDomainTypeTag<i>, Properties::FVGridGeometry>;
     template<std::size_t i> using NumEqVector = GetPropType<SubDomainTypeTag<i>, Properties::NumEqVector>;
 public:
@@ -79,14 +79,14 @@ class EmbeddedCouplingManagerBase
 {
     using ParentType = CouplingManager<MDTraits>;
     using Scalar = typename MDTraits::Scalar;
-    static constexpr auto bulkIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto lowDimIdx = typename MDTraits::template DomainIdx<1>();
+    static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
+    static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
     using SolutionVector = typename MDTraits::SolutionVector;
     using PointSourceData = typename PSTraits::PointSourceData;
 
     // the sub domain type tags
     template<std::size_t id> using PointSource = typename PSTraits::template PointSource<id>;
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
     template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
     template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
     template<std::size_t id> using FVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
diff --git a/dumux/multidomain/embedded/extendedsourcestencil.hh b/dumux/multidomain/embedded/extendedsourcestencil.hh
index ead7fb11a7a4117e1075689471edd16ee668bf6e..ad266597aaf8a22400942273c503a0105298634b 100644
--- a/dumux/multidomain/embedded/extendedsourcestencil.hh
+++ b/dumux/multidomain/embedded/extendedsourcestencil.hh
@@ -50,13 +50,13 @@ class ExtendedSourceStencil
     using MDTraits = typename CouplingManager::MultiDomainTraits;
     using Scalar = typename MDTraits::Scalar;
 
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
     template<std::size_t id> using FVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
     template<std::size_t id> using GridView = typename FVGridGeometry<id>::GridView;
     template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
 
-    static constexpr auto bulkIdx = typename MDTraits::template DomainIdx<0>();
-    static constexpr auto lowDimIdx = typename MDTraits::template DomainIdx<1>();
+    static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
+    static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
 
     template<std::size_t id>
     static constexpr bool isBox()
diff --git a/dumux/multidomain/embedded/pointsourcedata.hh b/dumux/multidomain/embedded/pointsourcedata.hh
index 088defdf6704c2aff382c1c8b5675eaf31794d5c..4bade7c0a900a0900e59b9040ba51ead970026f5 100644
--- a/dumux/multidomain/embedded/pointsourcedata.hh
+++ b/dumux/multidomain/embedded/pointsourcedata.hh
@@ -47,8 +47,8 @@ class PointSourceData
     using ShapeValues = typename std::vector<Dune::FieldVector<Scalar, 1> >;
 
     // obtain the type tags of the sub problems
-    using BulkTypeTag = typename MDTraits::template SubDomainTypeTag<0>;
-    using LowDimTypeTag = typename MDTraits::template SubDomainTypeTag<1>;
+    using BulkTypeTag = typename MDTraits::template SubDomain<0>::TypeTag;
+    using LowDimTypeTag = typename MDTraits::template SubDomain<1>::TypeTag;
 
     using BulkPrimaryVariables = GetPropType<BulkTypeTag, Properties::PrimaryVariables>;
     using LowDimPrimaryVariables = GetPropType<LowDimTypeTag, Properties::PrimaryVariables>;
@@ -155,8 +155,8 @@ class PointSourceDataCircleAverage : public PointSourceData<MDTraits>
     using ShapeValues = typename std::vector<Dune::FieldVector<Scalar, 1> >;
 
     // obtain the type tags of the sub problems
-    using BulkTypeTag = typename MDTraits::template SubDomainTypeTag<0>;
-    using LowDimTypeTag = typename MDTraits::template SubDomainTypeTag<1>;
+    using BulkTypeTag = typename MDTraits::template SubDomain<0>::TypeTag;
+    using LowDimTypeTag = typename MDTraits::template SubDomain<1>::TypeTag;
 
     using BulkPrimaryVariables = GetPropType<BulkTypeTag, Properties::PrimaryVariables>;
     using LowDimPrimaryVariables = GetPropType<LowDimTypeTag, Properties::PrimaryVariables>;
diff --git a/dumux/multidomain/facet/box/couplingmanager.hh b/dumux/multidomain/facet/box/couplingmanager.hh
index 6c5f2b9acd741240c696eec20311090d9f5ad883..c725306a2d4cfde725e46f9213e3e62c4e53c217 100644
--- a/dumux/multidomain/facet/box/couplingmanager.hh
+++ b/dumux/multidomain/facet/box/couplingmanager.hh
@@ -55,13 +55,13 @@ class FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, lowDimDomainI
     using ParentType = CouplingManager< MDTraits >;
 
     // convenience aliases and instances of the two domain ids
-    using BulkIdType = typename MDTraits::template DomainIdx<bulkDomainId>;
-    using LowDimIdType = typename MDTraits::template DomainIdx<lowDimDomainId>;
+    using BulkIdType = typename MDTraits::template SubDomain<bulkDomainId>::Index;
+    using LowDimIdType = typename MDTraits::template SubDomain<lowDimDomainId>::Index;
     static constexpr auto bulkId = BulkIdType();
     static constexpr auto lowDimId = LowDimIdType();
 
     // the sub-domain type tags
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
     // further types specific to the sub-problems
     template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
diff --git a/dumux/multidomain/facet/cellcentered/tpfa/couplingmanager.hh b/dumux/multidomain/facet/cellcentered/tpfa/couplingmanager.hh
index bc9bd8b214508f555968be3a4d9b059ed4b95890..b0b80259755b720fe8ff42a3b18a004e14b9e1be 100644
--- a/dumux/multidomain/facet/cellcentered/tpfa/couplingmanager.hh
+++ b/dumux/multidomain/facet/cellcentered/tpfa/couplingmanager.hh
@@ -55,13 +55,13 @@ class FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, lowDimDomainI
     using ParentType = CouplingManager< MDTraits >;
 
     // convenience aliases and instances of the two domain ids
-    using BulkIdType = typename MDTraits::template DomainIdx<bulkDomainId>;
-    using LowDimIdType = typename MDTraits::template DomainIdx<lowDimDomainId>;
+    using BulkIdType = typename MDTraits::template SubDomain<bulkDomainId>::Index;
+    using LowDimIdType = typename MDTraits::template SubDomain<lowDimDomainId>::Index;
     static constexpr auto bulkId = BulkIdType();
     static constexpr auto lowDimId = LowDimIdType();
 
     // the sub-domain type tags
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
     // further types specific to the sub-problems
     template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
diff --git a/dumux/multidomain/facet/couplingmanager.hh b/dumux/multidomain/facet/couplingmanager.hh
index 11c8261594ac4e4d67277f51910b31328ad36887..4668b2cefe4cadad4d645fb807796442afcd3999 100644
--- a/dumux/multidomain/facet/couplingmanager.hh
+++ b/dumux/multidomain/facet/couplingmanager.hh
@@ -90,7 +90,7 @@ template< class MDTraits,
           class CouplingMapper,
           std::size_t bulkDomainId = 0,
           std::size_t lowDimDomainId = 1,
-          DiscretizationMethod bulkDM = GetPropType<typename MDTraits::template SubDomainTypeTag<bulkDomainId>, Properties::FVGridGeometry>::discMethod >
+          DiscretizationMethod bulkDM = GetPropType<typename MDTraits::template SubDomain<bulkDomainId>::TypeTag, Properties::FVGridGeometry>::discMethod >
 class FacetCouplingManager;
 
 /*!
@@ -119,15 +119,15 @@ class FacetCouplingThreeDomainManager
     using FacetEdgeManager = FacetCouplingManager<MDTraits, CouplingMapper, facetDomainId, edgeDomainId>;
 
     // convenience aliases and instances of the domain ids
-    using BulkIdType = typename MDTraits::template DomainIdx<bulkDomainId>;
-    using FacetIdType = typename MDTraits::template DomainIdx<facetDomainId>;
-    using EdgeIdType = typename MDTraits::template DomainIdx<edgeDomainId>;
+    using BulkIdType = typename MDTraits::template SubDomain<bulkDomainId>::Index;
+    using FacetIdType = typename MDTraits::template SubDomain<facetDomainId>::Index;
+    using EdgeIdType = typename MDTraits::template SubDomain<edgeDomainId>::Index;
     static constexpr auto bulkId = BulkIdType();
     static constexpr auto facetId = FacetIdType();
     static constexpr auto edgeId = EdgeIdType();
 
     // the sub-domain type tags
-    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
     // further types specific to the sub-problems
     template<std::size_t id> using LocalResidual = GetPropType<SubDomainTypeTag<id>, Properties::LocalResidual>;
diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh
index 13e80fa458b1e8edd00e2b90c5d639096670954b..dfd40d2decc5c704687356fcdf966aedac1839dd 100644
--- a/dumux/multidomain/fvassembler.hh
+++ b/dumux/multidomain/fvassembler.hh
@@ -58,18 +58,25 @@ template<class MDTraits, class CMType, DiffMethod diffMethod, bool useImplicitAs
 class MultiDomainFVAssembler
 {
     template<std::size_t id>
-    using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
 public:
     using Traits = MDTraits;
 
     using Scalar = typename MDTraits::Scalar;
 
+    //! TODO get rid of this GetPropType
     template<std::size_t id>
     using LocalResidual = GetPropType<SubDomainTypeTag<id>, Properties::LocalResidual>;
 
     template<std::size_t id>
-    using GridVariables = typename std::tuple_element_t<id, typename MDTraits::GridVariablesTuple>::element_type;
+    using GridVariables = typename MDTraits::template SubDomain<id>::GridVariables;
+
+    template<std::size_t id>
+    using FVGridGeometry = typename MDTraits::template SubDomain<id>::FVGridGeometry;
+
+    template<std::size_t id>
+    using Problem = typename MDTraits::template SubDomain<id>::Problem;
 
     using JacobianMatrix = typename MDTraits::JacobianMatrix;
     using SolutionVector = typename MDTraits::SolutionVector;
@@ -85,9 +92,9 @@ public:
 
 private:
 
-    using ProblemTuple = typename MDTraits::ProblemTuple;
-    using FVGridGeometryTuple = typename MDTraits::FVGridGeometryTuple;
-    using GridVariablesTuple = typename MDTraits::GridVariablesTuple;
+    using ProblemTuple = typename MDTraits::template TupleOfSharedPtrConst<Problem>;
+    using FVGridGeometryTuple = typename MDTraits::template TupleOfSharedPtrConst<FVGridGeometry>;
+    using GridVariablesTuple = typename MDTraits::template TupleOfSharedPtr<GridVariables>;
 
     using TimeLoop = TimeLoopBase<Scalar>;
     using ThisType = MultiDomainFVAssembler<MDTraits, CouplingManager, diffMethod, isImplicit()>;
@@ -119,9 +126,6 @@ private:
         using type = SubDomainStaggeredLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
     };
 
-    template<std::size_t id>
-    using FVGridGeometry = typename std::tuple_element<id, FVGridGeometryTuple>::type::element_type;
-
     template<std::size_t id>
     using SubDomainAssembler = typename SubDomainAssemblerType<FVGridGeometry<id>::discMethod, id>::type;
 
diff --git a/dumux/multidomain/newtonsolver.hh b/dumux/multidomain/newtonsolver.hh
index 2f3b2a987cac094720214fa16236c47ea87c96e2..8b2b4aa0e42770dde5304186547c96498dfb134b 100644
--- a/dumux/multidomain/newtonsolver.hh
+++ b/dumux/multidomain/newtonsolver.hh
@@ -60,7 +60,7 @@ class MultiDomainNewtonSolver: public NewtonSolver<Assembler, LinearSolver, Reas
 
     template<std::size_t i>
     using PrivarSwitchPtr = std::unique_ptr<PrimaryVariableSwitch<i>>;
-    using PriVarSwitchPtrTuple = typename Assembler::Traits::template MakeTuple<PrivarSwitchPtr>;
+    using PriVarSwitchPtrTuple = typename Assembler::Traits::template Tuple<PrivarSwitchPtr>;
 
 public:
 
diff --git a/dumux/multidomain/staggeredcouplingmanager.hh b/dumux/multidomain/staggeredcouplingmanager.hh
index 60c90bdd0af23498d304903a938b3d14566a97cd..c463055362176b3c083d7da8e009e32fa1ef3d86 100644
--- a/dumux/multidomain/staggeredcouplingmanager.hh
+++ b/dumux/multidomain/staggeredcouplingmanager.hh
@@ -42,12 +42,12 @@ class StaggeredCouplingManagerBase: public CouplingManager<MDTraits>
 {
     using ParentType = CouplingManager<MDTraits>;
     template<std::size_t id>
-    using SubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<id>;
+    using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
     template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
 
-    using StaggeredSubDomainTypeTag = typename MDTraits::template SubDomainTypeTag<0>;
+    using StaggeredSubDomainTypeTag = typename MDTraits::template SubDomain<0>::TypeTag;
 
-    template<std::size_t id> using FVGridGeometry = typename std::tuple_element_t<id, typename MDTraits::FVGridGeometryTuple>::element_type;
+    template<std::size_t id> using FVGridGeometry = typename MDTraits::template SubDomain<id>::FVGridGeometry;
     template<std::size_t id> using GridView = typename FVGridGeometry<id>::GridView;
 
     using FVElementGeometry = typename FVGridGeometry<0>::LocalView;
@@ -205,8 +205,8 @@ public:
     decltype(auto) numericEpsilon(Dune::index_constant<i>,
                                   const std::string& paramGroup) const
     {
-        constexpr std::size_t numEqCellCenter = Traits::template PrimaryVariables<cellCenterIdx>::dimension;
-        constexpr std::size_t numEqFace = Traits::template PrimaryVariables<faceIdx>::dimension;
+        constexpr std::size_t numEqCellCenter = Traits::template SubDomain<cellCenterIdx>::PrimaryVariables::dimension;
+        constexpr std::size_t numEqFace = Traits::template SubDomain<faceIdx>::PrimaryVariables::dimension;
         constexpr bool isCellCenter = FVGridGeometry<i>::isCellCenter();
         constexpr std::size_t numEq = isCellCenter ? numEqCellCenter : numEqFace;
         constexpr auto prefix = isCellCenter ? "CellCenter" : "Face";
diff --git a/dumux/multidomain/staggeredtraits.hh b/dumux/multidomain/staggeredtraits.hh
index 96e9de31cd5a28dc0eaaeb5cd599ffaed8d23e7f..66b3af03f5e94e88b7370e21692b9ee6ed8a8d5c 100644
--- a/dumux/multidomain/staggeredtraits.hh
+++ b/dumux/multidomain/staggeredtraits.hh
@@ -43,6 +43,92 @@
 #include "traits.hh"
 
 namespace Dumux {
+namespace Detail {
+namespace Staggered {
+
+//////////////////////////////////////////////////////////
+template<template<std::size_t> class SubDomainTypeTag, std::size_t i>
+struct SubDomainFVGridGeometryImpl
+{ using type = GetPropType<SubDomainTypeTag<i>, Properties::FVGridGeometry>; };
+
+template<template<std::size_t> class SubDomainTypeTag>
+struct SubDomainFVGridGeometryImpl<SubDomainTypeTag, 0>
+{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::FVGridGeometry>::CellCenterFVGridGeometryType; };
+
+template<template<std::size_t> class SubDomainTypeTag>
+struct SubDomainFVGridGeometryImpl<SubDomainTypeTag, 1>
+{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::FVGridGeometry>::FaceFVGridGeometryType; };
+//////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////
+template<template<std::size_t> class SubDomainTypeTag, std::size_t i>
+struct SubDomainGridVariablesImpl
+{ using type = GetPropType<SubDomainTypeTag<i>, Properties::GridVariables>; };
+
+template<template<std::size_t> class SubDomainTypeTag>
+struct SubDomainGridVariablesImpl<SubDomainTypeTag, 0>
+{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::GridVariables>::CellCenterGridVariablesType; };
+
+template<template<std::size_t> class SubDomainTypeTag>
+struct SubDomainGridVariablesImpl<SubDomainTypeTag, 1>
+{ using type = typename GetPropType<SubDomainTypeTag<0>, Properties::GridVariables>::FaceGridVariablesType; };
+//////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////
+template<template<std::size_t> class SubDomainTypeTag, std::size_t i>
+struct SubDomainPrimaryVariablesImpl
+{ using type = GetPropType<SubDomainTypeTag<i>, Properties::PrimaryVariables>; };
+
+template<template<std::size_t> class SubDomainTypeTag>
+struct SubDomainPrimaryVariablesImpl<SubDomainTypeTag, 0>
+{ using type = GetPropType<SubDomainTypeTag<0>, Properties::CellCenterPrimaryVariables>; };
+
+template<template<std::size_t> class SubDomainTypeTag>
+struct SubDomainPrimaryVariablesImpl<SubDomainTypeTag, 1>
+{ using type = GetPropType<SubDomainTypeTag<0>, Properties::FacePrimaryVariables>; };
+//////////////////////////////////////////////////////////
+
+template<class Scalar, int numEq>
+struct JacobianTypeImpl
+{
+    private:
+        using MatrixBlock = typename Dune::FieldMatrix<Scalar, numEq, numEq>;
+    public:
+        using type = typename Dune::BCRSMatrix<MatrixBlock>;
+};
+
+//////////////////////////////////////////////////////////
+template<template<std::size_t> class SubDomainTypeTag, std::size_t i>
+struct SubDomainJacobianMatrixImpl
+{ using type = GetPropType<SubDomainTypeTag<i>, Properties::JacobianMatrix>; };
+
+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; };
+
+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; };
+//////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////
+template<template<std::size_t> class SubDomainTypeTag, std::size_t i>
+struct SubDomainSolutionVectorImpl
+{ using type = GetPropType<SubDomainTypeTag<i>, Properties::SolutionVector>; };
+
+template<template<std::size_t> class SubDomainTypeTag>
+struct SubDomainSolutionVectorImpl<SubDomainTypeTag, 0>
+{ using type = GetPropType<SubDomainTypeTag<0>, Properties::CellCenterSolutionVector>; };
+
+template<template<std::size_t> class SubDomainTypeTag>
+struct SubDomainSolutionVectorImpl<SubDomainTypeTag, 1>
+{ using type = GetPropType<SubDomainTypeTag<0>, Properties::FaceSolutionVector>; };
+//////////////////////////////////////////////////////////
+
+} // end namespace Staggered
+} // end namespace Detail
 
 /*
  * \ingroup MultiDomain
@@ -73,95 +159,83 @@ namespace Dumux {
 template<typename... SubDomainTypeTags>
 struct StaggeredMultiDomainTraits
 {
+    //! the number of subdomains
     static constexpr std::size_t numSubDomains = sizeof...(SubDomainTypeTags);
 
+private:
+
     //! the type tag of a sub domain problem
     template<std::size_t id>
     using SubDomainTypeTag = typename std::tuple_element_t<id, std::tuple<SubDomainTypeTags...>>;
 
-    //! the static domain indices
-    template<std::size_t id>
-    using DomainIdx = Dune::index_constant<id>;
-
-private:
+    //! helper alias to construct derived multidomain types like tuples
     using Indices = std::make_index_sequence<numSubDomains>;
 
-    template<std::size_t id>
-    using SolutionSubVector = std::conditional_t<(id < 2),
-                                                 std::conditional_t<(id == 0),
-                                                                    GetPropType<SubDomainTypeTag<0>, Properties::CellCenterSolutionVector>,
-                                                                    GetPropType<SubDomainTypeTag<0>, Properties::FaceSolutionVector>>,
-                                                 GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>>;
-
+    //! the scalar type of each sub domain
     template<std::size_t id>
     using SubDomainScalar = GetPropType<SubDomainTypeTag<id>, Properties::Scalar>;
 
     template<std::size_t id>
-    using SubDomainProblem = std::shared_ptr<const GetPropType<SubDomainTypeTag<id>, Properties::Problem>>;
+    using SubDomainJacobianMatrix = typename Detail::Staggered::SubDomainJacobianMatrixImpl<SubDomainTypeTag, id>::type;
 
     template<std::size_t id>
-    using SubDomainFVGridGeometry = std::shared_ptr<std::conditional_t<(id < 2),
-                                                                       std::conditional_t<(id == 0),
-                                                                                          typename GetPropType<SubDomainTypeTag<0>, Properties::FVGridGeometry>::CellCenterFVGridGeometryType,
-                                                                                          typename GetPropType<SubDomainTypeTag<0>, Properties::FVGridGeometry>::FaceFVGridGeometryType>,
-                                                                       GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>>>;
+    using SubDomainSolutionVector = typename Detail::Staggered::SubDomainSolutionVectorImpl<SubDomainTypeTag, id>::type;
+
+public:
+
+    /*
+     * \brief sub domain types
+     */
+    //\{
 
     template<std::size_t id>
-    using SubDomainGridVariables = std::shared_ptr<std::conditional_t<(id < 2),
-                                                                      std::conditional_t<(id == 0),
-                                                                                         typename GetPropType<SubDomainTypeTag<0>, Properties::GridVariables>::CellCenterGridVariablesType,
-                                                                                         typename GetPropType<SubDomainTypeTag<0>, Properties::GridVariables>::FaceGridVariablesType>,
-                                                                      GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>>>;
-
-    template<class Scalar, int numEq>
-    struct JacobianType
+    struct SubDomain
     {
-        private:
-            using MatrixBlock = typename Dune::FieldMatrix<Scalar, numEq, numEq>;
-        public:
-            using type = typename Dune::BCRSMatrix<MatrixBlock>;
+        using Index = Dune::index_constant<id>;
+        using TypeTag = SubDomainTypeTag<id>;
+        using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
+        using FVGridGeometry = typename Detail::Staggered::SubDomainFVGridGeometryImpl<SubDomainTypeTag, id>::type;
+        using GridVariables = typename Detail::Staggered::SubDomainGridVariablesImpl<SubDomainTypeTag, id>::type;
+        using SolutionVector = typename Detail::Staggered::SubDomainSolutionVectorImpl<SubDomainTypeTag, id>::type;
+        using PrimaryVariables = typename Detail::Staggered::SubDomainPrimaryVariablesImpl<SubDomainTypeTag, id>::type;
     };
 
+    //\}
 
-    template<std::size_t id>
-    using JacobianDiagBlock = std::conditional_t<(id < 2),
-                                                 std::conditional_t<(id == 0),
-                                                                    typename JacobianType<GetPropType<SubDomainTypeTag<0>, Properties::Scalar>, getPropValue<SubDomainTypeTag<0>, Properties::NumEqCellCenter>()>::type,
-                                                                    typename JacobianType<GetPropType<SubDomainTypeTag<0>, Properties::Scalar>, getPropValue<SubDomainTypeTag<0>, Properties::NumEqFace>()>::type>,
-                                                 GetPropType<SubDomainTypeTag<id>, Properties::JacobianMatrix>>;
-
-public:
+    /*
+     * \brief multi domain types
+     */
+    //\{
 
     //! the scalar type
     using Scalar = typename makeFromIndexedType<std::common_type_t, SubDomainScalar, Indices>::type;
 
     //! the solution vector type
-    using SolutionVector = typename makeFromIndexedType<Dune::MultiTypeBlockVector, SolutionSubVector, Indices>::type;
-
-    template<std::size_t id>
-    using PrimaryVariables = std::conditional_t<(id < 2),
-                                                 std::conditional_t<(id == 0), GetPropType<SubDomainTypeTag<0>, Properties::CellCenterPrimaryVariables>,
-                                                                               GetPropType<SubDomainTypeTag<0>, Properties::FacePrimaryVariables>>,
-                                                 GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>>;
-
-    template<typename... MatrixBlocks>
-    using createMatrixType = typename Detail::createMultiTypeBlockMatrixType<Scalar, MatrixBlocks...>::type::type;
+    using SolutionVector = typename makeFromIndexedType<Dune::MultiTypeBlockVector, SubDomainSolutionVector, Indices>::type;
 
     //! the jacobian type
-    using JacobianMatrix = typename makeFromIndexedType<createMatrixType, JacobianDiagBlock, Indices>::type;
+    using JacobianMatrix = typename Detail::MultiDomainMatrixType<SubDomainJacobianMatrix, Indices, Scalar>::type;
 
-    //! the tuple of problems
-    using ProblemTuple = typename makeFromIndexedType<std::tuple, SubDomainProblem, Indices>::type;
+    //\}
 
-    //! the tuple of fv grid geometries
-    using FVGridGeometryTuple = typename makeFromIndexedType<std::tuple, SubDomainFVGridGeometry, Indices>::type;
+    /*
+     * \brief helper aliases to contruct derived tuple types
+     */
+    //\{
 
-    //! the tuple of grid variables
-    using GridVariablesTuple = typename makeFromIndexedType<std::tuple, SubDomainGridVariables, Indices>::type;
+    //! helper alias to create tuple<...> from indexed type
+    template<template<std::size_t> class T>
+    using Tuple = typename makeFromIndexedType<std::tuple, T, Indices>::type;
 
-    //! convenience alias to create tuple from type
+    //! helper alias to create tuple<std::shared_ptr<...>> from indexed type
     template<template<std::size_t> class T>
-    using MakeTuple = typename makeFromIndexedType<std::tuple, T, Indices>::type;
+    using TupleOfSharedPtr = typename Detail::MultiDomainTupleSharedPtr<T, Indices>::type;
+
+    //! helper alias to create tuple<std::shared_ptr<const ...>> from indexed type
+    template<template<std::size_t> class T>
+    using TupleOfSharedPtrConst = typename Detail::MultiDomainTupleSharedPtrConst<T, Indices>::type;
+
+    //\}
 };
 
 } //end namespace Dumux
diff --git a/dumux/multidomain/traits.hh b/dumux/multidomain/traits.hh
index bf8b7806d1eef3c949c45f4da459def06ebcc600..3c4b0633cdc85e3204cc6da5242099b563a14231 100644
--- a/dumux/multidomain/traits.hh
+++ b/dumux/multidomain/traits.hh
@@ -28,6 +28,7 @@
 #include <type_traits>
 #include <tuple>
 #include <utility>
+#include <memory>
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/indices.hh>
@@ -83,6 +84,36 @@ public:
     using type = typename makeMatrix<Indices>::type;
 };
 
+//! helper alias to create a tuple of shared_ptr<...> from an indexed type
+template<template<std::size_t> class T, class Indices>
+struct MultiDomainTupleSharedPtr
+{
+    template<std::size_t i>
+    using PtrType = std::shared_ptr<T<i>>;
+
+    using type = typename makeFromIndexedType<std::tuple, PtrType, Indices>::type;
+};
+
+//! helper alias to create a tuple of shared_ptr<const ...> from an indexed type
+template<template<std::size_t> class T, class Indices>
+struct MultiDomainTupleSharedPtrConst
+{
+    template<std::size_t i>
+    using PtrType = std::shared_ptr<const T<i>>;
+
+    using type = typename makeFromIndexedType<std::tuple, PtrType, Indices>::type;
+};
+
+//! helper alias to create the JacobianMatrix type
+template<template<std::size_t> class SubDomainDiagBlocks, class Indices, class Scalar>
+struct MultiDomainMatrixType
+{
+    template<typename... MatrixBlocks>
+    using M = typename createMultiTypeBlockMatrixType<Scalar, MatrixBlocks...>::type::type;
+
+    using type = typename makeFromIndexedType<M, SubDomainDiagBlocks, Indices>::type;
+};
+
 } // end namespace Detail
 
 /*
@@ -113,67 +144,86 @@ public:
 template<typename... SubDomainTypeTags>
 struct MultiDomainTraits
 {
+    //! the number of subdomains
     static constexpr std::size_t numSubDomains = sizeof...(SubDomainTypeTags);
 
+private:
+
     //! the type tag of a sub domain problem
     template<std::size_t id>
     using SubDomainTypeTag = typename std::tuple_element_t<id, std::tuple<SubDomainTypeTags...>>;
 
-    //! the static domain indices
-    template<std::size_t id>
-    using DomainIdx = Dune::index_constant<id>;
-
-    //! the sub domain geometry
-    template<std::size_t id>
-    using SubDomainFVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
-
-private:
+    //! helper alias to construct derived multidomain types like tuples
     using Indices = std::make_index_sequence<numSubDomains>;
 
-    template<std::size_t id>
-    using SolutionSubVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>;
-
+    //! the scalar type of each sub domain
     template<std::size_t id>
     using SubDomainScalar = GetPropType<SubDomainTypeTag<id>, Properties::Scalar>;
 
+    //! the jacobian type of each sub domain
     template<std::size_t id>
-    using SubDomainProblem = std::shared_ptr<const GetPropType<SubDomainTypeTag<id>, Properties::Problem>>;
+    using SubDomainJacobianMatrix = GetPropType<SubDomainTypeTag<id>, Properties::JacobianMatrix>;
 
+    //! the solution type of each sub domain
     template<std::size_t id>
-    using SubDomainFVGridGeometryPtr = std::shared_ptr<const GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>>;
+    using SubDomainSolutionVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>;
 
-    template<std::size_t id>
-    using SubDomainGridVariables = std::shared_ptr<GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>>;
+public:
+
+    /*
+     * \brief sub domain types
+     */
+    //\{
 
     template<std::size_t id>
-    using JacobianDiagBlock = GetPropType<SubDomainTypeTag<id>, Properties::JacobianMatrix>;
+    struct SubDomain
+    {
+        using Index = Dune::index_constant<id>;
+        using TypeTag = SubDomainTypeTag<id>;
+        using Grid = GetPropType<SubDomainTypeTag<id>, Properties::Grid>;
+        using FVGridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::FVGridGeometry>;
+        using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
+        using GridVariables =GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
+        using IOFields = GetPropType<SubDomainTypeTag<id>, Properties::IOFields>;
+        using SolutionVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>;
+    };
 
-public:
+    //\}
+
+    /*
+     * \brief multi domain types
+     */
+    //\{
 
     //! the scalar type
     using Scalar = typename makeFromIndexedType<std::common_type_t, SubDomainScalar, Indices>::type;
 
     //! the solution vector type
-    using SolutionVector = typename makeFromIndexedType<Dune::MultiTypeBlockVector, SolutionSubVector, Indices>::type;
-
-    template<typename... MatrixBlocks>
-    using createMatrixType = typename Detail::createMultiTypeBlockMatrixType<Scalar, MatrixBlocks...>::type::type;
+    using SolutionVector = typename makeFromIndexedType<Dune::MultiTypeBlockVector, SubDomainSolutionVector, Indices>::type;
 
     //! the jacobian type
-    using JacobianMatrix = typename makeFromIndexedType<createMatrixType, JacobianDiagBlock, Indices>::type;
+    using JacobianMatrix = typename Detail::MultiDomainMatrixType<SubDomainJacobianMatrix, Indices, Scalar>::type;
 
-    //! the tuple of problems
-    using ProblemTuple = typename makeFromIndexedType<std::tuple, SubDomainProblem, Indices>::type;
+    //\}
 
-    //! the tuple of fv grid geometries
-    using FVGridGeometryTuple = typename makeFromIndexedType<std::tuple, SubDomainFVGridGeometryPtr, Indices>::type;
+    /*
+     * \brief helper aliases to contruct derived tuple types
+     */
+    //\{
 
-    //! the tuple of grid variables
-    using GridVariablesTuple = typename makeFromIndexedType<std::tuple, SubDomainGridVariables, Indices>::type;
+    //! helper alias to create tuple<...> from indexed type
+    template<template<std::size_t> class T>
+    using Tuple = typename makeFromIndexedType<std::tuple, T, Indices>::type;
+
+    //! helper alias to create tuple<std::shared_ptr<...>> from indexed type
+    template<template<std::size_t> class T>
+    using TupleOfSharedPtr = typename Detail::MultiDomainTupleSharedPtr<T, Indices>::type;
 
-    //! convenience alias to create tuple from type
+    //! helper alias to create tuple<std::shared_ptr<const ...>> from indexed type
     template<template<std::size_t> class T>
-    using MakeTuple = typename makeFromIndexedType<std::tuple, T, Indices>::type;
+    using TupleOfSharedPtrConst = typename Detail::MultiDomainTupleSharedPtrConst<T, Indices>::type;
+
+    //\}
 };
 
 } //end namespace Dumux
diff --git a/test/multidomain/boundary/darcydarcy/1p_1p/main.cc b/test/multidomain/boundary/darcydarcy/1p_1p/main.cc
index a3015814e858d524dd9bab7708c3cc5706b18378..f17ef17321489a5821986213f822e56aeee78d70 100644
--- a/test/multidomain/boundary/darcydarcy/1p_1p/main.cc
+++ b/test/multidomain/boundary/darcydarcy/1p_1p/main.cc
@@ -195,8 +195,8 @@ int main(int argc, char** argv) try
 
     // the mixed dimension type traits
     using Traits = MultiDomainTraits<SubTypeTag0, SubTypeTag1>;
-    constexpr auto domain0Idx = Traits::template DomainIdx<0>();
-    constexpr auto domain1Idx = Traits::template DomainIdx<1>();
+    constexpr auto domain0Idx = Traits::template SubDomain<0>::Index();
+    constexpr auto domain1Idx = Traits::template SubDomain<1>::Index();
 
     // the coupling manager
     using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
diff --git a/test/multidomain/boundary/darcydarcy/1p_2p/main.cc b/test/multidomain/boundary/darcydarcy/1p_2p/main.cc
index 71884942c06a20c518e7e550b9491862c3ea96e7..ac0d68bb7f57fa496529608ccb6acdf3adfda5c1 100644
--- a/test/multidomain/boundary/darcydarcy/1p_2p/main.cc
+++ b/test/multidomain/boundary/darcydarcy/1p_2p/main.cc
@@ -181,8 +181,8 @@ int main(int argc, char** argv) try
 
     // the mixed dimension type traits
     using Traits = MultiDomainTraits<SubTypeTag0, SubTypeTag1>;
-    constexpr auto domain0Idx = Traits::template DomainIdx<0>();
-    constexpr auto domain1Idx = Traits::template DomainIdx<1>();
+    constexpr auto domain0Idx = Traits::template SubDomain<0>::Index();
+    constexpr auto domain1Idx = Traits::template SubDomain<1>::Index();
 
     // the coupling manager
     using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
diff --git a/test/multidomain/embedded/1d3d/1p2c_richards2c/main.cc b/test/multidomain/embedded/1d3d/1p2c_richards2c/main.cc
index 14000499b587f370b0491c3f505cc93e9250f078..3fcce4bade8ad2d943aa8d8c237431c5734d4400 100644
--- a/test/multidomain/embedded/1d3d/1p2c_richards2c/main.cc
+++ b/test/multidomain/embedded/1d3d/1p2c_richards2c/main.cc
@@ -214,8 +214,8 @@ int main(int argc, char** argv) try
 
     // the mixed dimension type traits
     using Traits = MultiDomainTraits<BulkTypeTag, LowDimTypeTag>;
-    constexpr auto bulkIdx = Traits::template DomainIdx<0>();
-    constexpr auto lowDimIdx = Traits::template DomainIdx<1>();
+    constexpr auto bulkIdx = Traits::template SubDomain<0>::Index();
+    constexpr auto lowDimIdx = Traits::template SubDomain<1>::Index();
 
     // the coupling manager
     using CouplingManager = GetPropType<BulkTypeTag, Properties::CouplingManager>;
diff --git a/test/multidomain/embedded/1d3d/1p_1p/main.cc b/test/multidomain/embedded/1d3d/1p_1p/main.cc
index 56126b6cf3f258166f1cecb0094d2ab1e8eaea50..bec98efaf1a106d424e4dcc55e0ac6189bcf1d3b 100644
--- a/test/multidomain/embedded/1d3d/1p_1p/main.cc
+++ b/test/multidomain/embedded/1d3d/1p_1p/main.cc
@@ -131,8 +131,8 @@ int main(int argc, char** argv) try
 
     // the mixed dimension type traits
     using Traits = MultiDomainTraits<BulkTypeTag, LowDimTypeTag>;
-    constexpr auto bulkIdx = Traits::template DomainIdx<0>();
-    constexpr auto lowDimIdx = Traits::template DomainIdx<1>();
+    constexpr auto bulkIdx = Traits::template SubDomain<0>::Index();
+    constexpr auto lowDimIdx = Traits::template SubDomain<1>::Index();
 
     // the coupling manager
     using CouplingManager = GetPropType<BulkTypeTag, Properties::CouplingManager>;
diff --git a/test/multidomain/embedded/1d3d/1p_richards/main.cc b/test/multidomain/embedded/1d3d/1p_richards/main.cc
index 06aa6c6fb6502ee3cf0a943b2f5e397034e8f6f6..01dc0c0b8b8d8b9333df8265a84ae4a7960443c7 100644
--- a/test/multidomain/embedded/1d3d/1p_richards/main.cc
+++ b/test/multidomain/embedded/1d3d/1p_richards/main.cc
@@ -124,8 +124,8 @@ int main(int argc, char** argv) try
 
     // the mixed dimension type traits
     using Traits = MultiDomainTraits<BulkTypeTag, LowDimTypeTag>;
-    constexpr auto bulkIdx = Traits::template DomainIdx<0>();
-    constexpr auto lowDimIdx = Traits::template DomainIdx<1>();
+    constexpr auto bulkIdx = Traits::template SubDomain<0>::Index();
+    constexpr auto lowDimIdx = Traits::template SubDomain<1>::Index();
 
     // the coupling manager
     using CouplingManager = GetPropType<BulkTypeTag, Properties::CouplingManager>;
diff --git a/test/multidomain/embedded/2d3d/1p_1p/main.cc b/test/multidomain/embedded/2d3d/1p_1p/main.cc
index 8ab13d94355541386e028aacc32d35db76457496..0fca900c47796e94983e3ba8e2c9f052fdc84abb 100644
--- a/test/multidomain/embedded/2d3d/1p_1p/main.cc
+++ b/test/multidomain/embedded/2d3d/1p_1p/main.cc
@@ -179,8 +179,8 @@ int main(int argc, char** argv) try
 
     // the mixed dimension type traits
     using Traits = MultiDomainTraits<BulkTypeTag, LowDimTypeTag>;
-    constexpr auto bulkIdx = Traits::template DomainIdx<0>();
-    constexpr auto lowDimIdx = Traits::template DomainIdx<1>();
+    constexpr auto bulkIdx = Traits::template SubDomain<0>::Index();
+    constexpr auto lowDimIdx = Traits::template SubDomain<1>::Index();
 
     // the coupling manager
     using CouplingManager = GetPropType<BulkTypeTag, Properties::CouplingManager>;
diff --git a/test/multidomain/facet/1p_1p/analytical/main.cc b/test/multidomain/facet/1p_1p/analytical/main.cc
index 16afe1e2f5779e5a1a6ad6892fa2f95e2383d387..775d9ba5760d7e4eecea8116c56babf516921bd1 100644
--- a/test/multidomain/facet/1p_1p/analytical/main.cc
+++ b/test/multidomain/facet/1p_1p/analytical/main.cc
@@ -246,8 +246,8 @@ int main(int argc, char** argv) try
     using SolutionVector = typename MDTraits::SolutionVector;
     SolutionVector x;
 
-    static const auto bulkId = typename MDTraits::template DomainIdx<0>();
-    static const auto lowDimId = typename MDTraits::template DomainIdx<1>();
+    static const auto bulkId = typename MDTraits::template SubDomain<0>::Index();
+    static const auto lowDimId = typename MDTraits::template SubDomain<1>::Index();
     x[bulkId].resize(bulkFvGridGeometry->numDofs());
     x[lowDimId].resize(lowDimFvGridGeometry->numDofs());
     bulkProblem->applyInitialSolution(x[bulkId]);
diff --git a/test/multidomain/facet/1p_1p/threedomain/main.cc b/test/multidomain/facet/1p_1p/threedomain/main.cc
index b13c76e0ecac785a2d8f32a9a78a21df69664706..f24a0fc57e52bdcbe67bad48783ee8190db1dca1 100644
--- a/test/multidomain/facet/1p_1p/threedomain/main.cc
+++ b/test/multidomain/facet/1p_1p/threedomain/main.cc
@@ -145,12 +145,9 @@ int main(int argc, char** argv) try
     using SolutionVector = typename Traits::SolutionVector;
     SolutionVector x;
 
-    static const auto bulkId = Traits::template DomainIdx<0>();
-    static const auto facetId = Traits::template DomainIdx<1>();
-    static const auto edgeId = Traits::template DomainIdx<2>();
-    x[bulkId].resize(bulkFvGridGeometry->numDofs());
-    x[facetId].resize(facetFvGridGeometry->numDofs());
-    x[edgeId].resize(edgeFvGridGeometry->numDofs());
+    static const auto bulkId = Traits::template SubDomain<0>::Index();
+    static const auto facetId = Traits::template SubDomain<1>::Index();
+    static const auto edgeId = Traits::template SubDomain<2>::Index();
     bulkProblem->applyInitialSolution(x[bulkId]);
     facetProblem->applyInitialSolution(x[facetId]);
     edgeProblem->applyInitialSolution(x[edgeId]);
diff --git a/test/multidomain/facet/tracer_tracer/main.cc b/test/multidomain/facet/tracer_tracer/main.cc
index 45fd9a9cb318c2805848f52003ed6484b77d20a4..357ee5c56a942e739fa44c87b10ee4556d5281fd 100644
--- a/test/multidomain/facet/tracer_tracer/main.cc
+++ b/test/multidomain/facet/tracer_tracer/main.cc
@@ -250,8 +250,8 @@ int main(int argc, char** argv) try
     using OnePMDTraits = typename OnePTestTraits::MDTraits;
     using SolutionVector = typename OnePMDTraits::SolutionVector;
     using JacobianMatrix = typename OnePMDTraits::JacobianMatrix;
-    static const auto bulkId = typename OnePMDTraits::template DomainIdx<0>();
-    static const auto lowDimId = typename OnePMDTraits::template DomainIdx<1>();
+    static const auto bulkId = typename OnePMDTraits::template SubDomain<0>::Index();
+    static const auto lowDimId = typename OnePMDTraits::template SubDomain<1>::Index();
 
     ////////////////////////////////////////////////////////////
     // run stationary, simgle-phase problem on this grid
diff --git a/test/multidomain/poromechanics/el1p/main.cc b/test/multidomain/poromechanics/el1p/main.cc
index 226478ac8dbf37c815f8157ecec452a7057c3227..50175c4a693d63bf75f4be51788130b4b3f56d8a 100644
--- a/test/multidomain/poromechanics/el1p/main.cc
+++ b/test/multidomain/poromechanics/el1p/main.cc
@@ -130,8 +130,8 @@ int main(int argc, char** argv) try
     using SolutionVector = typename Traits::SolutionVector;
     SolutionVector x;
 
-    static const auto onePId = Traits::template DomainIdx<0>();
-    static const auto poroMechId = Traits::template DomainIdx<1>();
+    static const auto onePId = Traits::template SubDomain<0>::Index();
+    static const auto poroMechId = Traits::template SubDomain<1>::Index();
     x[onePId].resize(onePFvGridGeometry->numDofs());
     x[poroMechId].resize(poroMechFvGridGeometry->numDofs());
     onePProblem->applyInitialSolution(x[onePId]);
diff --git a/test/multidomain/poromechanics/el2p/main.cc b/test/multidomain/poromechanics/el2p/main.cc
index 65a170a88a420d6e41bef56b9857dff995658153..f833946d76c4cb85d6acd5ff2718fd7fb5eeb604 100644
--- a/test/multidomain/poromechanics/el2p/main.cc
+++ b/test/multidomain/poromechanics/el2p/main.cc
@@ -132,8 +132,8 @@ int main(int argc, char** argv) try
     using SolutionVector = typename Traits::SolutionVector;
     SolutionVector x;
 
-    static const auto twoPId = Traits::template DomainIdx<0>();
-    static const auto poroMechId = Traits::template DomainIdx<1>();
+    static const auto twoPId = Traits::template SubDomain<0>::Index();
+    static const auto poroMechId = Traits::template SubDomain<1>::Index();
     x[twoPId].resize(twoPFvGridGeometry->numDofs());
     x[poroMechId].resize(poroMechFvGridGeometry->numDofs());
     twoPProblem->applyInitialSolution(x[twoPId]);