diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh index 9069ebeff39e87f9634115bfdd2e6e38b936a617..01d8be6c03fbb1763b8549f722d3797dac7d7b74 100644 --- a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh +++ b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh @@ -141,13 +141,14 @@ public: const auto& problem = globalFluxVarsCache().problem_(); const auto& globalFvGeometry = problem.model().globalFvGeometry(); - const auto& neighborStencil = problem.model().stencils(element).neighborStencil(); - const auto& assemblyMap = problem.model().localJacobian().assemblyMap(); + const auto globalI = problem.elementMapper().index(element); + const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI]; // reserve memory - auto numNeighborScvfs = 0; - for (auto&& facesInNeighbor : assemblyMap[globalI]) numNeighborScvfs += facesInNeighbor.size(); + unsigned int numNeighborScvfs = 0; + for (auto&& dataJ : assemblyMapI) + numNeighborScvfs += dataJ.scvfsJ.size(); globalScvfIndices_.reserve(fvGeometry.numScvf() + numNeighborScvfs); // first add all the indices inside the element @@ -155,9 +156,9 @@ public: globalScvfIndices_.push_back(scvf.index()); // for the indices in the neighbors, use assembly map of the local jacobian - for (unsigned int j = 0; j < neighborStencil.size(); ++j) - for (auto fluxVarIdx : assemblyMap[globalI][j]) - globalScvfIndices_.push_back(fluxVarIdx); + for (auto&& dataJ : assemblyMapI) + for (auto scvfIdx : dataJ.scvfsJ) + globalScvfIndices_.push_back(scvfIdx); // prepare all the caches of the scvfs inside the corresponding interaction volumes using helper class fluxVarsCache_.resize(globalScvfIndices_.size()); @@ -166,21 +167,17 @@ public: FluxVariablesCacheFiller::fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, *this); // prepare the caches in the remaining neighbors - unsigned int j = 0; - for (auto globalJ : neighborStencil) + for (auto&& dataJ : assemblyMapI) { - for (auto fluxVarIdx : assemblyMap[globalI][j]) + for (auto scvfIdx : dataJ.scvfsJ) { - const auto& scvf = fvGeometry.scvf(fluxVarIdx); + const auto& scvf = fvGeometry.scvf(scvfIdx); if (!(*this)[scvf].isUpdated()) { - auto elementJ = globalFvGeometry.element(globalJ); + auto elementJ = globalFvGeometry.element(dataJ.globalJ); FluxVariablesCacheFiller::fillFluxVarCache(problem, elementJ, fvGeometry, elemVolVars, scvf, *this); } } - - // increment counter - j++; } } diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh index 3d20dfa2e659bec72f58d0369a570b3f04bcbb45..1046a21ce7d6da215674c4a7088a50973ca1dc6a 100644 --- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh @@ -121,8 +121,9 @@ public: const auto& globalFvGeometry = problem.model().globalFvGeometry(); // stencil information - const auto& neighborStencil = problem.model().stencils(element).neighborStencil(); - const auto numDofs = neighborStencil.size() + 1; + const auto globalI = problem.elementMapper().index(element); + const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI]; + const auto numDofs = assemblyMapI.size() + 1; // resize local containers to the required size (for internal elements) volumeVariables_.resize(numDofs); @@ -140,10 +141,10 @@ public: ++localIdx; // Update the volume variables of the neighboring elements - for (auto globalJ : neighborStencil) + for (auto&& dataJ : assemblyMapI) { - const auto& elementJ = globalFvGeometry.element(globalJ); - auto&& scvJ = fvGeometry.scv(globalJ); + const auto& elementJ = globalFvGeometry.element(dataJ.globalJ); + auto&& scvJ = fvGeometry.scv(dataJ.globalJ); volumeVariables_[localIdx].update(problem.model().elementSolution(elementJ, sol), problem, elementJ, diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh index b7e74e6b686f246f5f5aec62a14690e2f890070a..c5bf3269a435203990c94e5da7ab3f89e30f281c 100644 --- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh +++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh @@ -256,30 +256,23 @@ public: // get some references for convenience const auto& problem = globalFvGeometry().problem_(); - const auto& assemblyMap = problem.model().localJacobian().assemblyMap(); - const auto& neighborStencil = problem.model().stencils(element).neighborStencil(); const auto globalI = problem.elementMapper().index(element); + const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI]; // reserve memory - const auto numNeighbors = neighborStencil.size(); + const auto numNeighbors = assemblyMapI.size(); const auto estimate = neighborScvfEstimate_(element, numNeighbors); neighborScvs_.reserve(numNeighbors); neighborScvIndices_.reserve(numNeighbors); neighborScvfIndices_.reserve(estimate); neighborScvfs_.reserve(estimate); - // build scvfs in neighbors + // make neighbor geometries // use the assembly map to determine which faces are necessary - unsigned int j = 0; - for (auto globalJ : neighborStencil) - { - // make neighbor geometries - auto elementJ = globalFvGeometry().element(globalJ); - makeNeighborGeometries(elementJ, globalJ, assemblyMap[globalI][j]); - - // increment counter - j++; - } + for (auto&& dataJ : assemblyMapI) + makeNeighborGeometries(globalFvGeometry().element(dataJ.globalJ), + dataJ.globalJ, + dataJ.scvfsJ); // set up the scvf flip indices for network grids if (dim < dimWorld) diff --git a/dumux/discretization/cellcentered/stencils.hh b/dumux/discretization/cellcentered/stencils.hh deleted file mode 100644 index 49d70b767ad9958a3c2c7758307c6103c9ea0363..0000000000000000000000000000000000000000 --- a/dumux/discretization/cellcentered/stencils.hh +++ /dev/null @@ -1,259 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \file - * \brief Implements the notion of stencils for cell-centered models - */ -#ifndef DUMUX_DISCRETIZATION_CC_STENCILS_HH -#define DUMUX_DISCRETIZATION_CC_STENCILS_HH - -#include - -namespace Dumux -{ - -/*! - * \brief Element-related stencils for symmetric cc methods. - */ -template -class CCSymmetricElementStencils -{ - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using IndexType = typename GridView::IndexSet::IndexType; - using Element = typename GridView::template Codim<0>::Entity; - // TODO a separate stencil class all stencils can derive from? - using Stencil = std::vector; -public: - //! Update the stencil. - void update(const Problem& problem, - const Element& element) - { - // restrict the FvGeometry locally and bind to the element - auto fvGeometry = localView(problem.model().globalFvGeometry()); - fvGeometry.bindElement(element); // for TPFA bind element is enough - - elementStencil_.clear(); - - // loop over sub control faces - for (auto&& scvf : scvfs(fvGeometry)) - { - FluxVariables fluxVars; - const auto& stencil = fluxVars.computeStencil(problem, element, fvGeometry, scvf); - elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end()); - } - // make values in elementstencil unique - std::sort(elementStencil_.begin(), elementStencil_.end()); - elementStencil_.erase(std::unique(elementStencil_.begin(), elementStencil_.end()), elementStencil_.end()); - - auto globalI = problem.elementMapper().index(element); - neighborStencil_ = elementStencil_; - - // remove the element itself and possible ghost neighbors from the neighbor stencil - neighborStencil_.erase(std::remove_if(neighborStencil_.begin(), neighborStencil_.end(), - [globalI](auto i){ return (i == globalI); }), - neighborStencil_.end()); - } - - //! The full element stencil (all element this element is interacting with) - const Stencil& elementStencil() const - { - return elementStencil_; - } - - //! The full element stencil without this element - const Stencil& neighborStencil() const - { - return neighborStencil_; - } - -private: - Stencil elementStencil_; - Stencil neighborStencil_; -}; - -/*! - * \ingroup CCModel - * \brief The global stencil container class for symmetric cc methods - */ -template -class CCSymmetricStencilsVector -{ - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using StencilType = CCSymmetricElementStencils; - -public: - void update(Problem& problem) - { - problemPtr_ = &problem; - elementStencils_.resize(problem.gridView().size(0)); - for (const auto& element : elements(problem.gridView())) - { - auto eIdx = problem.elementMapper().index(element); - elementStencils_[eIdx].update(problem, element); - } - } - - //! overload for elements - template - typename std::enable_if::type - get(const Entity& entity) const - { - return elementStencils_[problemPtr_->elementMapper().index(entity)]; - } - -private: - std::vector elementStencils_; - const Problem* problemPtr_; -}; - -//! Forward declaration of the global stencil class for nonsymmetric cc methods -template class CCNonSymmetricStencilsVector; - -/*! - * \brief Element-related stencils for non-symmetric cc methods - */ -template -class CCNonSymmetricElementStencils -{ - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using IndexType = typename GridView::IndexSet::IndexType; - using Element = typename GridView::template Codim<0>::Entity; - // TODO a separate stencil class all stencils can derive from? - using Stencil = std::vector; -public: - //! Update the stencil. We expect a bound fvGeometry - void update(const Problem& problem, - const Element& element, - CCNonSymmetricStencilsVector& stencilsVector) - { - // the element index of our own here - auto eIdx = problem.elementMapper().index(element); - elementStencil_.clear(); - - // restrict the FvGeometry locally and bind to the element - auto fvGeometry = localView(problem.model().globalFvGeometry()); - fvGeometry.bindElement(element); - - // loop over sub control faces - for (auto&& scvf : scvfs(fvGeometry)) - { - FluxVariables fluxVars; - const auto& stencil = fluxVars.computeStencil(problem, element, fvGeometry, scvf); - - // insert stencil into the element stencil - elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end()); - - // insert our index in the neighbor stencils of the elements in the flux stencil - for (auto globalI : stencil) - { - if (globalI != eIdx) - stencilsVector[globalI].addNeighbor(eIdx); - } - } - - // make values in element stencil unique - std::sort(elementStencil_.begin(), elementStencil_.end()); - elementStencil_.erase(std::unique(elementStencil_.begin(), elementStencil_.end()), elementStencil_.end()); - } - - //! The full element stencil (all element this element is interacting with) - const Stencil& elementStencil() const - { - return elementStencil_; - } - - //! The full element stencil without this element - const Stencil& neighborStencil() const - { - return neighborStencil_; - } - - void addNeighbor(const IndexType nIdx) - { - neighborStencil_.push_back(nIdx); - } - - void makeNeighborStencilUnique() - { - // make values in neighbor stencil unique - std::sort(neighborStencil_.begin(), neighborStencil_.end()); - neighborStencil_.erase(std::unique(neighborStencil_.begin(), neighborStencil_.end()), neighborStencil_.end()); - } - -private: - Stencil elementStencil_; - Stencil neighborStencil_; -}; - -/*! - * \ingroup CCModel - * \brief The global stencil container class for non-symmetric cc methods - */ -template -class CCNonSymmetricStencilsVector -{ - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using IndexType = typename GridView::IndexSet::IndexType; - using StencilType = CCNonSymmetricElementStencils; - - // The stencil type of an element has to access bracket operator - friend StencilType; - -public: - void update(Problem& problem) - { - problemPtr_ = &problem; - elementStencils_.resize(problem.gridView().size(0)); - for (const auto& element : elements(problem.gridView())) - { - auto eIdx = problem.elementMapper().index(element); - elementStencils_[eIdx].update(problem, element, *this); - } - - for (auto&& stencil : elementStencils_) - stencil.makeNeighborStencilUnique(); - } - - //! overload for elements - template - typename std::enable_if::type - get(const Entity& entity) const - { - return elementStencils_[problemPtr_->elementMapper().index(entity)]; - } - -private: - StencilType& operator[] (const IndexType globalIdx) - { return elementStencils_[globalIdx]; } - - std::vector elementStencils_; - const Problem* problemPtr_; -}; - -} // end namespace - -#endif diff --git a/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh index 0636ba0478e4380f37ef6608aad0d2abd4e6dc99..343e7d13e7f24c10c755a36d30cfde0813cd7b41 100644 --- a/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh +++ b/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh @@ -145,39 +145,35 @@ public: const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) { - const auto globalI = globalFluxVarsCache().problem_().elementMapper().index(element); - const auto& neighborStencil = globalFluxVarsCache().problem_().model().stencils(element).neighborStencil(); - const auto numNeighbors = neighborStencil.size(); + const auto& problem = globalFluxVarsCache().problem_(); + const auto globalI = problem.elementMapper().index(element); + const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI]; + const auto numNeighbors = assemblyMapI.size(); // find the number of scv faces that need to be prepared auto numScvf = fvGeometry.numScvf(); - for (IndexType localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ) - { - const auto& fluxVarIndicesJ = globalFluxVarsCache().problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ]; - numScvf += fluxVarIndicesJ.size(); - } + for (unsigned int localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ) + numScvf += assemblyMapI[localIdxJ].scvfsJ.size(); // fill the containers with the data on the scv faces inside the actual element fluxVarsCache_.resize(numScvf); globalScvfIndices_.resize(numScvf); - IndexType localScvfIdx = 0; + unsigned int localScvfIdx = 0; for (auto&& scvf : scvfs(fvGeometry)) { - fluxVarsCache_[localScvfIdx].update(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf); + fluxVarsCache_[localScvfIdx].update(problem, element, fvGeometry, elemVolVars, scvf); globalScvfIndices_[localScvfIdx] = scvf.index(); localScvfIdx++; } // add required data on the scv faces in the neighboring elements - for (IndexType localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ) + for (unsigned int localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ) { - const auto& fluxVarIndicesJ = globalFluxVarsCache().problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ]; - const auto elementJ = fvGeometry.globalFvGeometry().element(neighborStencil[localIdxJ]); - - for (auto fluxVarIdx : fluxVarIndicesJ) + const auto elementJ = fvGeometry.globalFvGeometry().element(assemblyMapI[localIdxJ].globalJ); + for (auto scvfIdx : assemblyMapI[localIdxJ].scvfsJ) { - auto&& scvfJ = fvGeometry.scvf(fluxVarIdx); - fluxVarsCache_[localScvfIdx].update(globalFluxVarsCache().problem_(), elementJ, fvGeometry, elemVolVars, scvfJ); + auto&& scvfJ = fvGeometry.scvf(scvfIdx); + fluxVarsCache_[localScvfIdx].update(problem, elementJ, fvGeometry, elemVolVars, scvfJ); globalScvfIndices_[localScvfIdx] = scvfJ.index(); localScvfIdx++; } diff --git a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh index 9080810de0da51edff5bf601be56655e45b600dd..6cab6287d8c50ac5151b7839e557f0f22aba9463 100644 --- a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh @@ -117,11 +117,10 @@ public: const FVElementGeometry& fvGeometry, const SolutionVector& sol) { - auto eIdx = globalVolVars().problem_().elementMapper().index(element); - - // stencil information - const auto& neighborStencil = globalVolVars().problem_().model().stencils(element).neighborStencil(); - const auto numDofs = neighborStencil.size() + 1; + const auto& problem = globalVolVars().problem_(); + const auto globalI = problem.elementMapper().index(element); + const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI]; + const auto numDofs = assemblyMapI.size() + 1; // resize local containers to the required size (for internal elements) volumeVariables_.resize(numDofs); @@ -129,21 +128,21 @@ public: int localIdx = 0; // update the volume variables of the element at hand - auto&& scvI = fvGeometry.scv(eIdx); - volumeVariables_[localIdx].update(globalVolVars().problem_().model().elementSolution(element, sol), - globalVolVars().problem_(), + auto&& scvI = fvGeometry.scv(globalI); + volumeVariables_[localIdx].update(problem.model().elementSolution(element, sol), + problem, element, scvI); volVarIndices_[localIdx] = scvI.index(); ++localIdx; // Update the volume variables of the neighboring elements - for (auto globalJ : neighborStencil) + for (const auto& dataJ : assemblyMapI) { - const auto& elementJ = fvGeometry.globalFvGeometry().element(globalJ); - auto&& scvJ = fvGeometry.scv(globalJ); - volumeVariables_[localIdx].update(globalVolVars().problem_().model().elementSolution(elementJ, sol), - globalVolVars().problem_(), + const auto& elementJ = fvGeometry.globalFvGeometry().element(dataJ.globalJ); + auto&& scvJ = fvGeometry.scv(dataJ.globalJ); + volumeVariables_[localIdx].update(problem.model().elementSolution(elementJ, sol), + problem, elementJ, scvJ); volVarIndices_[localIdx] = scvJ.index(); @@ -158,15 +157,15 @@ public: continue; // check if boundary is a pure dirichlet boundary - const auto bcTypes = globalVolVars().problem_().boundaryTypes(element, scvf); + const auto bcTypes = problem.boundaryTypes(element, scvf); if (bcTypes.hasOnlyDirichlet()) { - const ElementSolution dirichletPriVars({globalVolVars().problem_().dirichlet(element, scvf)}); + const ElementSolution dirichletPriVars({problem.dirichlet(element, scvf)}); volumeVariables_.resize(localIdx+1); volVarIndices_.resize(localIdx+1); volumeVariables_[localIdx].update(dirichletPriVars, - globalVolVars().problem_(), + problem, element, scvI); volVarIndices_[localIdx] = scvf.outsideScvIdx(); diff --git a/dumux/implicit/assembler.hh b/dumux/implicit/assembler.hh index f54a127bbec200ca7c730627d7f5f968f9197f3b..4fc310e2bcd9750bf4bdadb7a40a479bfcaaec8b 100644 --- a/dumux/implicit/assembler.hh +++ b/dumux/implicit/assembler.hh @@ -176,27 +176,7 @@ private: // Construct the BCRS matrix for the global jacobian void createMatrix_() { - auto numDofs = problem_().model().numDofs(); - - // allocate raw matrix - matrix_ = std::make_shared(numDofs, numDofs, JacobianMatrix::random); - - // set the row sizes - asImp_().setRowSizes_(); - - // set the indices - asImp_().addIndices_(); - } - - //! Set the row sizes - void setRowSizes_() - { - DUNE_THROW(Dune::NotImplemented, "Actual implementation does not provide a setRowSizes_() method!"); - } - - void addIndices_() - { - DUNE_THROW(Dune::NotImplemented, "Actual implementation does not provide a addIndices_() method!"); + DUNE_THROW(Dune::NotImplemented, "Actual implementation does not provide a createMatrix_() method!"); } Implementation &asImp_() diff --git a/dumux/implicit/box/assembler.hh b/dumux/implicit/box/assembler.hh index 33a18bc25680505177cc1a57cd18a004d10b586f..fa34084202f19163adf2d7863de507ca4b384b7e 100644 --- a/dumux/implicit/box/assembler.hh +++ b/dumux/implicit/box/assembler.hh @@ -23,6 +23,8 @@ #ifndef DUMUX_BOX_ASSEMBLER_HH #define DUMUX_BOX_ASSEMBLER_HH +#include + #include #include @@ -35,40 +37,37 @@ namespace Dumux { template class BoxAssembler : public ImplicitAssembler { - typedef ImplicitAssembler ParentType; friend class ImplicitAssembler; - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - - typedef typename GridView::IndexSet::IndexType IndexType; + using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + static const int dim = GridView::dimension; - void setRowSizes_() + void createMatrix_() { - for (const auto& vertex : vertices(this->gridView_())) - { - // the global index of the element at hand - const auto globalI = this->vertexMapper_().index(vertex); - const auto& stencil = this->model_().stencils(vertex).vertexStencil(); + const auto numDofs = this->problem_().model().numDofs(); - this->matrix().setrowsize(globalI, stencil.size()); - } - this->matrix().endrowsizes(); + // allocate raw matrix + this->matrix_ = std::make_shared(numDofs, numDofs, JacobianMatrix::random); - } + // get occupation pattern of the matrix + Dune::MatrixIndexSet occupationPattern; + occupationPattern.resize(numDofs, numDofs); - void addIndices_() - { - for (const auto& vertex : vertices(this->gridView_())) + for (const auto& element : elements(this->problem_().gridView())) { - // the global index of the element at hand - const auto globalI = this->vertexMapper_().index(vertex); - const auto& stencil = this->model_().stencils(vertex).vertexStencil(); - - for (auto&& globalJ : stencil) - this->matrix().addindex(globalI, globalJ); + for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx) + { + const auto globalI = this->problem_().vertexMapper().subIndex(element, vIdx, dim); + for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(dim); ++vIdx2) + { + const auto globalJ = this->problem_().vertexMapper().subIndex(element, vIdx2, dim); + occupationPattern.add(globalI, globalJ); + occupationPattern.add(globalJ, globalI); + } + } } - this->matrix().endindices(); + // export pattern to matrix + occupationPattern.exportIdx(*this->matrix_); } }; diff --git a/dumux/implicit/box/propertydefaults.hh b/dumux/implicit/box/propertydefaults.hh index bc26b2565a5d18edad0dc7497631b835985feb3f..c95a65d75be728c9371637ccb9082bb1cd46bfe8 100644 --- a/dumux/implicit/box/propertydefaults.hh +++ b/dumux/implicit/box/propertydefaults.hh @@ -35,7 +35,6 @@ #include #include #include -#include #include #include @@ -53,7 +52,6 @@ namespace Dumux { // forward declarations template class BoxLocalResidual; template class BoxElementBoundaryTypes; -template class BoxStencilsVector; namespace Properties { //! Set the corresponding discretization method property @@ -101,9 +99,6 @@ SET_TYPE_PROP(BoxModel, ElementBoundaryTypes, BoxElementBoundaryTypes); //! Mapper for the degrees of freedoms. SET_TYPE_PROP(BoxModel, DofMapper, typename GET_PROP_TYPE(TypeTag, VertexMapper)); -//! The stencil container -SET_TYPE_PROP(BoxModel, StencilsVector, BoxStencilsVector); - //! The global volume variables vector class SET_TYPE_PROP(BoxModel, GlobalVolumeVariables, BoxGlobalVolumeVariables); diff --git a/dumux/implicit/cellcentered/assembler.hh b/dumux/implicit/cellcentered/assembler.hh index a2430d7cbd53139bd1ea52a0b55de4a03b8b92fd..d376a026a586effaa52486c9582d9b0c00d020d2 100644 --- a/dumux/implicit/cellcentered/assembler.hh +++ b/dumux/implicit/cellcentered/assembler.hh @@ -23,6 +23,8 @@ #ifndef DUMUX_CC_ASSEMBLER_HH #define DUMUX_CC_ASSEMBLER_HH +#include + #include #include @@ -35,40 +37,30 @@ namespace Dumux { template class CCAssembler : public ImplicitAssembler { - typedef ImplicitAssembler ParentType; friend class ImplicitAssembler; - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GridView::IndexSet::IndexType IndexType; + using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix); - void setRowSizes_() + void createMatrix_() { - for (const auto& element : elements(this->gridView_())) - { - // the global index of the element at hand - const auto globalI = this->elementMapper_().index(element); - const auto& stencil = this->model_().stencils(element).elementStencil(); + const auto numDofs = this->problem_().model().numDofs(); - this->matrix().setrowsize(globalI, stencil.size()); - } - this->matrix().endrowsizes(); - } + // allocate raw matrix + this->matrix_ = std::make_shared(numDofs, numDofs, JacobianMatrix::random); - void addIndices_() - { - for (const auto& element : elements(this->gridView_())) - { - // the global index of the element at hand - const auto globalI = this->elementMapper_().index(element); - const auto& stencil = this->model_().stencils(element).elementStencil(); + // get occupation pattern of the matrix + Dune::MatrixIndexSet occupationPattern; + occupationPattern.resize(numDofs, numDofs); - for (auto&& globalJ : stencil) - this->matrix().addindex(globalI, globalJ); + const auto& assemblyMap = this->problem_().model().localJacobian().assemblyMap(); + for (unsigned int globalI = 0; globalI < numDofs; ++globalI) + { + occupationPattern.add(globalI, globalI); + for (const auto& dataJ : assemblyMap[globalI]) + occupationPattern.add(dataJ.globalJ, globalI); } - this->matrix().endindices(); + + // export pattern to matrix + occupationPattern.exportIdx(*this->matrix_); } }; diff --git a/dumux/implicit/cellcentered/assemblymap.hh b/dumux/implicit/cellcentered/assemblymap.hh index 4bce9c22cb83d305c9910b79a3912e803f3d3ce0..6e8e1d78cafa6471bef1eefda9f48075346c89e9 100644 --- a/dumux/implicit/cellcentered/assemblymap.hh +++ b/dumux/implicit/cellcentered/assemblymap.hh @@ -19,26 +19,35 @@ /*! * \file * \brief Stores the face indices corresponding to the neighbors of an element - * that contribute to the derivative calculation + * that contribute to the derivative calculation. This is used for + * finite-volume schemes with symmetric sparsity pattern in the global matrix. */ -#ifndef DUMUX_CC_ASSEMBLY_MAP_HH -#define DUMUX_CC_ASSEMBLY_MAP_HH +#ifndef DUMUX_CC_SYMMETRIC_ASSEMBLY_MAP_HH +#define DUMUX_CC_SYMMETRIC_ASSEMBLY_MAP_HH -#include +#include +#include namespace Dumux { template -class CCAssemblyMap +class CCSymmetricAssemblyMap { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); using IndexType = typename GridView::IndexSet::IndexType; - using Map = std::vector< std::vector< std::vector > >; + + struct DataJ + { + IndexType globalJ; + std::vector scvfsJ; + }; + + using Map = std::vector>; public: @@ -52,58 +61,51 @@ public: map_.resize(problem.gridView().size(0)); for (const auto& element : elements(problem.gridView())) { - // get a local finite volume geometry object that is bindable - auto fvGeometryJ = localView(problem.model().globalFvGeometry()); + // We are looking for the elements I, for which this element J is in the flux stencil + auto globalJ = problem.elementMapper().index(element); - auto globalI = problem.elementMapper().index(element); - const auto& neighborStencil = problem.model().stencils(element).neighborStencil(); + auto fvGeometry = localView(problem.model().globalFvGeometry()); + fvGeometry.bindElement(element); - map_[globalI].reserve(neighborStencil.size()); - for (auto globalJ : neighborStencil) - { - const auto& elementJ = fvGeometryJ.globalFvGeometry().element(globalJ); + // obtain the data of J in elements I + std::vector> dataJForI; - // find the flux vars needed for the calculation of the flux into element - std::vector fluxVarIndices; + // loop over sub control faces + for (auto&& scvf : scvfs(fvGeometry)) + { + FluxVariables fluxVars; + const auto& stencil = fluxVars.computeStencil(problem, element, fvGeometry, scvf); - // only non-ghost neighbors (J) have to be considered, derivatives from non-ghost to ghost dofs - // are assembled when assembling the ghost element (I) - if (elementJ.partitionType() != Dune::GhostEntity) + // insert our index in the neighbor stencils of the elements in the flux stencil + for (auto globalI : stencil) { - fvGeometryJ.bindElement(elementJ); - for (auto&& scvFaceJ : scvfs(fvGeometryJ)) - { - auto fluxVarsIdx = scvFaceJ.index(); - - // if globalI is in flux var stencil, add to list - FluxVariables fluxVars; - const auto fluxStencil = fluxVars.computeStencil(problem, elementJ, fvGeometryJ, scvFaceJ); - - for (auto globalIdx : fluxStencil) - { - if (globalIdx == globalI) - { - fluxVarIndices.push_back(fluxVarsIdx); - break; - } - } - } + if (globalI == globalJ) + continue; + + auto it = std::find_if(dataJForI.begin(), + dataJForI.end(), + [globalI](const auto& pair) { return pair.first == globalI; }); + + if (it != dataJForI.end()) + it->second.scvfsJ.push_back(scvf.index()); + else + dataJForI.emplace_back(std::make_pair(globalI, + DataJ({globalJ, std::vector({scvf.index()})}))); } - map_[globalI].emplace_back(std::move(fluxVarIndices)); } + + for (auto&& pair : dataJForI) + map_[pair.first].emplace_back(std::move(pair.second)); } } - //! Some implementations of the assembly map can be solution-dependent and need to be updated - void update() {} - - const std::vector>& operator [] (const IndexType globalI) const + const std::vector& operator[] (const IndexType globalI) const { return map_[globalI]; } private: Map map_; }; -} +} // end namespace Dumux #endif diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index 6ac50bf3f912c34061f02a2e070772026a6e9bd8..7a5fafd63b768d86564c608da3cddc01abb2ebf6 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -211,8 +211,7 @@ private: const bool isGhost) { // get stencil informations - const auto& neighborStencil = this->model_().stencils(element).neighborStencil(); - const auto numNeighbors = neighborStencil.size(); + const auto numNeighbors = assemblyMap_[globalI_].size(); // the localresidual class used for the flux calculations LocalResidual localRes; @@ -227,14 +226,14 @@ private: origFlux = 0.0; unsigned int j = 0; - for (auto globalJ : neighborStencil) + for (const auto& dataJ : assemblyMap_[globalI_]) { - auto elementJ = fvGeometry.globalFvGeometry().element(globalJ); + auto elementJ = fvGeometry.globalFvGeometry().element(dataJ.globalJ); neighborElements.push_back(elementJ); - for (auto fluxVarIdx : assemblyMap_[globalI_][j]) + for (auto scvfIdx : dataJ.scvfsJ) { - auto&& scvf = fvGeometry.scvf(fluxVarIdx); + auto&& scvf = fvGeometry.scvf(scvfIdx); origFlux[j] += localRes.evalFlux_(elementJ, fvGeometry, curElemVolVars, @@ -295,9 +294,9 @@ private: // calculate the fluxes into element with the deflected primary variables for (std::size_t k = 0; k < numNeighbors; ++k) { - for (auto fluxVarIdx : assemblyMap_[globalI_][k]) + for (auto scvfIdx : assemblyMap_[globalI_][k].scvfsJ) { - auto&& scvf = fvGeometry.scvf(fluxVarIdx); + auto&& scvf = fvGeometry.scvf(scvfIdx); neighborDeriv[k] += localRes.evalFlux_(neighborElements[k], fvGeometry, curElemVolVars, @@ -346,9 +345,9 @@ private: // calculate the fluxes into element with the deflected primary variables for (std::size_t k = 0; k < numNeighbors; ++k) { - for (auto fluxVarIdx : assemblyMap_[globalI_][k]) + for (auto scvfIdx : assemblyMap_[globalI_][k].scvfsJ) { - auto&& scvf = fvGeometry.scvf(fluxVarIdx); + auto&& scvf = fvGeometry.scvf(scvfIdx); neighborDeriv[k] -= localRes.evalFlux_(neighborElements[k], fvGeometry, curElemVolVars, @@ -383,8 +382,8 @@ private: this->updateGlobalJacobian_(matrix, globalI_, globalI_, pvIdx, partialDeriv); j = 0; - for (auto globalJ : neighborStencil) - this->updateGlobalJacobian_(matrix, globalJ, globalI_, pvIdx, neighborDeriv[j++]); + for (const auto& dataJ : assemblyMap_[globalI_]) + this->updateGlobalJacobian_(matrix, dataJ.globalJ, globalI_, pvIdx, neighborDeriv[j++]); } } diff --git a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh index 9ff21ee167e7f90fa00b8933d457bc31eea491d6..3278dafaef82cb8bda520f5522d7cd82cacf0ff0 100644 --- a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh @@ -39,7 +39,6 @@ #include #include #include -#include #include #include #include @@ -93,9 +92,6 @@ SET_TYPE_PROP(CCMpfaModel, ElementVolumeVariables, Dumux::CCMpfaElementVolumeVar //! The local flux variables cache vector class SET_TYPE_PROP(CCMpfaModel, ElementFluxVariablesCache, Dumux::CCMpfaElementFluxVariablesCache); -//! The stencil container, some mpfa methods might lead to non-symmetric global systems -SET_TYPE_PROP(CCMpfaModel, StencilsVector, Dumux::CCMpfaStencilsVector); - SET_PROP(CCMpfaModel, SubControlVolumeFace) { private: diff --git a/dumux/implicit/cellcentered/propertydefaults.hh b/dumux/implicit/cellcentered/propertydefaults.hh index 03c46d81af62e4120b660dc7ee0f07a5b2ac4e24..20cadb848c32b332e8e8a0f1914571b823d2f078 100644 --- a/dumux/implicit/cellcentered/propertydefaults.hh +++ b/dumux/implicit/cellcentered/propertydefaults.hh @@ -76,8 +76,8 @@ SET_TYPE_PROP(CCModel, GlobalVolumeVariables, CCGlobalVolumeVariables); -//! Set the AssemblyMap to the CCAssemblyMap -SET_TYPE_PROP(CCModel, AssemblyMap, Dumux::CCAssemblyMap); +//! Set the AssemblyMap to the CCAssemblyMap (default: symmetric occupation pattern) +SET_TYPE_PROP(CCModel, AssemblyMap, Dumux::CCSymmetricAssemblyMap); //! indicate that this is no box discretization SET_BOOL_PROP(CCModel, ImplicitIsBox, false); diff --git a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh index f029f862d61f10951383a2f0b29437932c84169b..1422045377a834e5b1d48d561026c5d65093a619 100644 --- a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh @@ -35,7 +35,6 @@ #include #include #include -#include #include #include @@ -66,9 +65,6 @@ SET_TYPE_PROP(CCTpfaModel, ElementVolumeVariables, Dumux::CCTpfaElementVolumeVar //! The local flux variables cache vector class SET_TYPE_PROP(CCTpfaModel, ElementFluxVariablesCache, Dumux::CCTpfaElementFluxVariablesCache); -//! The stencil container, tpfa leads to a symmetric global matrix -SET_TYPE_PROP(CCTpfaModel, StencilsVector, Dumux::CCSymmetricStencilsVector); - SET_PROP(CCTpfaModel, SubControlVolumeFace) { private: diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index 5d846855c31c41298a829fa56e38e367cb5f9d8e..a067d55a698951c230053405bb4eccfca9c7f7be 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -71,7 +71,6 @@ class ImplicitModel using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); using GlobalFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GlobalFluxVariablesCache); - using StencilsVector = typename GET_PROP_TYPE(TypeTag, StencilsVector); enum { numEq = GET_PROP_VALUE(TypeTag, NumEq), @@ -129,9 +128,6 @@ public: // resize and update the volVars with the initial solution curGlobalVolVars_.update(problem, curSol()); - // update stencils - stencilsVector_.update(problem); - // initialize assembler and create matrix localJacobian_.init(problem); jacAsm_ = std::make_shared(); @@ -845,23 +841,12 @@ protected: // the flux variables cache vector vector GlobalFluxVariablesCache globalfluxVarsCache_; - // the stencils vector - StencilsVector stencilsVector_; - // the finite volume element geometries std::shared_ptr globalFvGeometryPtr_; // container to store the box volumes Dune::BlockVector > boxVolume_; -public: - //! Get stencils related to an entity of specified codimension - // Cell-centered discretizations typically only implement element stencils - // Vertex-centered discretizations both vertex and element stencils - template - auto stencils(const Entity& entity) const -> decltype(stencilsVector_.get(entity)) - { return stencilsVector_.get(entity); } - private: /*! * \brief Returns whether messages should be printed diff --git a/dumux/implicit/properties.hh b/dumux/implicit/properties.hh index d081d1503046f6cf096840ee0a11e76c411a1fec..44b5c023c88d5e41c96e4825b94c4d448cc702ad 100644 --- a/dumux/implicit/properties.hh +++ b/dumux/implicit/properties.hh @@ -103,9 +103,6 @@ NEW_PROP_TAG(SolutionDependentHeatConduction); //!< specifies if the parameters NEW_PROP_TAG(VtkAddVelocity); //!< specifies if an element velocity it reconstructed for the output NEW_PROP_TAG(VtkAddProcessRank); //!< specifies if the process rank should be added the output -// stencils -NEW_PROP_TAG(StencilsVector); //! The type of the global vector of stencils per element - // high level simulation control NEW_PROP_TAG(TimeManager); //!< Manages the simulation time NEW_PROP_TAG(NewtonMethod); //!< The type of the newton method diff --git a/dumux/porousmediumflow/implicit/velocityoutput.hh b/dumux/porousmediumflow/implicit/velocityoutput.hh index b9539a510ac8101d7b84d2769ce439361adbdaa0..3e75f61866e4e0e0bf4eee0fa58f56c0c0c4867c 100644 --- a/dumux/porousmediumflow/implicit/velocityoutput.hh +++ b/dumux/porousmediumflow/implicit/velocityoutput.hh @@ -65,7 +65,6 @@ class ImplicitVelocityOutput using Element = typename GridView::template Codim<0>::Entity; using IndexType = typename GridView::IndexSet::IndexType; using CoordScalar = typename GridView::ctype; - using Stencil = std::vector; using GlobalPosition = Dune::FieldVector; using ReferenceElements = Dune::ReferenceElements; @@ -89,8 +88,9 @@ public: // resize to the number of vertices of the grid cellNum_.assign(problem.gridView().size(dim), 0); - for (const auto& vertex : vertices(problem.gridView())) - cellNum_[problem.vertexMapper().index(vertex)] = getStencil(vertex).size(); + for (const auto& element : elements(problem.gridView())) + for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx) + ++cellNum_[problem.vertexMapper().subIndex(element, vIdx, dim)]; } } } @@ -111,8 +111,7 @@ public: // compiler error. However, that call is needed for calculating velocities // if the cell-centered discretization is used. By proceeding as in the // following lines, that call will only be compiled if cell-centered - // actually is used. For the same reason we also provide a isBox-specific - // implementation of the getStencil method below. + // actually is used. template typename std::enable_if::type problemBoundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const @@ -124,18 +123,6 @@ public: problemBoundaryTypes(const Element& element, const SubControlVolume& scv) const { return BoundaryTypes(); } - //! returns the elements connected to a vertex - template - const typename std::enable_if::type& - getStencil(const Vertex& vertex) const - { return problem_.model().stencils(vertex).elementIndices(); } - - //! we should never call this method for cc models - template - const typename std::enable_if::type - getStencil(const Vertex& vertex) const - { return Stencil(0); } - //! Calculate the velocities for the scvs in the element //! We assume the local containers to be bound to the complete stencil template