diff --git a/dumux/discretization/cellcentered/mpfa/CMakeLists.txt b/dumux/discretization/cellcentered/mpfa/CMakeLists.txt
index 07756ca4fd6ed7951b0f48a2645d219e78f1b9c0..5e4c4794d76ea53cbc355227c7fe81eb15aaa976 100644
--- a/dumux/discretization/cellcentered/mpfa/CMakeLists.txt
+++ b/dumux/discretization/cellcentered/mpfa/CMakeLists.txt
@@ -17,6 +17,7 @@ helper.hh
interactionvolumebase.hh
interactionvolumedatahandle.hh
localassemblerbase.hh
+localassemblerhelper.hh
localfacedata.hh
methods.hh
scvgradients.hh
diff --git a/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh
index 381466e123dcd109dad5d5374b7e159e580d4b08..b8ce045537b31d3a654e7fc1ce38050b51722abf 100644
--- a/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh
+++ b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh
@@ -26,23 +26,21 @@
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_BASE_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_BASE_HH
-#include <utility>
-#include <type_traits>
-
#include <dune/common/exceptions.hh>
#include <dune/common/reservedvector.hh>
#include <dumux/common/typetraits/isvalid.hh>
+#include "localassemblerhelper.hh"
+
+namespace Dumux {
-namespace Dumux
-{
/*!
* \ingroup CCMpfaDiscretization
* \brief Defines the general interface of the local assembler
* classes for the assembly of the interaction volume-local
* transmissibility matrix. Specializations have to be provided
* for the available interaction volume implementations. these
- * should derive from this base clases.
+ * should derive from this base class.
*
* \tparam P The problem type
* \tparam EG The element finite volume geometry
@@ -50,34 +48,13 @@ namespace Dumux
*/
template< class P, class EG, class EV >
class InteractionVolumeAssemblerBase
+: public InteractionVolumeAssemblerHelper
{
using Problem = P;
using FVElementGeometry = EG;
using ElementVolumeVariables = EV;
- // Helper structs to detect if matrix has resize function
- struct HasMatrixResize
- {
- template<class M>
- auto operator()(const M& m) -> decltype(std::declval<M>().resize(0, 0))
- {}
- };
-
- // Helper structs to detect if vector has resize function
- struct HasVectorResize
- {
- template<class V>
- auto operator()(const V& v) -> decltype(std::declval<V>().resize(0))
- {}
- };
-
- template<class Matrix>
- static constexpr auto matrixHasResizeFunction()
- { return decltype( isValid(HasMatrixResize())(std::declval<Matrix>()) )::value; }
-
- template<class Vector>
- static constexpr auto vectorHasResizeFunction()
- { return decltype( isValid(HasVectorResize())(std::declval<Vector>()) )::value; }
+ using Helper = InteractionVolumeAssemblerHelper;
public:
/*!
@@ -138,72 +115,6 @@ class InteractionVolumeAssemblerBase
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function for the cell/Dirichlet unknowns");
}
- /*!
- * \brief Assembles the vector of face unknowns within an interaction volume.
- * \note This requires the data handle to be fully assembled already.
- *
- * \param handle The data handle in which the vector is stored
- * \param iv The interaction volume
- */
- template< class DataHandle, class IV >
- static typename IV::Traits::MatVecTraits::FaceVector
- assembleFaceUnkowns(const DataHandle& handle, const IV& iv)
- {
- typename IV::Traits::MatVecTraits::FaceVector u;
- resizeVector_(u, iv.numFaces());
-
- handle.AB().mv(handle.uj(), u);
-
- // maybe add gravity terms
- if (handle.deltaG().size() == iv.numUnknowns())
- handle.AB().umv(handle.deltaG(), u);
-
- return u;
- }
-
- /*!
- * \brief Assembles the solution gradients in the
- * sub-control volumes within an interaction volume.
- * \note This requires the data handle to be fully assembled already.
- *
- * \param handle The data handle in which the vector is stored
- * \param iv The interaction volume
- */
- template< class DataHandle, class IV >
- static std::vector< typename IV::Traits::LocalScvType::GlobalCoordinate >
- assembleScvGradients(const DataHandle& handle, const IV& iv)
- {
- const auto u = assembleFaceUnkowns(handle, iv);
-
- using LocalScv = typename IV::Traits::LocalScvType;
- using Gradient = typename LocalScv::GlobalCoordinate;
-
- std::vector<Gradient> result; result.reserve(iv.numScvs());
- for (unsigned int scvIdx = 0; scvIdx < iv.numScvs(); ++scvIdx)
- {
- const auto& scv = iv.localScv(scvIdx);
-
- Gradient gradU(0.0);
- for (unsigned int dir = 0; dir < LocalScv::myDimension; ++dir)
- {
- auto nu = scv.nu(dir);
-
- // obtain face pressure
- const auto& scvf = iv.localScvf( scv.localScvfIndex(dir) );
- const auto faceU = !scvf.isDirichlet() ? u[scvf.localDofIndex()]
- : handle.uj()[scvf.localDofIndex()];
-
- nu *= faceU - handle.uj()[scv.localDofIndex()];
- gradU += nu;
- }
-
- gradU /= scv.detX();
- result.emplace_back( std::move(gradU) );
- }
-
- return result;
- }
-
/*!
* \brief Assembles the gravitational flux contributions on the scvfs within an
* interaction volume.
@@ -224,10 +135,10 @@ class InteractionVolumeAssemblerBase
auto& g = handle.g();
auto& deltaG = handle.deltaG();
auto& outsideG = handle.gOutside();
- resizeVector_(g, iv.numFaces());
- resizeVector_(deltaG, iv.numUnknowns());
+ Helper::resizeVector(g, iv.numFaces());
+ Helper::resizeVector(deltaG, iv.numUnknowns());
if (isSurfaceGrid)
- resizeVector_(outsideG, iv.numFaces());
+ Helper::resizeVector(outsideG, iv.numFaces());
//! For each face, we...
//! - arithmetically average the phase densities
@@ -260,7 +171,7 @@ class InteractionVolumeAssemblerBase
Scalar rho;
if (isSurfaceGrid)
- resizeVector_(outsideG[faceIdx], numOutsideFaces);
+ Helper::resizeVector(outsideG[faceIdx], numOutsideFaces);
if (!curLocalScvf.isDirichlet())
{
@@ -316,7 +227,7 @@ class InteractionVolumeAssemblerBase
{
using FaceVector = typename IV::Traits::MatVecTraits::FaceVector;
FaceVector AG;
- resizeVector_(AG, iv.numUnknowns());
+ Helper::resizeVector(AG, iv.numUnknowns());
handle.A().mv(deltaG, AG);
// compute gravitational accelerations
@@ -344,39 +255,6 @@ class InteractionVolumeAssemblerBase
}
}
-protected:
- //! resizes a matrix to the given sizes (specialization for dynamic matrix type)
- template< class Matrix,
- class size_type,
- std::enable_if_t<matrixHasResizeFunction<Matrix>(), int> = 0 >
- static void resizeMatrix_(Matrix& M, size_type rows, size_type cols)
- {
- M.resize(rows, cols);
- }
-
- //! resizes a matrix to the given sizes (specialization for static matrix type - do nothing)
- template< class Matrix,
- class size_type,
- std::enable_if_t<!matrixHasResizeFunction<Matrix>(), int> = 0 >
- static void resizeMatrix_(Matrix& M, size_type rows, size_type cols)
- {}
-
- //! resizes a vector to the given size (specialization for dynamic matrix type)
- template< class Vector,
- class size_type,
- std::enable_if_t<vectorHasResizeFunction<Vector>(), int> = 0 >
- static void resizeVector_(Vector& v, size_type size)
- {
- v.resize(size);
- }
-
- //! resizes a vector to the given size (specialization for static vector type - do nothing)
- template< class Vector,
- class size_type,
- std::enable_if_t<!vectorHasResizeFunction<Vector>(), int> = 0 >
- static void resizeVector_(Vector& v, size_type rows)
- {}
-
private:
// pointers to the data required for assembly
const Problem* problemPtr_;
diff --git a/dumux/discretization/cellcentered/mpfa/localassemblerhelper.hh b/dumux/discretization/cellcentered/mpfa/localassemblerhelper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0ce355049ea036114188f1cb15d9b4f370cdbfdb
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/localassemblerhelper.hh
@@ -0,0 +1,171 @@
+// -*- 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 CCMpfaDiscretization
+ * \brief A class that contains helper functions as well as functionality
+ * which is common to different mpfa schemes and which solely
+ * operate on the interaction volume.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_HELPER_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_HELPER_HH
+
+#include <vector>
+#include <utility>
+#include <type_traits>
+
+#include <dumux/common/typetraits/isvalid.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup CCMpfaDiscretization
+ * \brief A class that contains helper functions as well as functionality
+ * which is common to different mpfa schemes and which solely
+ * operate on the interaction volume.
+ */
+class InteractionVolumeAssemblerHelper
+{
+ // Helper structs to detect if matrix has resize function
+ struct HasMatrixResize
+ {
+ template<class M>
+ auto operator()(const M& m) -> decltype(std::declval<M>().resize(0, 0))
+ {}
+ };
+
+ // Helper structs to detect if vector has resize function
+ struct HasVectorResize
+ {
+ template<class V>
+ auto operator()(const V& v) -> decltype(std::declval<V>().resize(0))
+ {}
+ };
+
+ template<class Matrix>
+ static constexpr auto matrixHasResizeFunction()
+ { return decltype( isValid(HasMatrixResize())(std::declval<Matrix>()) )::value; }
+
+ template<class Vector>
+ static constexpr auto vectorHasResizeFunction()
+ { return decltype( isValid(HasVectorResize())(std::declval<Vector>()) )::value; }
+
+public:
+ /*!
+ * \brief Assembles the vector of face unknowns within an interaction volume.
+ * \note This requires the data handle to be fully assembled already.
+ *
+ * \param handle The data handle in which the vector is stored
+ * \param iv The interaction volume
+ */
+ template< class DataHandle, class IV >
+ static typename IV::Traits::MatVecTraits::FaceVector
+ assembleFaceUnkowns(const DataHandle& handle, const IV& iv)
+ {
+ typename IV::Traits::MatVecTraits::FaceVector u;
+ resizeVector(u, iv.numFaces());
+
+ handle.AB().mv(handle.uj(), u);
+
+ // maybe add gravity terms
+ if (handle.deltaG().size() == iv.numUnknowns())
+ handle.AB().umv(handle.deltaG(), u);
+
+ return u;
+ }
+
+ /*!
+ * \brief Assembles the solution gradients in the
+ * sub-control volumes within an interaction volume.
+ * \note This requires the data handle to be fully assembled already.
+ *
+ * \param handle The data handle in which the vector is stored
+ * \param iv The interaction volume
+ */
+ template< class DataHandle, class IV >
+ static std::vector< typename IV::Traits::LocalScvType::GlobalCoordinate >
+ assembleScvGradients(const DataHandle& handle, const IV& iv)
+ {
+ const auto u = assembleFaceUnkowns(handle, iv);
+
+ using LocalScv = typename IV::Traits::LocalScvType;
+ using Gradient = typename LocalScv::GlobalCoordinate;
+
+ std::vector<Gradient> result; result.reserve(iv.numScvs());
+ for (unsigned int scvIdx = 0; scvIdx < iv.numScvs(); ++scvIdx)
+ {
+ const auto& scv = iv.localScv(scvIdx);
+
+ Gradient gradU(0.0);
+ for (unsigned int dir = 0; dir < LocalScv::myDimension; ++dir)
+ {
+ auto nu = scv.nu(dir);
+
+ // obtain face pressure
+ const auto& scvf = iv.localScvf( scv.localScvfIndex(dir) );
+ const auto faceU = !scvf.isDirichlet() ? u[scvf.localDofIndex()]
+ : handle.uj()[scvf.localDofIndex()];
+
+ nu *= faceU - handle.uj()[scv.localDofIndex()];
+ gradU += nu;
+ }
+
+ gradU /= scv.detX();
+ result.emplace_back( std::move(gradU) );
+ }
+
+ return result;
+ }
+
+ //! resizes a matrix to the given sizes (specialization for dynamic matrix type)
+ template< class Matrix,
+ class size_type,
+ std::enable_if_t<matrixHasResizeFunction<Matrix>(), int> = 0 >
+ static void resizeMatrix(Matrix& M, size_type rows, size_type cols)
+ {
+ M.resize(rows, cols);
+ }
+
+ //! resizes a matrix to the given sizes (specialization for static matrix type - do nothing)
+ template< class Matrix,
+ class size_type,
+ std::enable_if_t<!matrixHasResizeFunction<Matrix>(), int> = 0 >
+ static void resizeMatrix(Matrix& M, size_type rows, size_type cols)
+ {}
+
+ //! resizes a vector to the given size (specialization for dynamic matrix type)
+ template< class Vector,
+ class size_type,
+ std::enable_if_t<vectorHasResizeFunction<Vector>(), int> = 0 >
+ static void resizeVector(Vector& v, size_type size)
+ {
+ v.resize(size);
+ }
+
+ //! resizes a vector to the given size (specialization for static vector type - do nothing)
+ template< class Vector,
+ class size_type,
+ std::enable_if_t<!vectorHasResizeFunction<Vector>(), int> = 0 >
+ static void resizeVector(Vector& v, size_type rows)
+ {}
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
index a4039645784954e00004fc1caf71c3dc756fedc6..0843c7f6f2f9320fda38f38122a418cf89bc7f4a 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
@@ -30,6 +30,7 @@
#include <dumux/common/math.hh>
#include <dumux/discretization/cellcentered/mpfa/methods.hh>
#include <dumux/discretization/cellcentered/mpfa/localassemblerbase.hh>
+#include <dumux/discretization/cellcentered/mpfa/localassemblerhelper.hh>
#include <dumux/discretization/cellcentered/mpfa/computetransmissibility.hh>
namespace Dumux
@@ -49,6 +50,7 @@ class MpfaOInteractionVolumeAssembler
: public InteractionVolumeAssemblerBase< P, EG, EV >
{
using ParentType = InteractionVolumeAssemblerBase< P, EG, EV >;
+ using Helper = InteractionVolumeAssemblerHelper;
public:
//! Pull up constructor of the base class
@@ -99,7 +101,7 @@ public:
tijOut[fIdx].resize(numOutsideFaces);
std::for_each(tijOut[fIdx].begin(),
tijOut[fIdx].end(),
- [&](auto& v) { this->resizeVector_(v, iv.numKnowns()); });
+ [&](auto& v) { Helper::resizeVector(v, iv.numKnowns()); });
}
// compute outside transmissibilities
@@ -155,7 +157,7 @@ public:
void assembleU(DataHandle& handle, const IV& iv, const GetU& getU)
{
auto& u = handle.uj();
- this->resizeVector_(u, iv.numKnowns());
+ Helper::resizeVector(u, iv.numKnowns());
// put the cell unknowns first, then Dirichlet values
typename IV::Traits::IndexSet::LocalIndexType i = 0;
@@ -202,14 +204,14 @@ private:
static constexpr int dimWorld = IV::Traits::GridView::dimensionworld;
// resize omegas
- this->resizeVector_(wijk, iv.numFaces());
+ Helper::resizeVector(wijk, iv.numFaces());
// if only Dirichlet faces are present in the iv,
// the matrices A, B & C are undefined and D = T
if (iv.numUnknowns() == 0)
{
// resize & reset D matrix
- this->resizeMatrix_(D, iv.numFaces(), iv.numKnowns()); D = 0.0;
+ Helper::resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0;
// Loop over all the faces, in this case these are all dirichlet boundaries
for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
@@ -226,7 +228,7 @@ private:
const auto tensor = getT(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
// the omega factors of the "positive" sub volume
- this->resizeVector_(wijk[faceIdx], /*no outside scvs present*/1);
+ Helper::resizeVector(wijk[faceIdx], /*no outside scvs present*/1);
wijk[faceIdx][0] = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor());
const auto posScvLocalDofIdx = posLocalScv.localDofIndex();
@@ -242,10 +244,10 @@ private:
else
{
// resize & reset matrices
- this->resizeMatrix_(A, iv.numUnknowns(), iv.numUnknowns()); A = 0.0;
- this->resizeMatrix_(B, iv.numUnknowns(), iv.numKnowns()); B = 0.0;
- this->resizeMatrix_(C, iv.numFaces(), iv.numUnknowns()); C = 0.0;
- this->resizeMatrix_(D, iv.numFaces(), iv.numKnowns()); D = 0.0;
+ Helper::resizeMatrix(A, iv.numUnknowns(), iv.numUnknowns()); A = 0.0;
+ Helper::resizeMatrix(B, iv.numUnknowns(), iv.numKnowns()); B = 0.0;
+ Helper::resizeMatrix(C, iv.numFaces(), iv.numUnknowns()); C = 0.0;
+ Helper::resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0;
for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
{
@@ -263,7 +265,7 @@ private:
const auto tensor = getT(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv);
// the omega factors of the "positive" sub volume
- this->resizeVector_(wijk[faceIdx], neighborScvIndices.size());
+ Helper::resizeVector(wijk[faceIdx], neighborScvIndices.size());
wijk[faceIdx][0] = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor());
// go over the coordinate directions in the positive sub volume