From a878e4c819113888e5e827b90e760eed58a0a3f9 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 9 Jan 2018 20:22:26 +0100
Subject: [PATCH] [assembly] Make getJacobianPattern a free function
independent of assembler
---
dumux/assembly/CMakeLists.txt | 1 +
dumux/assembly/fvassembler.hh | 66 +----------
dumux/assembly/jacobianpattern.hh | 110 ++++++++++++++++++
dumux/discretization/box/fvgridgeometry.hh | 7 ++
.../cellcentered/mpfa/fvgridgeometry.hh | 7 ++
.../cellcentered/tpfa/fvgridgeometry.hh | 7 ++
6 files changed, 134 insertions(+), 64 deletions(-)
create mode 100644 dumux/assembly/jacobianpattern.hh
diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt
index f7ea8bb601..a0e99b95c8 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 4310071b2a..abeeb48e59 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 0000000000..1847c11621
--- /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 e3b457bebe..71ea05f596 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 dfba3592e2..248c0129d9 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 f71a30ea60..9d38e45f97 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)
--
GitLab