diff --git a/dumux/common/math.hh b/dumux/common/math.hh
index 8222a2ca9d92d9dfd9c333aae517891c0673931c..cc5eb1d58ac87311357f16dc56cd35e37b8b9eb6 100644
--- a/dumux/common/math.hh
+++ b/dumux/common/math.hh
@@ -553,6 +553,26 @@ Dune::FieldMatrix<Scalar, n, m> getTransposed(const Dune::FieldMatrix<Scalar, m,
     return T;
 }
 
+/*!
+ * \brief Transpose a DynamicMatrix
+ *
+ * \param M The matrix to be transposed
+ */
+template <class Scalar>
+Dune::DynamicMatrix<Scalar> getTransposed(const Dune::DynamicMatrix<Scalar>& M)
+{
+    std::size_t rows_T = M.M();
+    std::size_t cols_T = M.N();
+
+    Dune::DynamicMatrix<Scalar> M_T(rows_T, cols_T, 0.0);
+
+    for (std::size_t i = 0; i < rows_T; ++i)
+        for (std::size_t j = 0; j < cols_T; ++j)
+            M_T[i][j] = M[j][i];
+
+    return M_T;
+}
+
 /*!
  * \brief Multiply two dynamic matrices
  *
@@ -577,6 +597,27 @@ Dune::DynamicMatrix<Scalar> multiplyMatrices(const Dune::DynamicMatrix<Scalar> &
     return result;
 }
 
+/*!
+ * \brief Trace of dynamic matrix
+ *
+ * \param M The dynamic matrix
+ */
+template <class Scalar>
+Scalar trace(const Dune::DynamicMatrix<Scalar>& M)
+{
+    std::size_t rows_T = M.M();
+    std::size_t cols_T = M.N();
+
+    DUNE_ASSERT_BOUNDS(rows_T == cols_T);
+
+    Scalar trace = 0.0;
+
+    for (std::size_t i = 0; i < rows_T; ++i)
+        trace += M[i][i];
+
+    return trace;
+}
+
 } // end namespace Dumux
 
 #endif
diff --git a/dumux/freeflow/staggered/indices.hh b/dumux/freeflow/staggered/indices.hh
index 9ff3adb4c10c970915da31a37fd9fc41c3a350d4..704590dac4c5db4a160bf7a1ea16c0bd31838cb7 100644
--- a/dumux/freeflow/staggered/indices.hh
+++ b/dumux/freeflow/staggered/indices.hh
@@ -50,7 +50,7 @@ struct NavierStokesCommonIndices
     static constexpr int momentumXBalanceIdx = momentumBalanceIdx; //!< Index of the momentum balance equation
     static constexpr int momentumYBalanceIdx = momentumBalanceIdx + 1; //!< Index of the momentum balance equation
     static constexpr int momentumZBalanceIdx = momentumBalanceIdx + 2; //!< Index of the momentum balance equation
-    static constexpr int velocityIdx = PVOffset + momentumBalanceIdx; //!< Index of the velocity in a solution vector (NOTE: This PV lives in a different vector than the pressure)
+    static constexpr int velocityIdx = momentumBalanceIdx; //!< Index of the velocity in a solution vector (NOTE: This PV lives in a different vector than the pressure)
     static constexpr int velocityXIdx = velocityIdx; //!< Index of the velocity in a solution vector (NOTE: This PV lives in a different vector than the pressure)
     static constexpr int velocityYIdx = velocityIdx + 1; //!< Index of the velocity in a solution vector (NOTE: This PV lives in a different vector than the pressure)
     static constexpr int velocityZIdx = velocityIdx + 2; //!< Index of the velocity in a solution vector (NOTE: This PV lives in a different vector than the pressure)
diff --git a/dumux/implicit/staggered/localresidual.hh b/dumux/implicit/staggered/localresidual.hh
index a2d11b36c628f627be5fc3d58eb9164ae5864aa1..9353b71b16ca81674d2f584cdec457fee1221c66 100644
--- a/dumux/implicit/staggered/localresidual.hh
+++ b/dumux/implicit/staggered/localresidual.hh
@@ -313,7 +313,7 @@ protected:
                             const ElementFluxVariablesCache& elemFluxVarsCache)
     {
         if(!scvf.boundary())
-            faceResiduals_[scvf.localFaceIdx()] += asImp_().computeFluxForFace(element, scvf, fvGeometry, elemVolVars, globalFaceVars);
+            faceResiduals_[scvf.localFaceIdx()] += asImp_().computeFluxForFace(element, scvf, fvGeometry, elemVolVars, globalFaceVars, elemFluxVarsCache);
     }
 
      /*!
diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/implicit/staggered/propertydefaults.hh
index 3917d92c883135b3ed37806e885dd6c21d9f310b..5e5f215a59dba5b9dec31177d9d94d049268fe5e 100644
--- a/dumux/implicit/staggered/propertydefaults.hh
+++ b/dumux/implicit/staggered/propertydefaults.hh
@@ -49,6 +49,8 @@
 
 #include <dumux/linear/seqsolverbackend.hh>
 
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+
 #include "assembler.hh"
 #include "localresidual.hh"
 #include "localjacobian.hh"
@@ -246,6 +248,7 @@ SET_BOOL_PROP(StaggeredModel, VtkWriteFaceData, true);
 //! For compatibility
 SET_BOOL_PROP(StaggeredModel, EnableInteriorBoundaries, false);
 
+SET_TYPE_PROP(StaggeredModel, VtkOutputModule, StaggeredVtkOutputModule<TypeTag>);
 
 } // namespace Properties