[assembly] Make getJacobianPattern a free function independent of assembler

dumux/assembly/CMakeLists.txt

@@ -7,6 +7,7 @@ cclocalresidual.hh
 diffmethod.hh
 fvassembler.hh
 fvlocalresidual.hh
+jacobianpattern.hh
 staggeredfvassembler.hh
 staggeredlocalassembler.hh
 staggeredlocalresidual.hh

dumux/assembly/fvassembler.hh

@@ -33,6 +33,7 @@
 #include <dumux/discretization/methods.hh>
 #include <dumux/parallel/vertexhandles.hh>
 
+#include "jacobianpattern.hh"
 #include "diffmethod.hh"
 #include "boxlocalassembler.hh"
 #include "cclocalassembler.hh"
@@ -227,11 +228,7 @@ public:
         jacobian_->setSize(numDofs, numDofs);
 
         // create occupation pattern of the jacobian
-        Dune::MatrixIndexSet occupationPattern;
-        occupationPattern.resize(numDofs, numDofs);
-
-        // set the jacobian pattern depending on space and time discretization
-        setJacobianPattern_(occupationPattern, numDofs);
+        const auto occupationPattern = getJacobianPattern<isImplicit>(fvGridGeometry());
 
         // export pattern to jacobian
         occupationPattern.exportIdx(*jacobian_);
@@ -336,65 +333,6 @@ private:
             DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
     }
 
-    //! Implicit jacobian pattern for cell-centered fv schemes
-    template<class T = TypeTag, typename std::enable_if_t<GET_PROP_VALUE(T, DiscretizationMethod) != DiscretizationMethods::Box, int> = 0>
-    void setJacobianPattern_(Dune::MatrixIndexSet& pattern, std::size_t numDofs)
-    {
-        // matrix pattern for implicit Jacobians
-        if (isImplicit)
-        {
-            for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
-            {
-                pattern.add(globalI, globalI);
-                for (const auto& dataJ : fvGridGeometry().connectivityMap()[globalI])
-                    pattern.add(dataJ.globalJ, globalI);
-
-                // reserve index for additional user defined DOF dependencies
-                // const auto& additionalDofDependencies = problem().getAdditionalDofDependencies(globalI);
-                // for (auto globalJ : additionalDofDependencies)
-                //     pattern.add(globalI, globalJ);
-            }
-        }
-
-        // matrix pattern for explicit Jacobians -> diagonal matrix
-        else
-        {
-            for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
-                pattern.add(globalI, globalI);
-        }
-    }
-
-    //! Implicit jacobian pattern for vertex-centered fv schemes
-    template<class T = TypeTag, typename std::enable_if_t<GET_PROP_VALUE(T, DiscretizationMethod) == DiscretizationMethods::Box, int> = 0>
-    void setJacobianPattern_(Dune::MatrixIndexSet& pattern, std::size_t numDofs)
-    {
-        // matrix pattern for implicit Jacobians
-        if (isImplicit)
-        {
-            static constexpr int dim = GridView::dimension;
-            for (const auto& element : elements(fvGridGeometry().gridView()))
-            {
-                for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
-                {
-                    const auto globalI = fvGridGeometry().vertexMapper().subIndex(element, vIdx, dim);
-                    for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(dim); ++vIdx2)
-                    {
-                        const auto globalJ = fvGridGeometry().vertexMapper().subIndex(element, vIdx2, dim);
-                        pattern.add(globalI, globalJ);
-                        pattern.add(globalJ, globalI);
-                    }
-                }
-            }
-        }
-
-        // matrix pattern for explicit Jacobians -> diagonal matrix
-        else
-        {
-            for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
-                pattern.add(globalI, globalI);
-        }
-    }
-
     /*!
      * \brief A method assembling something per element
      * \note Handles exceptions for parallel runs

dumux/assembly/jacobianpattern.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief Helper function to generate Jacobian pattern for different discretization methods
+ */
+#ifndef DUMUX_JACOBIAN_PATTERN_HH
+#define DUMUX_JACOBIAN_PATTERN_HH
+
+#include <type_traits>
+#include <dune/istl/matrixindexset.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Assembly
+ * \brief Helper function to generate Jacobian pattern for the box method
+ */
+template<bool isImplicit, class GridGeometry,
+         typename std::enable_if_t<(GridGeometry::discretizationMethod == DiscretizationMethods::Box), int> = 0>
+Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
+{
+    const auto numDofs = gridGeometry.numDofs();
+    Dune::MatrixIndexSet pattern;
+    pattern.resize(numDofs, numDofs);
+
+    // matrix pattern for implicit Jacobians
+    if (isImplicit)
+    {
+        static constexpr int dim = std::decay_t<decltype(gridGeometry.gridView())>::dimension;
+        for (const auto& element : elements(gridGeometry.gridView()))
+        {
+            for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
+            {
+                const auto globalI = gridGeometry.vertexMapper().subIndex(element, vIdx, dim);
+                for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(dim); ++vIdx2)
+                {
+                    const auto globalJ = gridGeometry.vertexMapper().subIndex(element, vIdx2, dim);
+                    pattern.add(globalI, globalJ);
+                    pattern.add(globalJ, globalI);
+                }
+            }
+        }
+    }
+
+    // matrix pattern for explicit Jacobians -> diagonal matrix
+    else
+    {
+        for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+            pattern.add(globalI, globalI);
+    }
+
+    return pattern;
+}
+
+/*!
+ * \ingroup Assembly
+ * \brief Helper function to generate Jacobian pattern for cell-centered methods
+ */
+template<bool isImplicit, class GridGeometry,
+         typename std::enable_if_t<( (GridGeometry::discretizationMethod == DiscretizationMethods::CCTpfa)
+                                     || (GridGeometry::discretizationMethod == DiscretizationMethods::CCMpfa) ), int> = 0>
+Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
+{
+    const auto numDofs = gridGeometry.numDofs();
+    Dune::MatrixIndexSet pattern;
+    pattern.resize(numDofs, numDofs);
+
+    // matrix pattern for implicit Jacobians
+    if (isImplicit)
+    {
+        for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+        {
+            pattern.add(globalI, globalI);
+            for (const auto& dataJ : gridGeometry.connectivityMap()[globalI])
+                pattern.add(dataJ.globalJ, globalI);
+        }
+    }
+
+    // matrix pattern for explicit Jacobians -> diagonal matrix
+    else
+    {
+        for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+            pattern.add(globalI, globalI);
+    }
+
+    return pattern;
+}
+
+} // namespace Dumux
+
+#endif

