diff --git a/dune/gfe/assemblers/mixedgfeassembler.hh b/dune/gfe/assemblers/mixedgfeassembler.hh
index 21203f0d04c279280297f651170778fba9d2f45f..ebfe80f2f2da579d4de8239d6351c6f1c9c5d70b 100644
--- a/dune/gfe/assemblers/mixedgfeassembler.hh
+++ b/dune/gfe/assemblers/mixedgfeassembler.hh
@@ -227,16 +227,16 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
         auto col = localView.index(localIndexCol);
 
         if (row[0]==0 and col[0]==0)
-          hessian[_0][_0][row[1]][col[1]] += localStiffness_->A00_[i][j];
+          hessian[_0][_0][row[1]][col[1]] += localStiffness_->A_[_0][_0][i][j];
 
         if (row[0]==0 and col[0]==1)
-          hessian[_0][_1][row[1]][col[1]] += localStiffness_->A01_[i][j-nDofs0];
+          hessian[_0][_1][row[1]][col[1]] += localStiffness_->A_[_0][_1][i][j-nDofs0];
 
         if (row[0]==1 and col[0]==0)
-          hessian[_1][_0][row[1]][col[1]] += localStiffness_->A10_[i-nDofs0][j];
+          hessian[_1][_0][row[1]][col[1]] += localStiffness_->A_[_1][_0][i-nDofs0][j];
 
         if (row[0]==1 and col[0]==1)
-          hessian[_1][_1][row[1]][col[1]] += localStiffness_->A11_[i-nDofs0][j-nDofs0];
+          hessian[_1][_1][row[1]][col[1]] += localStiffness_->A_[_1][_1][i-nDofs0][j-nDofs0];
       }
 
       // Add local gradient to global gradient
diff --git a/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh b/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
index dab88c545df7d68b7d61a6d816ac450b0882b62d..b74d94fb58fac620e27fbf04e2aa07eb150bb3fd 100644
--- a/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
+++ b/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
@@ -2,7 +2,9 @@
 #define DUNE_GFE_MIXEDLOCALGEODESICFESTIFFNESS_HH
 
 #include <dune/common/fmatrix.hh>
+
 #include <dune/istl/matrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
 
 
 /** \brief Abstract base class for second-order energy approximations on one grid element
@@ -21,8 +23,8 @@ class MixedLocalGeodesicFEStiffness
 public:
 
   //! Dimension of a tangent space
-  constexpr static int deformationBlocksize = DeformationTargetSpace::TangentVector::dimension;
-  constexpr static int orientationBlocksize = OrientationTargetSpace::TangentVector::dimension;
+  constexpr static int blocksize0 = DeformationTargetSpace::TangentVector::dimension;
+  constexpr static int blocksize1 = OrientationTargetSpace::TangentVector::dimension;
 
   /** \brief Assemble the local stiffness matrix at the current position
    */
@@ -40,12 +42,13 @@ public:
                      const std::vector<DeformationTargetSpace>& localDeformationConfiguration,
                      const std::vector<OrientationTargetSpace>& localOrientationConfiguration) const = 0;
 
-  // assembled data
-  Dune::Matrix<Dune::FieldMatrix<RT, deformationBlocksize, deformationBlocksize> > A00_;
-  Dune::Matrix<Dune::FieldMatrix<RT, deformationBlocksize, orientationBlocksize> > A01_;
-  Dune::Matrix<Dune::FieldMatrix<RT, orientationBlocksize, deformationBlocksize> > A10_;
-  Dune::Matrix<Dune::FieldMatrix<RT, orientationBlocksize, orientationBlocksize> > A11_;
+  // assembled tangent matrix
+  using Row0 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<RT, blocksize0, blocksize0> >,
+      Dune::Matrix<Dune::FieldMatrix<RT, blocksize0, blocksize1> > >;
+  using Row1 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<RT, blocksize1, blocksize0> >,
+      Dune::Matrix<Dune::FieldMatrix<RT, blocksize1, blocksize1> > >;
 
+  Dune::MultiTypeBlockMatrix<Row0, Row1> A_;
 };
 
 #endif
diff --git a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
index 7fb2f4c3498a99fac4d6d1148b0477d6d8a41f8c..603cb8119a93b5511aa7e334051fd8761793117f 100644
--- a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
@@ -351,7 +351,9 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
   // For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look
   // at the Riemannian Hessian".
 
-  this->A00_.setSize(nDofs0,nDofs0);
+  using namespace Dune::Indices;
+
+  this->A_[_0][_0].setSize(nDofs0,nDofs0);
 
   for (size_t col=0; col<nDofs0; col++) {
 
@@ -365,14 +367,14 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
         embeddedHessian00[row][col].mv(z,semiEmbeddedProduct);
 
         for (int subRow=0; subRow<blocksize0; subRow++)
-          this->A00_[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+          this->A_[_0][_0][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
       }
 
     }
 
   }
 
-  this->A01_.setSize(nDofs0,nDofs1);
+  this->A_[_0][_1].setSize(nDofs0,nDofs1);
 
   for (size_t col=0; col<nDofs1; col++) {
 
@@ -386,14 +388,14 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
         embeddedHessian01[row][col].mv(z,semiEmbeddedProduct);
 
         for (int subRow=0; subRow<blocksize0; subRow++)
-          this->A01_[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+          this->A_[_0][_1][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
       }
 
     }
 
   }
 
-  this->A10_.setSize(nDofs1,nDofs0);
+  this->A_[_1][_0].setSize(nDofs1,nDofs0);
 
   for (size_t col=0; col<nDofs0; col++) {
 
@@ -407,14 +409,14 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
         embeddedHessian10[row][col].mv(z,semiEmbeddedProduct);
 
         for (int subRow=0; subRow<blocksize1; subRow++)
-          this->A10_[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+          this->A_[_1][_0][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
       }
 
     }
 
   }
 
-  this->A11_.setSize(nDofs1,nDofs1);
+  this->A_[_1][_1].setSize(nDofs1,nDofs1);
 
   for (size_t col=0; col<nDofs1; col++) {
 
@@ -428,7 +430,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
         embeddedHessian11[row][col].mv(z,semiEmbeddedProduct);
 
         for (int subRow=0; subRow<blocksize1; subRow++)
-          this->A11_[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+          this->A_[_1][_1][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
       }
 
     }
@@ -460,7 +462,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
       typename TargetSpace0::TangentVector tmp2;
       orthonormalFrame0[row].mv(tmp1,tmp2);
 
-      this->A00_[row][row][subRow] += tmp2;
+      this->A_[_0][_0][row][row][subRow] += tmp2;
     }
 
   }
@@ -475,7 +477,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
       typename TargetSpace1::TangentVector tmp2;
       orthonormalFrame1[row].mv(tmp1,tmp2);
 
-      this->A11_[row][row][subRow] += tmp2;
+      this->A_[_1][_1][row][row][subRow] += tmp2;
     }
 
   }