diff --git a/dumux/linear/matrixconverter.hh b/dumux/linear/matrixconverter.hh
index 5fb0a395ddd201d6342d4efea7478b593c1d214f..ea398feb3e4c000a92e7c22e1f50ad049d9be769 100644
--- a/dumux/linear/matrixconverter.hh
+++ b/dumux/linear/matrixconverter.hh
@@ -25,13 +25,12 @@
 #define DUMUX_MATRIX_CONVERTER
 
 #include <cmath>
+#include <utility>
 #include <dune/common/indices.hh>
 #include <dune/common/hybridutilities.hh>
 #include <dune/istl/bvector.hh>
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixindexset.hh>
-#include <dune/istl/multitypeblockvector.hh>
-#include <dune/istl/multitypeblockmatrix.hh>
 
 #include <dumux/common/parameters.hh>
 
@@ -87,7 +86,7 @@ private:
         occupationPattern.resize(numRows, numRows);
 
         // lambda function to fill the occupation pattern
-        auto addIndices = [&occupationPattern](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol)
+        auto addIndices = [&](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol)
         {
             using std::abs;
             static const Scalar eps = getParam<Scalar>("MatrixConverter.DeletePatternEntriesBelowAbsThreshold", -1.0);
@@ -104,11 +103,12 @@ private:
         };
 
         // fill the pattern
+        using namespace Dune::Hybrid;
         std::size_t rowIndex = 0;
-        Dune::Hybrid::forEach(A, [&addIndices, &rowIndex, numRows](const auto& rowOfMultiTypeMatrix)
+        forEach(std::make_index_sequence<MultiTypeBlockMatrix::N()>(), [&A, &addIndices, &rowIndex, numRows](const auto i)
         {
             std::size_t colIndex = 0;
-            Dune::Hybrid::forEach(rowOfMultiTypeMatrix, [&addIndices, &colIndex, &rowIndex, numRows](const auto& subMatrix)
+            forEach(A[i], [&](const auto& subMatrix)
             {
                 addIndices(subMatrix, rowIndex, colIndex);
 
@@ -137,7 +137,7 @@ private:
         const auto numRows = M.N();
 
         // lambda function to copy the values
-        auto copyValues = [&M](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol)
+        auto copyValues = [&](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol)
         {
             using std::abs;
             static const Scalar eps = getParam<Scalar>("MatrixConverter.DeletePatternEntriesBelowAbsThreshold", -1.0);
@@ -154,11 +154,12 @@ private:
 
         };
 
+        using namespace Dune::Hybrid;
         std::size_t rowIndex = 0;
-        Dune::Hybrid::forEach(A, [&copyValues, &rowIndex, numRows](const auto& rowOfMultiTypeMatrix)
+        forEach(std::make_index_sequence<MultiTypeBlockMatrix::N()>(), [&A, &copyValues, &rowIndex, numRows](const auto i)
         {
             std::size_t colIndex = 0;
-            Dune::Hybrid::forEach(rowOfMultiTypeMatrix, [&copyValues, &colIndex, &rowIndex, numRows](const auto& subMatrix)
+            forEach(A[i], [&](const auto& subMatrix)
             {
                 copyValues(subMatrix, rowIndex, colIndex);
 
@@ -217,7 +218,7 @@ public:
         bTmp.resize(b.dim());
 
         std::size_t startIndex = 0;
-        Dune::Hybrid::forEach(b, [&bTmp, &startIndex](const auto& subVector)
+        Dune::Hybrid::forEach(b, [&](const auto& subVector)
         {
             const auto numEq = std::decay_t<decltype(subVector)>::block_type::size();
 
@@ -240,7 +241,7 @@ public:
     static void retrieveValues(MultiTypeBlockVector& x, const BlockVector& y)
     {
         std::size_t startIndex = 0;
-        Dune::Hybrid::forEach(x, [&y, &startIndex](auto& subVector)
+        Dune::Hybrid::forEach(x, [&](auto& subVector)
         {
             const auto numEq = std::decay_t<decltype(subVector)>::block_type::size();
 
diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh
index 5040a1e15e6bf99201100f2552dbd8b97a7b3f37..401ad03d74cccb75bfc352af1e9cc7d8e5ef6d21 100644
--- a/dumux/multidomain/fvassembler.hh
+++ b/dumux/multidomain/fvassembler.hh
@@ -188,7 +188,7 @@ public:
         resetResidual_();
 
         using namespace Dune::Hybrid;
-        forEach(integralRange(Dune::Hybrid::size(*jacobian_)), [&](const auto domainId)
+        forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainId)
         {
             auto& jacRow = (*jacobian_)[domainId];
             auto& subRes = (*residual_)[domainId];
@@ -269,9 +269,9 @@ public:
     void setJacobianBuildMode(JacobianMatrix& jac) const
     {
         using namespace Dune::Hybrid;
-        forEach(jac, [](auto& jacRow)
+        forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto i)
         {
-            forEach(jacRow, [](auto& jacBlock)
+            forEach(jac[i], [&](auto& jacBlock)
             {
                 using BlockType = std::decay_t<decltype(jacBlock)>;
                 if (jacBlock.buildMode() == BlockType::BuildMode::unknown)
@@ -288,7 +288,7 @@ public:
     void setJacobianPattern(JacobianMatrix& jac) const
     {
         using namespace Dune::Hybrid;
-        forEach(integralRange(Dune::Hybrid::size(jac)), [&](const auto domainI)
+        forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainI)
         {
             forEach(integralRange(Dune::Hybrid::size(jac[domainI])), [&](const auto domainJ)
             {