dumux/discretization/box/fvgridgeometry.hh

@@ -28,6 +28,7 @@
 #include <dune/geometry/referenceelements.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
+#include <dumux/discretization/methods.hh>
 #include <dumux/discretization/basefvgridgeometry.hh>
 #include <dumux/discretization/box/boxgeometryhelper.hh>
 #include <dumux/discretization/box/fvelementgeometry.hh>
@@ -71,6 +72,9 @@ class BoxFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
     using GeometryHelper = BoxGeometryHelper<GridView, dim, SubControlVolume, SubControlVolumeFace>;
 
 public:
+    //! export discretization method
+    static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::Box;
+
     //! Constructor
     BoxFVGridGeometry(const GridView gridView)
     : ParentType(gridView) {}
@@ -263,6 +267,9 @@ class BoxFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
     using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 
 public:
+    //! export discretization method
+    static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::Box;
+
     //! Constructor
     BoxFVGridGeometry(const GridView gridView)
     : ParentType(gridView) {}

dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh

@@ -32,6 +32,7 @@
 #include <dumux/common/properties.hh>
 #include <dumux/common/parameters.hh>
 
+#include <dumux/discretization/methods.hh>
 #include <dumux/discretization/basefvgridgeometry.hh>
 #include <dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh>
 #include <dumux/discretization/cellcentered/mpfa/connectivitymap.hh>
@@ -87,6 +88,9 @@ class CCMpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
     using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 
 public:
+    //! export discretization method
+    static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::CCMpfa;
+
     using SecondaryIvIndicatorType = std::function<bool(const Element&, const Intersection&, bool)>;
 
     //! Constructor without indicator function for secondary interaction volumes
@@ -429,6 +433,9 @@ class CCMpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
     using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 
 public:
+    //! export discretization method
+    static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::CCMpfa;
+
     using SecondaryIvIndicator = std::function<bool(const Element&, const Intersection&, bool)>;
 
     //! Constructor without indicator function for secondary interaction volumes

dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh

@@ -28,6 +28,7 @@
 
 #include <dune/common/version.hh>
 
+#include <dumux/discretization/methods.hh>
 #include <dumux/discretization/basefvgridgeometry.hh>
 #include <dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh>
 #include <dumux/discretization/cellcentered/connectivitymap.hh>
@@ -74,6 +75,9 @@ class CCTpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
     friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
 
 public:
+    //! export discretization method
+    static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::CCTpfa;
+
     //! Constructor
     CCTpfaFVGridGeometry(const GridView& gridView)
     : ParentType(gridView)
@@ -331,6 +335,9 @@ class CCTpfaFVGridGeometry<TypeTag, false>  : public BaseFVGridGeometry<TypeTag>
     using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
 
 public:
+    //! export discretization method
+    static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::CCTpfa;
+
     //! Constructor
     CCTpfaFVGridGeometry(const GridView& gridView)
     : ParentType(gridView)