[staggered][model] Slight modifications of staggered model

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

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)

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);
     }
 
      /*!

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