[stokes_solver] Fix mass matrix for linear CVFE schemes

dumux/linear/stokes_solver.hh

@@ -35,6 +35,7 @@
 #include <dumux/linear/parallelmatrixadapter.hh>
 #include <dumux/discretization/extrusion.hh>
 #include <dumux/linear/symmetrize_constraints.hh>
+#include <dumux/assembly/jacobianpattern.hh>
 
 namespace Dumux::Detail {
 
@@ -400,36 +401,76 @@ private:
         return std::make_shared<LinearOperator>(massMatrix);
     }
 
-    template<class M>
+    template<class M, class DM = typename PressureGG::DiscretizationMethod>
     std::shared_ptr<M> createMassMatrix_()
     {
         auto massMatrix = std::make_shared<M>();
         massMatrix->setBuildMode(M::random);
         const auto numDofs = pGridGeometry_->numDofs();
 
-        Dune::MatrixIndexSet pattern;
-        pattern.resize(numDofs, numDofs);
-        for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
-            pattern.add(globalI, globalI);
-        pattern.exportIdx(*massMatrix);
+        if constexpr (DM{} == DiscretizationMethods::cctpfa || DM{} == DiscretizationMethods::ccmpfa)
+        {
+            Dune::MatrixIndexSet pattern;
+            pattern.resize(numDofs, numDofs);
+            for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+                pattern.add(globalI, globalI);
+            pattern.exportIdx(*massMatrix);
+
+            const auto& gv = pGridGeometry_->gridView();
+            auto fvGeometry = localView(*pGridGeometry_);
+            for (const auto& element : elements(gv))
+            {
+                fvGeometry.bindElement(element);
+                for (const auto& scv : scvs(fvGeometry))
+                {
+                    using Extrusion = Extrusion_t<PressureGG>;
+                    const auto dofIndex = scv.dofIndex();
+                    if (element.partitionType() == Dune::GhostEntity) // do not modify ghosts
+                        (*massMatrix)[dofIndex][dofIndex] = 1.0;
+                    else
+                        (*massMatrix)[dofIndex][dofIndex] += weight_*Extrusion::volume(fvGeometry, scv)/(2.0*viscosity_);
+                }
+            }
 
-        const auto& gv = pGridGeometry_->gridView();
-        auto fvGeometry = localView(*pGridGeometry_);
-        for (const auto& element : elements(gv))
+            return massMatrix;
+        }
+        else if constexpr (DM{} == DiscretizationMethods::box || DM{} == DiscretizationMethods::fcdiamond)
         {
-            fvGeometry.bindElement(element);
-            for (const auto& scv : scvs(fvGeometry))
+            getJacobianPattern<true>(*pGridGeometry_).exportIdx(*massMatrix);
+            (*massMatrix) = 0.0;
+
+            const auto& gv = pGridGeometry_->gridView();
+            auto fvGeometry = localView(*pGridGeometry_);
+            std::vector<Dune::FieldVector<double, 1>> values;
+
+            for (const auto& element : elements(gv))
             {
-                using Extrusion = Extrusion_t<PressureGG>;
-                const auto dofIndex = scv.dofIndex();
-                if (element.partitionType() == Dune::GhostEntity) // do not modify ghosts
-                    (*massMatrix)[dofIndex][dofIndex] = 1.0;
-                else
-                    (*massMatrix)[dofIndex][dofIndex] += weight_*Extrusion::volume(fvGeometry, scv)/(2.0*viscosity_);
+                const auto geometry = element.geometry();
+                const auto& localBasis = pGridGeometry_->feCache().get(element.type()).localBasis();
+                const auto numLocalDofs = localBasis.size();
+                values.resize(numLocalDofs);
+
+                fvGeometry.bindElement(element);
+                for (const auto& scvJ : scvs(fvGeometry))
+                {
+                    // Use mid-point rule (only works for linear ansatz functions)
+                    const auto globalJ = scvJ.dofIndex();
+                    const auto qWeightJ = PressureGG::Extrusion::volume(fvGeometry, scvJ);
+                    const auto quadPos = geometry.local(scvJ.center());
+                    localBasis.evaluateFunction(quadPos, values);
+
+                    for (const auto& scvI : scvs(fvGeometry))
+                    {
+                        const auto valueIJ = values[scvI.localDofIndex()]*qWeightJ/(2.0*viscosity_);
+                        (*massMatrix)[scvI.dofIndex()][globalJ][0][0] += valueIJ;
+                    }
+                }
             }
+
+            return massMatrix;
         }
 
-        return massMatrix;
+        DUNE_THROW(Dune::NotImplemented, "Mass matrix for discretization method not implemented");
     }
 
     double density_, viscosity_, weight_;