diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt
index f7ea8bb601bd42ca26054537cdfdf28f24bcaf70..a0e99b95c8475e766e32b282fbbbf0b92af6504c 100644
--- a/dumux/assembly/CMakeLists.txt
+++ b/dumux/assembly/CMakeLists.txt
@@ -7,6 +7,7 @@ cclocalresidual.hh
diffmethod.hh
fvassembler.hh
fvlocalresidual.hh
+jacobianpattern.hh
staggeredfvassembler.hh
staggeredlocalassembler.hh
staggeredlocalresidual.hh
diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh
index 4310071b2a7d97005dc21e045ac3cd38a330c8a6..abeeb48e597ae5335e52e7f881cee5f92451dfa3 100644
--- a/dumux/assembly/fvassembler.hh
+++ b/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
diff --git a/dumux/assembly/jacobianpattern.hh b/dumux/assembly/jacobianpattern.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1847c116213aafef95713ded2028032301323ec7
--- /dev/null
+++ b/dumux/assembly/jacobianpattern.hh
@@ -0,0 +1,110 @@
+// -*- 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
diff --git a/dumux/discretization/box/fvgridgeometry.hh b/dumux/discretization/box/fvgridgeometry.hh
index e3b457bebea5fa0230ab8c757d40cc9884fdf58f..71ea05f596b0285cb446132659b122a2beb95e32 100644
--- a/dumux/discretization/box/fvgridgeometry.hh
+++ b/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) {}
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index dfba3592e209700e45d4642e397f517b5433f180..248c0129d9739c625ab49c728ab6f79c0414d3da 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/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
diff --git a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
index f71a30ea60989698801ff931a22dad27b616e5a4..9d38e45f97d540192eec3a8dcf194fdddbce4ba2 100644
--- a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
+++ b/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)