From df2b6156cb12760b199317581ef3c4af47cd1bc1 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch <bernd@iws.uni-stuttgart.de> Date: Wed, 20 Jan 2016 12:27:49 +0100 Subject: [PATCH] [folder structure] rename files from implicit/cellcentered, adapt includes --- dumux/implicit/1p/1pproperties.hh | 2 +- dumux/implicit/1p2c/1p2cproperties.hh | 2 +- dumux/implicit/2p/2pproperties.hh | 2 +- dumux/implicit/2p2c/2p2cproperties.hh | 2 +- dumux/implicit/2pdfm/2pdfmproperties.hh | 2 +- dumux/implicit/2pnc/2pncproperties.hh | 2 +- dumux/implicit/3p/3pproperties.hh | 2 +- dumux/implicit/3p3c/3p3cproperties.hh | 2 +- dumux/implicit/cellcentered/assembler.hh | 253 +++++++++++++++ dumux/implicit/cellcentered/ccassembler.hh | 253 +-------------- .../cellcentered/ccelementboundarytypes.hh | 140 +-------- .../cellcentered/ccelementvolumevariables.hh | 142 +-------- .../cellcentered/ccfvelementgeometry.hh | 235 +------------- .../implicit/cellcentered/cclocalresidual.hh | 288 +----------------- dumux/implicit/cellcentered/ccproperties.hh | 52 +--- .../cellcentered/ccpropertydefaults.hh | 72 +---- .../cellcentered/elementboundarytypes.hh | 140 +++++++++ .../cellcentered/elementvolumevariables.hh | 142 +++++++++ .../cellcentered/fvelementgeometry.hh | 235 ++++++++++++++ dumux/implicit/cellcentered/localresidual.hh | 288 ++++++++++++++++++ dumux/implicit/cellcentered/properties.hh | 52 ++++ .../implicit/cellcentered/propertydefaults.hh | 72 +++++ .../cornerpoint/cpelementvolumevariables.hh | 2 +- dumux/implicit/mpnc/mpncproperties.hh | 2 +- dumux/implicit/richards/richardsproperties.hh | 2 +- dumux/linear/amgparallelhelpers.hh | 2 +- dumux/linear/amgproperties.hh | 2 +- .../2p/implicit/cc2pcornerpointproblem.hh | 2 +- .../2p/implicit/lensproblem.hh | 2 +- 29 files changed, 1225 insertions(+), 1169 deletions(-) create mode 100644 dumux/implicit/cellcentered/assembler.hh create mode 100644 dumux/implicit/cellcentered/elementboundarytypes.hh create mode 100644 dumux/implicit/cellcentered/elementvolumevariables.hh create mode 100644 dumux/implicit/cellcentered/fvelementgeometry.hh create mode 100644 dumux/implicit/cellcentered/localresidual.hh create mode 100644 dumux/implicit/cellcentered/properties.hh create mode 100644 dumux/implicit/cellcentered/propertydefaults.hh diff --git a/dumux/implicit/1p/1pproperties.hh b/dumux/implicit/1p/1pproperties.hh index e511d700cf..a563291924 100644 --- a/dumux/implicit/1p/1pproperties.hh +++ b/dumux/implicit/1p/1pproperties.hh @@ -28,7 +28,7 @@ #define DUMUX_1P_PROPERTIES_HH #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> #include <dumux/implicit/nonisothermal/niproperties.hh> namespace Dumux diff --git a/dumux/implicit/1p2c/1p2cproperties.hh b/dumux/implicit/1p2c/1p2cproperties.hh index 52bc7eaf9e..64cde96026 100644 --- a/dumux/implicit/1p2c/1p2cproperties.hh +++ b/dumux/implicit/1p2c/1p2cproperties.hh @@ -31,7 +31,7 @@ #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> #include <dumux/implicit/nonisothermal/niproperties.hh> namespace Dumux diff --git a/dumux/implicit/2p/2pproperties.hh b/dumux/implicit/2p/2pproperties.hh index ccc78dd982..c746920ba1 100644 --- a/dumux/implicit/2p/2pproperties.hh +++ b/dumux/implicit/2p/2pproperties.hh @@ -31,7 +31,7 @@ #define DUMUX_2P_PROPERTIES_HH #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> #include <dumux/implicit/nonisothermal/niproperties.hh> namespace Dumux diff --git a/dumux/implicit/2p2c/2p2cproperties.hh b/dumux/implicit/2p2c/2p2cproperties.hh index d1bbee33a0..64a9f1e9be 100644 --- a/dumux/implicit/2p2c/2p2cproperties.hh +++ b/dumux/implicit/2p2c/2p2cproperties.hh @@ -30,7 +30,7 @@ #define DUMUX_2P2C_PROPERTIES_HH #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> #include<dumux/implicit/nonisothermal/niproperties.hh> namespace Dumux diff --git a/dumux/implicit/2pdfm/2pdfmproperties.hh b/dumux/implicit/2pdfm/2pdfmproperties.hh index b95724ca86..839e49fa46 100644 --- a/dumux/implicit/2pdfm/2pdfmproperties.hh +++ b/dumux/implicit/2pdfm/2pdfmproperties.hh @@ -28,7 +28,7 @@ #define DUMUX_MODELS_2PDFM_PROPERTIES_HH #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> namespace Dumux { diff --git a/dumux/implicit/2pnc/2pncproperties.hh b/dumux/implicit/2pnc/2pncproperties.hh index e62eeb4f4b..021646ff55 100644 --- a/dumux/implicit/2pnc/2pncproperties.hh +++ b/dumux/implicit/2pnc/2pncproperties.hh @@ -30,7 +30,7 @@ #define DUMUX_2PNC_PROPERTIES_HH #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> #include <dumux/implicit/nonisothermal/niproperties.hh> namespace Dumux diff --git a/dumux/implicit/3p/3pproperties.hh b/dumux/implicit/3p/3pproperties.hh index bb1619d988..aaa53a3491 100644 --- a/dumux/implicit/3p/3pproperties.hh +++ b/dumux/implicit/3p/3pproperties.hh @@ -28,7 +28,7 @@ #define DUMUX_3P_PROPERTIES_HH #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> #include <dumux/implicit/nonisothermal/niproperties.hh> namespace Dumux diff --git a/dumux/implicit/3p3c/3p3cproperties.hh b/dumux/implicit/3p3c/3p3cproperties.hh index 7d8fb1884b..98104b81c6 100644 --- a/dumux/implicit/3p3c/3p3cproperties.hh +++ b/dumux/implicit/3p3c/3p3cproperties.hh @@ -31,7 +31,7 @@ #define DUMUX_3P3C_PROPERTIES_HH #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> #include <dumux/implicit/nonisothermal/niproperties.hh> diff --git a/dumux/implicit/cellcentered/assembler.hh b/dumux/implicit/cellcentered/assembler.hh new file mode 100644 index 0000000000..bdfcc15295 --- /dev/null +++ b/dumux/implicit/cellcentered/assembler.hh @@ -0,0 +1,253 @@ +// -*- 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 + * \brief An assembler for the global Jacobian matrix for models using the cell centered discretization. + */ +#ifndef DUMUX_CC_ASSEMBLER_HH +#define DUMUX_CC_ASSEMBLER_HH + +#include <dumux/implicit/assembler.hh> + +namespace Dumux { + +/*! + * \ingroup CCModel + * \brief An assembler for the global Jacobian matrix for models using the cell centered discretization. + */ +template<class TypeTag> +class CCAssembler : public ImplicitAssembler<TypeTag> +{ + typedef ImplicitAssembler<TypeTag> ParentType; + friend class ImplicitAssembler<TypeTag>; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix; + typedef typename GridView::template Codim<0>::Entity Element; + +public: + CCAssembler(): ParentType() {} + +private: + // copying the jacobian assembler is not a good idea + CCAssembler(const CCAssembler &); + + /*! + * \brief Determine the colors of elements for partial + * reassembly given a relative tolerance. + * + * The following approach is used: + * + * - Set all elements to 'green' + * - Mark all elements as 'red' which exhibit an relative error above + * the tolerance + * - Mark all neighbors of 'red' elements also 'red' + * + * \param relTol The relative error below which an element won't be + * reassembled. Note that this specifies the + * worst-case relative error between the last + * linearization point and the current solution and + * _not_ the delta vector of the Newton iteration! + */ + void computeColors_(Scalar relTol) + { + if (!this->enablePartialReassemble_()) + return; + + // mark the red elements and update the tolerance of the + // linearization which actually will get achieved + this->nextReassembleAccuracy_ = 0; + for (const auto& element : Dune::elements(this->gridView_())) { + int eIdx = this->elementMapper_().index(element); + + if (this->delta_[eIdx] > relTol) + { + // mark element as red if discrepancy is larger than + // the relative tolerance + this->elementColor_[eIdx] = ParentType::Red; + } + else + { + this->elementColor_[eIdx] = ParentType::Green; + this->nextReassembleAccuracy_ = + std::max(this->nextReassembleAccuracy_, this->delta_[eIdx]); + } + } + + // mark the neighbors also red + for (const auto& element : Dune::elements(this->gridView_())) { + int eIdx = this->elementMapper_().index(element); + + if (this->elementColor_[eIdx] == ParentType::Red) + continue; // element is red already! + + if (this->delta_[eIdx] > relTol) + { + // also mark the neighbors + for (const auto& intersection : Dune::intersections(this->gridView_(), element)) + { + if (intersection.neighbor()) + { + int neighborIdx = this->elementMapper_().index(intersection.outside()); + + this->elementColor_[neighborIdx] = ParentType::Red; + } + } + } + } + + // set the discrepancy of the red elements to zero + for (unsigned int i = 0; i < this->delta_.size(); i++) + if (this->elementColor_[i] == ParentType::Red) + this->delta_[i] = 0; + } + + // Construct the BCRS matrix for the global jacobian + void createMatrix_() + { + int numElements = this->gridView_().size(0); + + // allocate raw matrix + this->matrix_ = std::make_shared<JacobianMatrix>(numElements, numElements, JacobianMatrix::random); + + // find out the global indices of the neighboring elements of + // each element + typedef std::set<int> NeighborSet; + std::vector<NeighborSet> neighbors(numElements); + for (const auto& element : Dune::elements(this->gridView_())) { + + int globalI = this->elementMapper_().index(element); + + neighbors[globalI].insert(globalI); + + // if the element is ghost, + // all dofs just contain main-diagonal entries + //if (element.partitionType() == Dune::GhostEntity) + // continue; + + // loop over all neighbors + for (const auto& intersection : Dune::intersections(this->gridView_(), element)) + { + if (intersection.neighbor()) + { + int globalJ = this->elementMapper_().index(intersection.outside()); + + neighbors[globalI].insert(globalJ); + } + } + } + + // allocate space for the rows of the matrix + for (int i = 0; i < numElements; ++i) { + this->matrix_->setrowsize(i, neighbors[i].size()); + } + this->matrix_->endrowsizes(); + + // fill the rows with indices. each element talks to all of its + // neighbors and itself. + for (int i = 0; i < numElements; ++i) { + typename NeighborSet::iterator nIt = neighbors[i].begin(); + typename NeighborSet::iterator nEndIt = neighbors[i].end(); + for (; nIt != nEndIt; ++nIt) { + this->matrix_->addindex(i, *nIt); + } + } + this->matrix_->endindices(); + } + + // assemble a non-ghost element + void assembleElement_(const Element &element) + { + if (this->enablePartialReassemble_()) { + int eIdxGlobal = this->model_().elementMapper().index(element); + + if (this->elementColor_[eIdxGlobal] == ParentType::Green) { + ++this->greenElems_; + + assembleGreenElement_(element); + return; + } + } + + this->model_().localJacobian().assemble(element); + + int globalI = this->elementMapper_().index(element); + + + // update the right hand side + this->residual_[globalI] = this->model_().localJacobian().residual(0); + for (int j = 0; j < this->residual_[globalI].dimension; ++j) + assert(std::isfinite(this->residual_[globalI][j])); + if (this->enableJacobianRecycling_()) { + this->storageTerm_[globalI] += + this->model_().localJacobian().storageTerm(0); + } + + if (this->enableJacobianRecycling_()) + this->storageJacobian_[globalI] += + this->model_().localJacobian().storageJacobian(0); + + // update the diagonal entry + (*this->matrix_)[globalI][globalI] = this->model_().localJacobian().mat(0,0); + + int j = 0; + for (const auto& intersection : Dune::intersections(this->gridView_(), element)) + { + if (intersection.neighbor()) + { + int globalJ = this->elementMapper_().index(intersection.outside()); + + (*this->matrix_)[globalI][globalJ] = this->model_().localJacobian().mat(0,++j); + } + } + } + + // "assemble" a green element. green elements only get the + // residual updated, but the jacobian is left alone... + void assembleGreenElement_(const Element &element) + { + this->model_().localResidual().eval(element); + + int globalI = this->elementMapper_().index(element); + + // update the right hand side + this->residual_[globalI] += this->model_().localResidual().residual(0); + if (this->enableJacobianRecycling_()) + this->storageTerm_[globalI] += this->model_().localResidual().storageTerm(0); + } + + // "assemble" a ghost element + void assembleGhostElement_(const Element &element) + { + int globalI = this->elementMapper_().index(element); + + // update the right hand side + this->residual_[globalI] = 0.0; + + // update the diagonal entry + typedef typename JacobianMatrix::block_type BlockType; + BlockType &J = (*this->matrix_)[globalI][globalI]; + for (int j = 0; j < BlockType::rows; ++j) + J[j][j] = 1.0; + } +}; + +} // namespace Dumux + +#endif diff --git a/dumux/implicit/cellcentered/ccassembler.hh b/dumux/implicit/cellcentered/ccassembler.hh index bdfcc15295..805132a908 100644 --- a/dumux/implicit/cellcentered/ccassembler.hh +++ b/dumux/implicit/cellcentered/ccassembler.hh @@ -1,253 +1,8 @@ -// -*- 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 - * \brief An assembler for the global Jacobian matrix for models using the cell centered discretization. - */ -#ifndef DUMUX_CC_ASSEMBLER_HH -#define DUMUX_CC_ASSEMBLER_HH +#ifndef DUMUX_CC_ASSEMBLER_HH_OLD +#define DUMUX_CC_ASSEMBLER_HH_OLD -#include <dumux/implicit/assembler.hh> +#warning this header is deprecated, use dumux/implicit/cellcentered/assembler.hh instead -namespace Dumux { - -/*! - * \ingroup CCModel - * \brief An assembler for the global Jacobian matrix for models using the cell centered discretization. - */ -template<class TypeTag> -class CCAssembler : public ImplicitAssembler<TypeTag> -{ - typedef ImplicitAssembler<TypeTag> ParentType; - friend class ImplicitAssembler<TypeTag>; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix; - typedef typename GridView::template Codim<0>::Entity Element; - -public: - CCAssembler(): ParentType() {} - -private: - // copying the jacobian assembler is not a good idea - CCAssembler(const CCAssembler &); - - /*! - * \brief Determine the colors of elements for partial - * reassembly given a relative tolerance. - * - * The following approach is used: - * - * - Set all elements to 'green' - * - Mark all elements as 'red' which exhibit an relative error above - * the tolerance - * - Mark all neighbors of 'red' elements also 'red' - * - * \param relTol The relative error below which an element won't be - * reassembled. Note that this specifies the - * worst-case relative error between the last - * linearization point and the current solution and - * _not_ the delta vector of the Newton iteration! - */ - void computeColors_(Scalar relTol) - { - if (!this->enablePartialReassemble_()) - return; - - // mark the red elements and update the tolerance of the - // linearization which actually will get achieved - this->nextReassembleAccuracy_ = 0; - for (const auto& element : Dune::elements(this->gridView_())) { - int eIdx = this->elementMapper_().index(element); - - if (this->delta_[eIdx] > relTol) - { - // mark element as red if discrepancy is larger than - // the relative tolerance - this->elementColor_[eIdx] = ParentType::Red; - } - else - { - this->elementColor_[eIdx] = ParentType::Green; - this->nextReassembleAccuracy_ = - std::max(this->nextReassembleAccuracy_, this->delta_[eIdx]); - } - } - - // mark the neighbors also red - for (const auto& element : Dune::elements(this->gridView_())) { - int eIdx = this->elementMapper_().index(element); - - if (this->elementColor_[eIdx] == ParentType::Red) - continue; // element is red already! - - if (this->delta_[eIdx] > relTol) - { - // also mark the neighbors - for (const auto& intersection : Dune::intersections(this->gridView_(), element)) - { - if (intersection.neighbor()) - { - int neighborIdx = this->elementMapper_().index(intersection.outside()); - - this->elementColor_[neighborIdx] = ParentType::Red; - } - } - } - } - - // set the discrepancy of the red elements to zero - for (unsigned int i = 0; i < this->delta_.size(); i++) - if (this->elementColor_[i] == ParentType::Red) - this->delta_[i] = 0; - } - - // Construct the BCRS matrix for the global jacobian - void createMatrix_() - { - int numElements = this->gridView_().size(0); - - // allocate raw matrix - this->matrix_ = std::make_shared<JacobianMatrix>(numElements, numElements, JacobianMatrix::random); - - // find out the global indices of the neighboring elements of - // each element - typedef std::set<int> NeighborSet; - std::vector<NeighborSet> neighbors(numElements); - for (const auto& element : Dune::elements(this->gridView_())) { - - int globalI = this->elementMapper_().index(element); - - neighbors[globalI].insert(globalI); - - // if the element is ghost, - // all dofs just contain main-diagonal entries - //if (element.partitionType() == Dune::GhostEntity) - // continue; - - // loop over all neighbors - for (const auto& intersection : Dune::intersections(this->gridView_(), element)) - { - if (intersection.neighbor()) - { - int globalJ = this->elementMapper_().index(intersection.outside()); - - neighbors[globalI].insert(globalJ); - } - } - } - - // allocate space for the rows of the matrix - for (int i = 0; i < numElements; ++i) { - this->matrix_->setrowsize(i, neighbors[i].size()); - } - this->matrix_->endrowsizes(); - - // fill the rows with indices. each element talks to all of its - // neighbors and itself. - for (int i = 0; i < numElements; ++i) { - typename NeighborSet::iterator nIt = neighbors[i].begin(); - typename NeighborSet::iterator nEndIt = neighbors[i].end(); - for (; nIt != nEndIt; ++nIt) { - this->matrix_->addindex(i, *nIt); - } - } - this->matrix_->endindices(); - } - - // assemble a non-ghost element - void assembleElement_(const Element &element) - { - if (this->enablePartialReassemble_()) { - int eIdxGlobal = this->model_().elementMapper().index(element); - - if (this->elementColor_[eIdxGlobal] == ParentType::Green) { - ++this->greenElems_; - - assembleGreenElement_(element); - return; - } - } - - this->model_().localJacobian().assemble(element); - - int globalI = this->elementMapper_().index(element); - - - // update the right hand side - this->residual_[globalI] = this->model_().localJacobian().residual(0); - for (int j = 0; j < this->residual_[globalI].dimension; ++j) - assert(std::isfinite(this->residual_[globalI][j])); - if (this->enableJacobianRecycling_()) { - this->storageTerm_[globalI] += - this->model_().localJacobian().storageTerm(0); - } - - if (this->enableJacobianRecycling_()) - this->storageJacobian_[globalI] += - this->model_().localJacobian().storageJacobian(0); - - // update the diagonal entry - (*this->matrix_)[globalI][globalI] = this->model_().localJacobian().mat(0,0); - - int j = 0; - for (const auto& intersection : Dune::intersections(this->gridView_(), element)) - { - if (intersection.neighbor()) - { - int globalJ = this->elementMapper_().index(intersection.outside()); - - (*this->matrix_)[globalI][globalJ] = this->model_().localJacobian().mat(0,++j); - } - } - } - - // "assemble" a green element. green elements only get the - // residual updated, but the jacobian is left alone... - void assembleGreenElement_(const Element &element) - { - this->model_().localResidual().eval(element); - - int globalI = this->elementMapper_().index(element); - - // update the right hand side - this->residual_[globalI] += this->model_().localResidual().residual(0); - if (this->enableJacobianRecycling_()) - this->storageTerm_[globalI] += this->model_().localResidual().storageTerm(0); - } - - // "assemble" a ghost element - void assembleGhostElement_(const Element &element) - { - int globalI = this->elementMapper_().index(element); - - // update the right hand side - this->residual_[globalI] = 0.0; - - // update the diagonal entry - typedef typename JacobianMatrix::block_type BlockType; - BlockType &J = (*this->matrix_)[globalI][globalI]; - for (int j = 0; j < BlockType::rows; ++j) - J[j][j] = 1.0; - } -}; - -} // namespace Dumux +#include <dumux/implicit/cellcentered/assembler.hh> #endif diff --git a/dumux/implicit/cellcentered/ccelementboundarytypes.hh b/dumux/implicit/cellcentered/ccelementboundarytypes.hh index ae33fdffdf..9d304de6f6 100644 --- a/dumux/implicit/cellcentered/ccelementboundarytypes.hh +++ b/dumux/implicit/cellcentered/ccelementboundarytypes.hh @@ -1,140 +1,8 @@ -// -*- 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 - * \brief Boundary types gathered on an element - */ -#ifndef DUMUX_CC_ELEMENT_BOUNDARY_TYPES_HH -#define DUMUX_CC_ELEMENT_BOUNDARY_TYPES_HH +#ifndef DUMUX_CC_ELEMENT_BOUNDARY_TYPES_HH_OLD +#define DUMUX_CC_ELEMENT_BOUNDARY_TYPES_HH_OLD -#include "ccproperties.hh" +#warning this header is deprecated, use dumux/implicit/cellcentered/elementboundarytypes.hh instead -#include <dumux/common/valgrind.hh> - -namespace Dumux -{ - -/*! - * \ingroup CCModel - * \ingroup ImplicitBoundaryTypes - * \brief This class stores an array of BoundaryTypes objects - */ -template<class TypeTag> -class CCElementBoundaryTypes : public std::vector<typename GET_PROP_TYPE(TypeTag, BoundaryTypes) > -{ - typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; - typedef std::vector<BoundaryTypes> ParentType; - - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; - - typedef typename GridView::template Codim<0>::Entity Element; - -public: - /*! - * \brief Copy constructor. - * - * Copying a the boundary types of an element should be explicitly - * requested - */ - explicit CCElementBoundaryTypes(const CCElementBoundaryTypes &v) - : ParentType(v) - {} - - /*! - * \brief Default constructor. - */ - CCElementBoundaryTypes() - { - hasDirichlet_ = false; - hasNeumann_ = false; - hasOutflow_ = false; - } - - /*! - * \brief Update the boundary types for all vertices of an element. - * - * \param problem The problem object which needs to be simulated - * \param element The DUNE Codim<0> entity for which the boundary - * types should be collected - */ - void update(const Problem &problem, - const Element &element) - { - this->resize(1); - - hasDirichlet_ = false; - hasNeumann_ = false; - hasOutflow_ = false; - - (*this)[0].reset(); - - if (!problem.model().onBoundary(element)) - return; - - for (const auto& intersection : Dune::intersections(problem.gridView(), element)) { - if (!intersection.boundary()) - continue; - - problem.boundaryTypes((*this)[0], intersection); - - hasDirichlet_ = hasDirichlet_ || (*this)[0].hasDirichlet(); - hasNeumann_ = hasNeumann_ || (*this)[0].hasNeumann(); - hasOutflow_ = hasOutflow_ || (*this)[0].hasOutflow(); - } - } - - void update(const Problem &problem, - const Element &element, - const FVElementGeometry &fvGeometry) - { - update(problem, element); - } - - /*! - * \brief Returns whether the element has a vertex which contains - * a Dirichlet value. - */ - bool hasDirichlet() const - { return hasDirichlet_; } - - /*! - * \brief Returns whether the element potentially features a - * Neumann boundary segment. - */ - bool hasNeumann() const - { return hasNeumann_; } - - /*! - * \brief Returns whether the element potentially features an - * outflow boundary segment. - */ - bool hasOutflow() const - { return hasOutflow_; } - -protected: - bool hasDirichlet_; - bool hasNeumann_; - bool hasOutflow_; -}; - -} // namespace Dumux +#include <dumux/implicit/cellcentered/elementboundarytypes.hh> #endif diff --git a/dumux/implicit/cellcentered/ccelementvolumevariables.hh b/dumux/implicit/cellcentered/ccelementvolumevariables.hh index 13f28c42d4..a7258d246c 100644 --- a/dumux/implicit/cellcentered/ccelementvolumevariables.hh +++ b/dumux/implicit/cellcentered/ccelementvolumevariables.hh @@ -1,142 +1,8 @@ -// -*- 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 - * \brief Volume variables gathered on an element - */ -#ifndef DUMUX_CC_ELEMENT_VOLUME_VARIABLES_HH -#define DUMUX_CC_ELEMENT_VOLUME_VARIABLES_HH +#ifndef DUMUX_CC_ELEMENT_VOLUME_VARIABLES_HH_OLD +#define DUMUX_CC_ELEMENT_VOLUME_VARIABLES_HH_OLD -#include "ccproperties.hh" +#warning this header is deprecated, use dumux/implicit/cellcentered/elementvolumevariables.hh instead -namespace Dumux -{ - -/*! - * \ingroup CCModel - * \brief This class stores an array of VolumeVariables objects, one - * volume variables object for each of the element's vertices - */ -template<class TypeTag> -class CCElementVolumeVariables : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables) > -{ - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; - typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; - typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; - typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; - - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GridView::template Codim<0>::Entity Element; - -public: - /*! - * \brief The constructor. - */ - CCElementVolumeVariables() - { } - - /*! - * \brief Construct the volume variables for all of vertices of an element. - * - * \param problem The problem which needs to be simulated. - * \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated - * \param fvGeometry The finite volume geometry of the element - * \param oldSol Tells whether the model's previous or current solution should be used. - */ - void update(const Problem &problem, - const Element &element, - const FVElementGeometry &fvGeometry, - bool oldSol) - { - const SolutionVector &globalSol = - oldSol? - problem.model().prevSol(): - problem.model().curSol(); - - int numNeighbors = fvGeometry.numNeighbors; - this->resize(numNeighbors); - - for (int i = 0; i < numNeighbors; i++) - { - const Element& neighbor = fvGeometry.neighbors[i]; - - const PrimaryVariables &solI - = globalSol[problem.elementMapper().index(neighbor)]; - - FVElementGeometry neighborFVGeom; - neighborFVGeom.updateInner(neighbor); - - (*this)[i].update(solI, - problem, - neighbor, - neighborFVGeom, - /*scvIdx=*/0, - oldSol); - } - - // only treat boundary if current solution is evaluated - if (!oldSol) - { - // check if element intersects with the boundary - ElementBoundaryTypes elemBCTypes; - elemBCTypes.update(problem, element); - if (elemBCTypes.hasDirichlet() - || elemBCTypes.hasNeumann() - || elemBCTypes.hasOutflow()) - { - this->resize(numNeighbors + element.subEntities(1)); - - // add volume variables for the boundary faces - for (const auto& intersection : Dune::intersections(problem.gridView(), element)) { - if (!intersection.boundary()) - continue; - - BoundaryTypes bcTypes; - problem.boundaryTypes(bcTypes, intersection); - - int fIdx = intersection.indexInInside(); - int indexInVariables = numNeighbors + fIdx; - - if (bcTypes.hasDirichlet()) - { - PrimaryVariables dirichletValues; - problem.dirichlet(dirichletValues, intersection); - - (*this)[indexInVariables].update(dirichletValues, - problem, - element, - fvGeometry, - /*scvIdx=*/0, - oldSol); - } - else - { - (*this)[indexInVariables] = (*this)[0]; - } - } - } - } - } -}; - -} // namespace Dumux +#include <dumux/implicit/cellcentered/elementvolumevariables.hh> #endif diff --git a/dumux/implicit/cellcentered/ccfvelementgeometry.hh b/dumux/implicit/cellcentered/ccfvelementgeometry.hh index 0603825c95..dad20f4340 100644 --- a/dumux/implicit/cellcentered/ccfvelementgeometry.hh +++ b/dumux/implicit/cellcentered/ccfvelementgeometry.hh @@ -1,235 +1,8 @@ -// -*- 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 - * \brief Represents the finite volume geometry of a single element in - * the cell-centered fv scheme. - */ -#ifndef DUMUX_CC_FV_ELEMENTGEOMETRY_HH -#define DUMUX_CC_FV_ELEMENTGEOMETRY_HH +#ifndef DUMUX_CC_FV_ELEMENTGEOMETRY_HH_OLD +#define DUMUX_CC_FV_ELEMENTGEOMETRY_HH_OLD -#include <dune/geometry/affinegeometry.hh> -#include <dune/geometry/referenceelements.hh> -#include <dune/grid/common/intersectioniterator.hh> +#warning this header is deprecated, use dumux/implicit/cellcentered/fvelementgeometry.hh instead -#include <dumux/common/propertysystem.hh> - -namespace Dumux -{ -namespace Properties -{ -NEW_PROP_TAG(GridView); -NEW_PROP_TAG(Scalar); -} - -/*! - * \ingroup CCModel - * \brief Represents the finite volume geometry of a single element in - * the cell-centered fv scheme. - */ -template<class TypeTag> -class CCFVElementGeometry -{ - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - enum{dim = GridView::dimension}; - enum{dimWorld = GridView::dimensionworld}; - - enum{maxNFAP = 2}; //! maximum number of flux approximation points (two-point flux) - enum{maxNE = 2*dim*(1 << (dim - 1))}; //! maximum number of neighbors (works for one hanging node per face) - enum{maxBF = 2*dim}; //! maximum number of boundary faces - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GridView::ctype CoordScalar; - typedef typename GridView::Traits::template Codim<0>::Entity Element; - typedef typename Element::Geometry Geometry; - typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition; - typedef Dune::FieldVector<CoordScalar,dim> LocalPosition; - typedef typename Dune::AffineGeometry<CoordScalar, dim, dimWorld> AffineGeometry; - -public: - struct SubControlVolume //! FV intersected with element - { - LocalPosition local; //!< local position - GlobalPosition global; //!< global position - Scalar volume; //!< volume of scv - bool inner; - }; - - struct SubControlVolumeFace //! interior face of a sub control volume - { - int i,j; //!< scvf seperates corner i and j of elem - LocalPosition ipLocal; //!< integration point in local coords - GlobalPosition ipGlobal; //!< integration point in global coords - GlobalPosition normal; //!< normal on face pointing to CV j or outward of the domain with length equal to |scvf| - Scalar area; //!< area of face - Dune::FieldVector<GlobalPosition, maxNFAP> grad; //!< derivatives of shape functions at ip - Dune::FieldVector<Scalar, maxNFAP> shapeValue; //!< value of shape functions at ip - Dune::FieldVector<int, maxNFAP> fapIndices; //!< indices w.r.t.neighbors of the flux approximation points - unsigned numFap; //!< number of flux approximation points - unsigned fIdx; //!< the index (w.r.t. the element) of the face (codim 1 entity) that the scvf is part of - }; - - typedef SubControlVolumeFace BoundaryFace; //!< compatibility typedef - - LocalPosition elementLocal; //!< local coordinate of element center - GlobalPosition elementGlobal; //!< global coordinate of element center - Scalar elementVolume; //!< element volume - SubControlVolume subContVol[1]; //!< data of the sub control volumes - std::vector<AffineGeometry> subContVolGeometries; //!< geometries of the subcontrol volumes - SubControlVolumeFace subContVolFace[maxNE]; //!< data of the sub control volume faces - BoundaryFace boundaryFace[maxBF]; //!< data of the boundary faces - int numScv; //!< number of subcontrol volumes - int numScvf; //!< number of inner-domain subcontrolvolume faces - int numNeighbors; //!< number of neighboring elements including the element itself - std::vector<Element> neighbors; //!< stores the neighboring elements - - void updateInner(const Element& element) - { - const Geometry geometry = element.geometry(); - - elementVolume = geometry.volume(); - elementGlobal = geometry.center(); - elementLocal = geometry.local(elementGlobal); - - numScv = 1; - numScvf = 0; - - subContVol[0].local = elementLocal; - subContVol[0].global = elementGlobal; - subContVol[0].inner = true; - subContVol[0].volume = elementVolume; - - // initialize neighbors list with self: - numNeighbors = 1; - neighbors.clear(); - neighbors.reserve(maxNE); - neighbors.push_back(element); - } - - void update(const GridView& gridView, const Element& element) - { - updateInner(element); - - const Geometry geometry = element.geometry(); - - bool onBoundary = false; - - // fill neighbor information and control volume face data: - for (const auto& intersection : Dune::intersections(gridView, element)) - { - const auto isGeometry = intersection.geometry(); - - // neighbor information and inner cvf data: - if (intersection.neighbor()) - { - numNeighbors++; - neighbors.push_back(intersection.outside()); - - int scvfIdx = numNeighbors - 2; - SubControlVolumeFace& scvFace = subContVolFace[scvfIdx]; - - scvFace.i = 0; - scvFace.j = scvfIdx + 1; - - scvFace.ipGlobal = isGeometry.center(); - scvFace.ipLocal = geometry.local(scvFace.ipGlobal); - scvFace.normal = intersection.centerUnitOuterNormal(); - Scalar volume = isGeometry.volume(); - scvFace.normal *= volume; - scvFace.area = volume; - - GlobalPosition distVec = elementGlobal - - neighbors[scvfIdx+1].geometry().center(); - distVec /= distVec.two_norm2(); - - // gradients using a two-point flux approximation - scvFace.numFap = 2; - for (unsigned int fapIdx = 0; fapIdx < scvFace.numFap; fapIdx++) - { - scvFace.grad[fapIdx] = distVec; - scvFace.shapeValue[fapIdx] = 0.5; - } - scvFace.grad[1] *= -1.0; - - scvFace.fapIndices[0] = scvFace.i; - scvFace.fapIndices[1] = scvFace.j; - - scvFace.fIdx = intersection.indexInInside(); - } - - // boundary cvf data - if (intersection.boundary()) - { - onBoundary = true; - int bfIdx = intersection.indexInInside(); - SubControlVolumeFace& bFace = boundaryFace[bfIdx]; - - bFace.ipGlobal = isGeometry.center(); - bFace.ipLocal = geometry.local(bFace.ipGlobal); - bFace.normal = intersection.centerUnitOuterNormal(); - Scalar volume = isGeometry.volume(); - bFace.normal *= volume; - bFace.area = volume; - bFace.i = 0; - bFace.j = 0; - - GlobalPosition distVec = elementGlobal - bFace.ipGlobal; - distVec /= distVec.two_norm2(); - - // gradients using a two-point flux approximation - bFace.numFap = 2; - for (unsigned int fapIdx = 0; fapIdx < bFace.numFap; fapIdx++) - { - bFace.grad[fapIdx] = distVec; - bFace.shapeValue[fapIdx] = 0.5; - } - bFace.grad[1] *= -1.0; - - bFace.fapIndices[0] = bFace.i; - bFace.fapIndices[1] = bFace.j; - } - } - - // set the number of inner-domain subcontrolvolume faces - numScvf = numNeighbors - 1; - - // treat elements on the boundary - if (onBoundary) - { - for (int bfIdx = 0; bfIdx < element.subEntities(1); bfIdx++) - { - SubControlVolumeFace& bFace = boundaryFace[bfIdx]; - bFace.j = numNeighbors + bfIdx; - bFace.fapIndices[1] = bFace.j; - neighbors.push_back(element); - } - } - } - - /*! - * \brief For compatibilty with the box element geometry - */ - int boundaryFaceIndex(const int fIdx, const int vIdxInFace) const - { - return fIdx; - } -}; - -} +#include <dumux/implicit/cellcentered/fvelementgeometry.hh> #endif diff --git a/dumux/implicit/cellcentered/cclocalresidual.hh b/dumux/implicit/cellcentered/cclocalresidual.hh index cca821bc06..88d766295e 100644 --- a/dumux/implicit/cellcentered/cclocalresidual.hh +++ b/dumux/implicit/cellcentered/cclocalresidual.hh @@ -1,288 +1,8 @@ -// -*- 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 - * \brief Calculates the residual of models based on the box scheme element-wise. - */ -#ifndef DUMUX_CC_LOCAL_RESIDUAL_HH -#define DUMUX_CC_LOCAL_RESIDUAL_HH +#ifndef DUMUX_CC_LOCAL_RESIDUAL_HH_OLD +#define DUMUX_CC_LOCAL_RESIDUAL_HH_OLD -#include <dune/istl/matrix.hh> +#warning this header is deprecated, use dumux/implicit/cellcentered/localresidual.hh instead -#include <dumux/common/valgrind.hh> -#include <dumux/implicit/localresidual.hh> - -#include "ccproperties.hh" - -namespace Dumux -{ -/*! - * \ingroup CCModel - * \ingroup CCLocalResidual - * \brief Element-wise calculation of the residual for models - * based on the fully implicit cell-centered scheme. - * - * \todo Please doc me more! - */ -template<class TypeTag> -class CCLocalResidual : public ImplicitLocalResidual<TypeTag> -{ - typedef ImplicitLocalResidual<TypeTag> ParentType; - friend class ImplicitLocalResidual<TypeTag>; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - - enum { - numEq = GET_PROP_VALUE(TypeTag, NumEq) - }; - - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; - typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; - typedef typename FVElementGeometry::SubControlVolumeFace SCVFace; - - // copying the local residual class is not a good idea - CCLocalResidual(const CCLocalResidual &); - -public: - CCLocalResidual() : ParentType() - { } - -protected: - - /*! - * \brief Add all fluxes resulting from Neumann, outflow and pure Dirichlet - * boundary conditions to the local residual. - */ - void evalBoundaryFluxes_() - { - for (const auto& intersection : Dune::intersections(this->gridView_(), this->element_())) - { - // handle only faces on the boundary - if (!intersection.boundary()) - continue; - - BoundaryTypes bcTypes; - this->problem_().boundaryTypes(bcTypes, intersection); - - // evaluate the Neumann conditions at the boundary face - if (bcTypes.hasNeumann()) - this->asImp_().evalNeumannSegment_(&intersection, bcTypes); - - // evaluate the outflow conditions at the boundary face - if (bcTypes.hasOutflow()) - this->asImp_().evalOutflowSegment_(&intersection, bcTypes); - - // evaluate the pure Dirichlet conditions at the boundary face - if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann()) - this->asImp_().evalDirichletSegment_(&intersection, bcTypes); - } - } - - /*! - * \brief Evaluate Dirichlet conditions on faces that have mixed - * Dirichlet/Neumann boundary conditions. - */ - void evalDirichlet_() - { - for (const auto& intersection : Dune::intersections(this->gridView_(), this->element_())) - { - // handle only faces on the boundary - if (!intersection.boundary()) - continue; - - BoundaryTypes bcTypes; - this->problem_().boundaryTypes(bcTypes, intersection); - - if (bcTypes.hasDirichlet() && bcTypes.hasNeumann()) - this->asImp_().evalDirichletSegmentMixed_(&intersection, bcTypes); - } - } - - /*! - * \brief Add Neumann boundary conditions for a single intersection - */ - template <class IntersectionIterator> - void evalNeumannSegment_(const IntersectionIterator &isIt, - const BoundaryTypes &bcTypes) - { - // temporary vector to store the neumann boundary fluxes - PrimaryVariables values; - Valgrind::SetUndefined(values); - - unsigned bfIdx = isIt->indexInInside(); - this->problem_().solDependentNeumann(values, - this->element_(), - this->fvGeometry_(), - *isIt, - /*scvIdx=*/0, - bfIdx, - this->curVolVars_()); - values *= isIt->geometry().volume() - * this->curVolVars_(0).extrusionFactor(); - - // add fluxes to the residual - Valgrind::CheckDefined(values); - for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { - if (bcTypes.isNeumann(eqIdx)) - this->residual_[0][eqIdx] += values[eqIdx]; - } - } - - /*! - * \brief Add outflow boundary conditions for a single intersection - */ - template <class IntersectionIterator> - void evalOutflowSegment_(const IntersectionIterator &isIt, - const BoundaryTypes &bcTypes) - { - if (this->element_().geometry().type().isCube() == false) - DUNE_THROW(Dune::InvalidStateException, - "for cell-centered models, outflow BCs only work for cubes."); - - // store pointer to the current FVElementGeometry - const FVElementGeometry *oldFVGeometryPtr = this->fvElemGeomPtr_; - - // copy the current FVElementGeometry to a local variable - // and set the pointer to this local variable - FVElementGeometry fvGeometry = this->fvGeometry_(); - this->fvElemGeomPtr_ = &fvGeometry; - - // get the index of the boundary face - unsigned bfIdx = isIt->indexInInside(); - unsigned oppositeIdx = bfIdx^1; - - // manipulate the corresponding subcontrolvolume face - SCVFace& boundaryFace = fvGeometry.boundaryFace[bfIdx]; - - // set the second flux approximation index for the boundary face - for (int nIdx = 0; nIdx < fvGeometry.numNeighbors-1; nIdx++) - { - // check whether the two faces are opposite of each other - if (fvGeometry.subContVolFace[nIdx].fIdx == oppositeIdx) - { - boundaryFace.j = nIdx+1; - break; - } - } - boundaryFace.fapIndices[1] = boundaryFace.j; - boundaryFace.grad[0] *= -0.5; - boundaryFace.grad[1] *= -0.5; - - // temporary vector to store the outflow boundary fluxes - PrimaryVariables values; - Valgrind::SetUndefined(values); - - this->asImp_().computeFlux(values, bfIdx, true); - values *= this->curVolVars_(0).extrusionFactor(); - - // add fluxes to the residual - Valgrind::CheckDefined(values); - for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) - { - if (bcTypes.isOutflow(eqIdx)) - this->residual_[0][eqIdx] += values[eqIdx]; - } - - // restore the pointer to the FVElementGeometry - this->fvElemGeomPtr_ = oldFVGeometryPtr; - } - - /*! - * \brief Treat Dirichlet boundary conditions in a weak sense for a single - * intersection that only has Dirichlet boundary conditions - */ - template <class IntersectionIterator> - void evalDirichletSegment_(const IntersectionIterator &isIt, - const BoundaryTypes &bcTypes) - { - // temporary vector to store the Dirichlet boundary fluxes - PrimaryVariables values; - Valgrind::SetUndefined(values); - - unsigned bfIdx = isIt->indexInInside(); - - this->asImp_().computeFlux(values, bfIdx, true); - values *= this->curVolVars_(0).extrusionFactor(); - - // add fluxes to the residual - Valgrind::CheckDefined(values); - for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) - { - if (bcTypes.isDirichlet(eqIdx)) - this->residual_[0][eqIdx] += values[eqIdx]; - } - } - - /*! - * \brief Treat Dirichlet boundary conditions in a strong sense for a - * single intersection that has mixed D/N boundary conditions - */ - template <class IntersectionIterator> - void evalDirichletSegmentMixed_(const IntersectionIterator &isIt, - const BoundaryTypes &bcTypes) - { - // temporary vector to store the Dirichlet boundary fluxes - PrimaryVariables values; - Valgrind::SetUndefined(values); - - this->problem_().dirichlet(values, *isIt); - Valgrind::CheckDefined(values); - - // set Dirichlet conditions in a strong sense - for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) - { - if (bcTypes.isDirichlet(eqIdx)) - { - int pvIdx = bcTypes.eqToDirichletIndex(eqIdx); - this->residual_[0][eqIdx] - = this->curPriVar_(0, pvIdx) - values[pvIdx]; - } - } - } - - /*! - * \brief Add the flux terms to the local residual of the current element - */ - void evalFluxes_() - { - // calculate the mass flux over the faces and subtract - // it from the local rates - int fIdx = -1; - for (const auto& intersection : Dune::intersections(this->gridView_(), this->element_())) { - if (!intersection.neighbor()) - continue; - - fIdx++; - PrimaryVariables flux; - - Valgrind::SetUndefined(flux); - this->asImp_().computeFlux(flux, fIdx); - Valgrind::CheckDefined(flux); - - flux *= this->curVolVars_(0).extrusionFactor(); - - this->residual_[0] += flux; - } - } -}; - -} +#include <dumux/implicit/cellcentered/localresidual.hh> #endif diff --git a/dumux/implicit/cellcentered/ccproperties.hh b/dumux/implicit/cellcentered/ccproperties.hh index 28c8c32b02..c53d509333 100644 --- a/dumux/implicit/cellcentered/ccproperties.hh +++ b/dumux/implicit/cellcentered/ccproperties.hh @@ -1,52 +1,8 @@ -// -*- 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/>. * - *****************************************************************************/ -#ifndef DUMUX_CC_PROPERTIES_HH -#define DUMUX_CC_PROPERTIES_HH +#ifndef DUMUX_CC_PROPERTIES_HH_OLD +#define DUMUX_CC_PROPERTIES_HH_OLD -#include <dumux/implicit/properties.hh> +#warning this header is deprecated, use dumux/implicit/cellcentered/properties.hh instead -/*! - * \ingroup Properties - * \ingroup ImplicitProperties - * \ingroup CCModel - * \file - * \brief Specify the shape functions, operator assemblers, etc - * used for the CCModel. - */ -namespace Dumux -{ - -namespace Properties -{ -// \{ - -////////////////////////////////////////////////////////////////// -// Type tags -////////////////////////////////////////////////////////////////// - -//! The type tag for models based on the cell-centered scheme -NEW_TYPE_TAG(CCModel, INHERITS_FROM(ImplicitBase)); -} -} - -// \} - -#include "ccpropertydefaults.hh" +#include <dumux/implicit/cellcentered/properties.hh> #endif diff --git a/dumux/implicit/cellcentered/ccpropertydefaults.hh b/dumux/implicit/cellcentered/ccpropertydefaults.hh index 25ce3dcef2..3244780714 100644 --- a/dumux/implicit/cellcentered/ccpropertydefaults.hh +++ b/dumux/implicit/cellcentered/ccpropertydefaults.hh @@ -1,72 +1,8 @@ -// -*- 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/>. * - *****************************************************************************/ -/*! - * \ingroup Properties - * \ingroup CCProperties - * \ingroup CCModel - * \file - * - * \brief Default properties for cell centered models - */ -#ifndef DUMUX_CC_PROPERTY_DEFAULTS_HH -#define DUMUX_CC_PROPERTY_DEFAULTS_HH +#ifndef DUMUX_CC_PROPERTY_DEFAULTS_HH_OLD +#define DUMUX_CC_PROPERTY_DEFAULTS_HH_OLD -#include <dumux/implicit/propertydefaults.hh> -#include "ccassembler.hh" -#include "ccfvelementgeometry.hh" -#include "ccelementboundarytypes.hh" -#include "cclocalresidual.hh" -#include "ccelementvolumevariables.hh" -#include "ccproperties.hh" +#warning this header is deprecated, use dumux/implicit/cellcentered/propertydefaults.hh instead -namespace Dumux { - -// forward declarations -template<class TypeTag> class CCModel; -template<class TypeTag> class CCLocalResidual; -template<class TypeTag> class CCElementBoundaryTypes; -template<class TypeTag> class CCElementVolumeVariables; -template<class TypeTag> class CCFVElementGeometry; - -namespace Properties { -//! Set the default for the FVElementGeometry -SET_TYPE_PROP(CCModel, FVElementGeometry, Dumux::CCFVElementGeometry<TypeTag>); - -//! Set the default for the ElementBoundaryTypes -SET_TYPE_PROP(CCModel, ElementBoundaryTypes, Dumux::CCElementBoundaryTypes<TypeTag>); - -//! Mapper for the degrees of freedoms. -SET_TYPE_PROP(CCModel, DofMapper, typename GET_PROP_TYPE(TypeTag, ElementMapper)); - -//! Set the BaseLocalResidual to CCLocalResidual -SET_TYPE_PROP(CCModel, BaseLocalResidual, Dumux::CCLocalResidual<TypeTag>); - -//! An array of secondary variable containers -SET_TYPE_PROP(CCModel, ElementVolumeVariables, Dumux::CCElementVolumeVariables<TypeTag>); - -//! Assembler for the global jacobian matrix -SET_TYPE_PROP(CCModel, JacobianAssembler, Dumux::CCAssembler<TypeTag>); - -//! indicate that this is no box discretization -SET_BOOL_PROP(CCModel, ImplicitIsBox, false); - -} // namespace Properties -} // namespace Dumux +#include <dumux/implicit/cellcentered/propertydefaults.hh> #endif diff --git a/dumux/implicit/cellcentered/elementboundarytypes.hh b/dumux/implicit/cellcentered/elementboundarytypes.hh new file mode 100644 index 0000000000..c6248e9d1c --- /dev/null +++ b/dumux/implicit/cellcentered/elementboundarytypes.hh @@ -0,0 +1,140 @@ +// -*- 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 + * \brief Boundary types gathered on an element + */ +#ifndef DUMUX_CC_ELEMENT_BOUNDARY_TYPES_HH +#define DUMUX_CC_ELEMENT_BOUNDARY_TYPES_HH + +#include "properties.hh" + +#include <dumux/common/valgrind.hh> + +namespace Dumux +{ + +/*! + * \ingroup CCModel + * \ingroup ImplicitBoundaryTypes + * \brief This class stores an array of BoundaryTypes objects + */ +template<class TypeTag> +class CCElementBoundaryTypes : public std::vector<typename GET_PROP_TYPE(TypeTag, BoundaryTypes) > +{ + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef std::vector<BoundaryTypes> ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + + typedef typename GridView::template Codim<0>::Entity Element; + +public: + /*! + * \brief Copy constructor. + * + * Copying a the boundary types of an element should be explicitly + * requested + */ + explicit CCElementBoundaryTypes(const CCElementBoundaryTypes &v) + : ParentType(v) + {} + + /*! + * \brief Default constructor. + */ + CCElementBoundaryTypes() + { + hasDirichlet_ = false; + hasNeumann_ = false; + hasOutflow_ = false; + } + + /*! + * \brief Update the boundary types for all vertices of an element. + * + * \param problem The problem object which needs to be simulated + * \param element The DUNE Codim<0> entity for which the boundary + * types should be collected + */ + void update(const Problem &problem, + const Element &element) + { + this->resize(1); + + hasDirichlet_ = false; + hasNeumann_ = false; + hasOutflow_ = false; + + (*this)[0].reset(); + + if (!problem.model().onBoundary(element)) + return; + + for (const auto& intersection : Dune::intersections(problem.gridView(), element)) { + if (!intersection.boundary()) + continue; + + problem.boundaryTypes((*this)[0], intersection); + + hasDirichlet_ = hasDirichlet_ || (*this)[0].hasDirichlet(); + hasNeumann_ = hasNeumann_ || (*this)[0].hasNeumann(); + hasOutflow_ = hasOutflow_ || (*this)[0].hasOutflow(); + } + } + + void update(const Problem &problem, + const Element &element, + const FVElementGeometry &fvGeometry) + { + update(problem, element); + } + + /*! + * \brief Returns whether the element has a vertex which contains + * a Dirichlet value. + */ + bool hasDirichlet() const + { return hasDirichlet_; } + + /*! + * \brief Returns whether the element potentially features a + * Neumann boundary segment. + */ + bool hasNeumann() const + { return hasNeumann_; } + + /*! + * \brief Returns whether the element potentially features an + * outflow boundary segment. + */ + bool hasOutflow() const + { return hasOutflow_; } + +protected: + bool hasDirichlet_; + bool hasNeumann_; + bool hasOutflow_; +}; + +} // namespace Dumux + +#endif diff --git a/dumux/implicit/cellcentered/elementvolumevariables.hh b/dumux/implicit/cellcentered/elementvolumevariables.hh new file mode 100644 index 0000000000..54201aa192 --- /dev/null +++ b/dumux/implicit/cellcentered/elementvolumevariables.hh @@ -0,0 +1,142 @@ +// -*- 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 + * \brief Volume variables gathered on an element + */ +#ifndef DUMUX_CC_ELEMENT_VOLUME_VARIABLES_HH +#define DUMUX_CC_ELEMENT_VOLUME_VARIABLES_HH + +#include "properties.hh" + +namespace Dumux +{ + +/*! + * \ingroup CCModel + * \brief This class stores an array of VolumeVariables objects, one + * volume variables object for each of the element's vertices + */ +template<class TypeTag> +class CCElementVolumeVariables : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables) > +{ + typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GridView::template Codim<0>::Entity Element; + +public: + /*! + * \brief The constructor. + */ + CCElementVolumeVariables() + { } + + /*! + * \brief Construct the volume variables for all of vertices of an element. + * + * \param problem The problem which needs to be simulated. + * \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated + * \param fvGeometry The finite volume geometry of the element + * \param oldSol Tells whether the model's previous or current solution should be used. + */ + void update(const Problem &problem, + const Element &element, + const FVElementGeometry &fvGeometry, + bool oldSol) + { + const SolutionVector &globalSol = + oldSol? + problem.model().prevSol(): + problem.model().curSol(); + + int numNeighbors = fvGeometry.numNeighbors; + this->resize(numNeighbors); + + for (int i = 0; i < numNeighbors; i++) + { + const Element& neighbor = fvGeometry.neighbors[i]; + + const PrimaryVariables &solI + = globalSol[problem.elementMapper().index(neighbor)]; + + FVElementGeometry neighborFVGeom; + neighborFVGeom.updateInner(neighbor); + + (*this)[i].update(solI, + problem, + neighbor, + neighborFVGeom, + /*scvIdx=*/0, + oldSol); + } + + // only treat boundary if current solution is evaluated + if (!oldSol) + { + // check if element intersects with the boundary + ElementBoundaryTypes elemBCTypes; + elemBCTypes.update(problem, element); + if (elemBCTypes.hasDirichlet() + || elemBCTypes.hasNeumann() + || elemBCTypes.hasOutflow()) + { + this->resize(numNeighbors + element.subEntities(1)); + + // add volume variables for the boundary faces + for (const auto& intersection : Dune::intersections(problem.gridView(), element)) { + if (!intersection.boundary()) + continue; + + BoundaryTypes bcTypes; + problem.boundaryTypes(bcTypes, intersection); + + int fIdx = intersection.indexInInside(); + int indexInVariables = numNeighbors + fIdx; + + if (bcTypes.hasDirichlet()) + { + PrimaryVariables dirichletValues; + problem.dirichlet(dirichletValues, intersection); + + (*this)[indexInVariables].update(dirichletValues, + problem, + element, + fvGeometry, + /*scvIdx=*/0, + oldSol); + } + else + { + (*this)[indexInVariables] = (*this)[0]; + } + } + } + } + } +}; + +} // namespace Dumux + +#endif diff --git a/dumux/implicit/cellcentered/fvelementgeometry.hh b/dumux/implicit/cellcentered/fvelementgeometry.hh new file mode 100644 index 0000000000..0603825c95 --- /dev/null +++ b/dumux/implicit/cellcentered/fvelementgeometry.hh @@ -0,0 +1,235 @@ +// -*- 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 + * \brief Represents the finite volume geometry of a single element in + * the cell-centered fv scheme. + */ +#ifndef DUMUX_CC_FV_ELEMENTGEOMETRY_HH +#define DUMUX_CC_FV_ELEMENTGEOMETRY_HH + +#include <dune/geometry/affinegeometry.hh> +#include <dune/geometry/referenceelements.hh> +#include <dune/grid/common/intersectioniterator.hh> + +#include <dumux/common/propertysystem.hh> + +namespace Dumux +{ +namespace Properties +{ +NEW_PROP_TAG(GridView); +NEW_PROP_TAG(Scalar); +} + +/*! + * \ingroup CCModel + * \brief Represents the finite volume geometry of a single element in + * the cell-centered fv scheme. + */ +template<class TypeTag> +class CCFVElementGeometry +{ + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + enum{dim = GridView::dimension}; + enum{dimWorld = GridView::dimensionworld}; + + enum{maxNFAP = 2}; //! maximum number of flux approximation points (two-point flux) + enum{maxNE = 2*dim*(1 << (dim - 1))}; //! maximum number of neighbors (works for one hanging node per face) + enum{maxBF = 2*dim}; //! maximum number of boundary faces + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GridView::ctype CoordScalar; + typedef typename GridView::Traits::template Codim<0>::Entity Element; + typedef typename Element::Geometry Geometry; + typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition; + typedef Dune::FieldVector<CoordScalar,dim> LocalPosition; + typedef typename Dune::AffineGeometry<CoordScalar, dim, dimWorld> AffineGeometry; + +public: + struct SubControlVolume //! FV intersected with element + { + LocalPosition local; //!< local position + GlobalPosition global; //!< global position + Scalar volume; //!< volume of scv + bool inner; + }; + + struct SubControlVolumeFace //! interior face of a sub control volume + { + int i,j; //!< scvf seperates corner i and j of elem + LocalPosition ipLocal; //!< integration point in local coords + GlobalPosition ipGlobal; //!< integration point in global coords + GlobalPosition normal; //!< normal on face pointing to CV j or outward of the domain with length equal to |scvf| + Scalar area; //!< area of face + Dune::FieldVector<GlobalPosition, maxNFAP> grad; //!< derivatives of shape functions at ip + Dune::FieldVector<Scalar, maxNFAP> shapeValue; //!< value of shape functions at ip + Dune::FieldVector<int, maxNFAP> fapIndices; //!< indices w.r.t.neighbors of the flux approximation points + unsigned numFap; //!< number of flux approximation points + unsigned fIdx; //!< the index (w.r.t. the element) of the face (codim 1 entity) that the scvf is part of + }; + + typedef SubControlVolumeFace BoundaryFace; //!< compatibility typedef + + LocalPosition elementLocal; //!< local coordinate of element center + GlobalPosition elementGlobal; //!< global coordinate of element center + Scalar elementVolume; //!< element volume + SubControlVolume subContVol[1]; //!< data of the sub control volumes + std::vector<AffineGeometry> subContVolGeometries; //!< geometries of the subcontrol volumes + SubControlVolumeFace subContVolFace[maxNE]; //!< data of the sub control volume faces + BoundaryFace boundaryFace[maxBF]; //!< data of the boundary faces + int numScv; //!< number of subcontrol volumes + int numScvf; //!< number of inner-domain subcontrolvolume faces + int numNeighbors; //!< number of neighboring elements including the element itself + std::vector<Element> neighbors; //!< stores the neighboring elements + + void updateInner(const Element& element) + { + const Geometry geometry = element.geometry(); + + elementVolume = geometry.volume(); + elementGlobal = geometry.center(); + elementLocal = geometry.local(elementGlobal); + + numScv = 1; + numScvf = 0; + + subContVol[0].local = elementLocal; + subContVol[0].global = elementGlobal; + subContVol[0].inner = true; + subContVol[0].volume = elementVolume; + + // initialize neighbors list with self: + numNeighbors = 1; + neighbors.clear(); + neighbors.reserve(maxNE); + neighbors.push_back(element); + } + + void update(const GridView& gridView, const Element& element) + { + updateInner(element); + + const Geometry geometry = element.geometry(); + + bool onBoundary = false; + + // fill neighbor information and control volume face data: + for (const auto& intersection : Dune::intersections(gridView, element)) + { + const auto isGeometry = intersection.geometry(); + + // neighbor information and inner cvf data: + if (intersection.neighbor()) + { + numNeighbors++; + neighbors.push_back(intersection.outside()); + + int scvfIdx = numNeighbors - 2; + SubControlVolumeFace& scvFace = subContVolFace[scvfIdx]; + + scvFace.i = 0; + scvFace.j = scvfIdx + 1; + + scvFace.ipGlobal = isGeometry.center(); + scvFace.ipLocal = geometry.local(scvFace.ipGlobal); + scvFace.normal = intersection.centerUnitOuterNormal(); + Scalar volume = isGeometry.volume(); + scvFace.normal *= volume; + scvFace.area = volume; + + GlobalPosition distVec = elementGlobal + - neighbors[scvfIdx+1].geometry().center(); + distVec /= distVec.two_norm2(); + + // gradients using a two-point flux approximation + scvFace.numFap = 2; + for (unsigned int fapIdx = 0; fapIdx < scvFace.numFap; fapIdx++) + { + scvFace.grad[fapIdx] = distVec; + scvFace.shapeValue[fapIdx] = 0.5; + } + scvFace.grad[1] *= -1.0; + + scvFace.fapIndices[0] = scvFace.i; + scvFace.fapIndices[1] = scvFace.j; + + scvFace.fIdx = intersection.indexInInside(); + } + + // boundary cvf data + if (intersection.boundary()) + { + onBoundary = true; + int bfIdx = intersection.indexInInside(); + SubControlVolumeFace& bFace = boundaryFace[bfIdx]; + + bFace.ipGlobal = isGeometry.center(); + bFace.ipLocal = geometry.local(bFace.ipGlobal); + bFace.normal = intersection.centerUnitOuterNormal(); + Scalar volume = isGeometry.volume(); + bFace.normal *= volume; + bFace.area = volume; + bFace.i = 0; + bFace.j = 0; + + GlobalPosition distVec = elementGlobal - bFace.ipGlobal; + distVec /= distVec.two_norm2(); + + // gradients using a two-point flux approximation + bFace.numFap = 2; + for (unsigned int fapIdx = 0; fapIdx < bFace.numFap; fapIdx++) + { + bFace.grad[fapIdx] = distVec; + bFace.shapeValue[fapIdx] = 0.5; + } + bFace.grad[1] *= -1.0; + + bFace.fapIndices[0] = bFace.i; + bFace.fapIndices[1] = bFace.j; + } + } + + // set the number of inner-domain subcontrolvolume faces + numScvf = numNeighbors - 1; + + // treat elements on the boundary + if (onBoundary) + { + for (int bfIdx = 0; bfIdx < element.subEntities(1); bfIdx++) + { + SubControlVolumeFace& bFace = boundaryFace[bfIdx]; + bFace.j = numNeighbors + bfIdx; + bFace.fapIndices[1] = bFace.j; + neighbors.push_back(element); + } + } + } + + /*! + * \brief For compatibilty with the box element geometry + */ + int boundaryFaceIndex(const int fIdx, const int vIdxInFace) const + { + return fIdx; + } +}; + +} + +#endif diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh new file mode 100644 index 0000000000..3f361adc0e --- /dev/null +++ b/dumux/implicit/cellcentered/localresidual.hh @@ -0,0 +1,288 @@ +// -*- 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 + * \brief Calculates the residual of models based on the box scheme element-wise. + */ +#ifndef DUMUX_CC_LOCAL_RESIDUAL_HH +#define DUMUX_CC_LOCAL_RESIDUAL_HH + +#include <dune/istl/matrix.hh> + +#include <dumux/common/valgrind.hh> +#include <dumux/implicit/localresidual.hh> + +#include "properties.hh" + +namespace Dumux +{ +/*! + * \ingroup CCModel + * \ingroup CCLocalResidual + * \brief Element-wise calculation of the residual for models + * based on the fully implicit cell-centered scheme. + * + * \todo Please doc me more! + */ +template<class TypeTag> +class CCLocalResidual : public ImplicitLocalResidual<TypeTag> +{ + typedef ImplicitLocalResidual<TypeTag> ParentType; + friend class ImplicitLocalResidual<TypeTag>; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + enum { + numEq = GET_PROP_VALUE(TypeTag, NumEq) + }; + + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename FVElementGeometry::SubControlVolumeFace SCVFace; + + // copying the local residual class is not a good idea + CCLocalResidual(const CCLocalResidual &); + +public: + CCLocalResidual() : ParentType() + { } + +protected: + + /*! + * \brief Add all fluxes resulting from Neumann, outflow and pure Dirichlet + * boundary conditions to the local residual. + */ + void evalBoundaryFluxes_() + { + for (const auto& intersection : Dune::intersections(this->gridView_(), this->element_())) + { + // handle only faces on the boundary + if (!intersection.boundary()) + continue; + + BoundaryTypes bcTypes; + this->problem_().boundaryTypes(bcTypes, intersection); + + // evaluate the Neumann conditions at the boundary face + if (bcTypes.hasNeumann()) + this->asImp_().evalNeumannSegment_(&intersection, bcTypes); + + // evaluate the outflow conditions at the boundary face + if (bcTypes.hasOutflow()) + this->asImp_().evalOutflowSegment_(&intersection, bcTypes); + + // evaluate the pure Dirichlet conditions at the boundary face + if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann()) + this->asImp_().evalDirichletSegment_(&intersection, bcTypes); + } + } + + /*! + * \brief Evaluate Dirichlet conditions on faces that have mixed + * Dirichlet/Neumann boundary conditions. + */ + void evalDirichlet_() + { + for (const auto& intersection : Dune::intersections(this->gridView_(), this->element_())) + { + // handle only faces on the boundary + if (!intersection.boundary()) + continue; + + BoundaryTypes bcTypes; + this->problem_().boundaryTypes(bcTypes, intersection); + + if (bcTypes.hasDirichlet() && bcTypes.hasNeumann()) + this->asImp_().evalDirichletSegmentMixed_(&intersection, bcTypes); + } + } + + /*! + * \brief Add Neumann boundary conditions for a single intersection + */ + template <class IntersectionIterator> + void evalNeumannSegment_(const IntersectionIterator &isIt, + const BoundaryTypes &bcTypes) + { + // temporary vector to store the neumann boundary fluxes + PrimaryVariables values; + Valgrind::SetUndefined(values); + + unsigned bfIdx = isIt->indexInInside(); + this->problem_().solDependentNeumann(values, + this->element_(), + this->fvGeometry_(), + *isIt, + /*scvIdx=*/0, + bfIdx, + this->curVolVars_()); + values *= isIt->geometry().volume() + * this->curVolVars_(0).extrusionFactor(); + + // add fluxes to the residual + Valgrind::CheckDefined(values); + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { + if (bcTypes.isNeumann(eqIdx)) + this->residual_[0][eqIdx] += values[eqIdx]; + } + } + + /*! + * \brief Add outflow boundary conditions for a single intersection + */ + template <class IntersectionIterator> + void evalOutflowSegment_(const IntersectionIterator &isIt, + const BoundaryTypes &bcTypes) + { + if (this->element_().geometry().type().isCube() == false) + DUNE_THROW(Dune::InvalidStateException, + "for cell-centered models, outflow BCs only work for cubes."); + + // store pointer to the current FVElementGeometry + const FVElementGeometry *oldFVGeometryPtr = this->fvElemGeomPtr_; + + // copy the current FVElementGeometry to a local variable + // and set the pointer to this local variable + FVElementGeometry fvGeometry = this->fvGeometry_(); + this->fvElemGeomPtr_ = &fvGeometry; + + // get the index of the boundary face + unsigned bfIdx = isIt->indexInInside(); + unsigned oppositeIdx = bfIdx^1; + + // manipulate the corresponding subcontrolvolume face + SCVFace& boundaryFace = fvGeometry.boundaryFace[bfIdx]; + + // set the second flux approximation index for the boundary face + for (int nIdx = 0; nIdx < fvGeometry.numNeighbors-1; nIdx++) + { + // check whether the two faces are opposite of each other + if (fvGeometry.subContVolFace[nIdx].fIdx == oppositeIdx) + { + boundaryFace.j = nIdx+1; + break; + } + } + boundaryFace.fapIndices[1] = boundaryFace.j; + boundaryFace.grad[0] *= -0.5; + boundaryFace.grad[1] *= -0.5; + + // temporary vector to store the outflow boundary fluxes + PrimaryVariables values; + Valgrind::SetUndefined(values); + + this->asImp_().computeFlux(values, bfIdx, true); + values *= this->curVolVars_(0).extrusionFactor(); + + // add fluxes to the residual + Valgrind::CheckDefined(values); + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) + { + if (bcTypes.isOutflow(eqIdx)) + this->residual_[0][eqIdx] += values[eqIdx]; + } + + // restore the pointer to the FVElementGeometry + this->fvElemGeomPtr_ = oldFVGeometryPtr; + } + + /*! + * \brief Treat Dirichlet boundary conditions in a weak sense for a single + * intersection that only has Dirichlet boundary conditions + */ + template <class IntersectionIterator> + void evalDirichletSegment_(const IntersectionIterator &isIt, + const BoundaryTypes &bcTypes) + { + // temporary vector to store the Dirichlet boundary fluxes + PrimaryVariables values; + Valgrind::SetUndefined(values); + + unsigned bfIdx = isIt->indexInInside(); + + this->asImp_().computeFlux(values, bfIdx, true); + values *= this->curVolVars_(0).extrusionFactor(); + + // add fluxes to the residual + Valgrind::CheckDefined(values); + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) + { + if (bcTypes.isDirichlet(eqIdx)) + this->residual_[0][eqIdx] += values[eqIdx]; + } + } + + /*! + * \brief Treat Dirichlet boundary conditions in a strong sense for a + * single intersection that has mixed D/N boundary conditions + */ + template <class IntersectionIterator> + void evalDirichletSegmentMixed_(const IntersectionIterator &isIt, + const BoundaryTypes &bcTypes) + { + // temporary vector to store the Dirichlet boundary fluxes + PrimaryVariables values; + Valgrind::SetUndefined(values); + + this->problem_().dirichlet(values, *isIt); + Valgrind::CheckDefined(values); + + // set Dirichlet conditions in a strong sense + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) + { + if (bcTypes.isDirichlet(eqIdx)) + { + int pvIdx = bcTypes.eqToDirichletIndex(eqIdx); + this->residual_[0][eqIdx] + = this->curPriVar_(0, pvIdx) - values[pvIdx]; + } + } + } + + /*! + * \brief Add the flux terms to the local residual of the current element + */ + void evalFluxes_() + { + // calculate the mass flux over the faces and subtract + // it from the local rates + int fIdx = -1; + for (const auto& intersection : Dune::intersections(this->gridView_(), this->element_())) { + if (!intersection.neighbor()) + continue; + + fIdx++; + PrimaryVariables flux; + + Valgrind::SetUndefined(flux); + this->asImp_().computeFlux(flux, fIdx); + Valgrind::CheckDefined(flux); + + flux *= this->curVolVars_(0).extrusionFactor(); + + this->residual_[0] += flux; + } + } +}; + +} + +#endif diff --git a/dumux/implicit/cellcentered/properties.hh b/dumux/implicit/cellcentered/properties.hh new file mode 100644 index 0000000000..75d3d85ace --- /dev/null +++ b/dumux/implicit/cellcentered/properties.hh @@ -0,0 +1,52 @@ +// -*- 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/>. * + *****************************************************************************/ +#ifndef DUMUX_CC_PROPERTIES_HH +#define DUMUX_CC_PROPERTIES_HH + +#include <dumux/implicit/properties.hh> + +/*! + * \ingroup Properties + * \ingroup ImplicitProperties + * \ingroup CCModel + * \file + * \brief Specify the shape functions, operator assemblers, etc + * used for the CCModel. + */ +namespace Dumux +{ + +namespace Properties +{ +// \{ + +////////////////////////////////////////////////////////////////// +// Type tags +////////////////////////////////////////////////////////////////// + +//! The type tag for models based on the cell-centered scheme +NEW_TYPE_TAG(CCModel, INHERITS_FROM(ImplicitBase)); +} +} + +// \} + +#include "propertydefaults.hh" + +#endif diff --git a/dumux/implicit/cellcentered/propertydefaults.hh b/dumux/implicit/cellcentered/propertydefaults.hh new file mode 100644 index 0000000000..72f88bcca2 --- /dev/null +++ b/dumux/implicit/cellcentered/propertydefaults.hh @@ -0,0 +1,72 @@ +// -*- 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/>. * + *****************************************************************************/ +/*! + * \ingroup Properties + * \ingroup CCProperties + * \ingroup CCModel + * \file + * + * \brief Default properties for cell centered models + */ +#ifndef DUMUX_CC_PROPERTY_DEFAULTS_HH +#define DUMUX_CC_PROPERTY_DEFAULTS_HH + +#include <dumux/implicit/propertydefaults.hh> +#include "assembler.hh" +#include "fvelementgeometry.hh" +#include "elementboundarytypes.hh" +#include "localresidual.hh" +#include "elementvolumevariables.hh" +#include "properties.hh" + +namespace Dumux { + +// forward declarations +template<class TypeTag> class CCModel; +template<class TypeTag> class CCLocalResidual; +template<class TypeTag> class CCElementBoundaryTypes; +template<class TypeTag> class CCElementVolumeVariables; +template<class TypeTag> class CCFVElementGeometry; + +namespace Properties { +//! Set the default for the FVElementGeometry +SET_TYPE_PROP(CCModel, FVElementGeometry, Dumux::CCFVElementGeometry<TypeTag>); + +//! Set the default for the ElementBoundaryTypes +SET_TYPE_PROP(CCModel, ElementBoundaryTypes, Dumux::CCElementBoundaryTypes<TypeTag>); + +//! Mapper for the degrees of freedoms. +SET_TYPE_PROP(CCModel, DofMapper, typename GET_PROP_TYPE(TypeTag, ElementMapper)); + +//! Set the BaseLocalResidual to CCLocalResidual +SET_TYPE_PROP(CCModel, BaseLocalResidual, Dumux::CCLocalResidual<TypeTag>); + +//! An array of secondary variable containers +SET_TYPE_PROP(CCModel, ElementVolumeVariables, Dumux::CCElementVolumeVariables<TypeTag>); + +//! Assembler for the global jacobian matrix +SET_TYPE_PROP(CCModel, JacobianAssembler, Dumux::CCAssembler<TypeTag>); + +//! indicate that this is no box discretization +SET_BOOL_PROP(CCModel, ImplicitIsBox, false); + +} // namespace Properties +} // namespace Dumux + +#endif diff --git a/dumux/implicit/cornerpoint/cpelementvolumevariables.hh b/dumux/implicit/cornerpoint/cpelementvolumevariables.hh index 00913e6ad6..7a8da16acb 100644 --- a/dumux/implicit/cornerpoint/cpelementvolumevariables.hh +++ b/dumux/implicit/cornerpoint/cpelementvolumevariables.hh @@ -23,7 +23,7 @@ #ifndef DUMUX_CP_ELEMENT_VOLUME_VARIABLES_HH #define DUMUX_CP_ELEMENT_VOLUME_VARIABLES_HH -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> namespace Dumux { diff --git a/dumux/implicit/mpnc/mpncproperties.hh b/dumux/implicit/mpnc/mpncproperties.hh index c6d99cf100..3943d328ca 100644 --- a/dumux/implicit/mpnc/mpncproperties.hh +++ b/dumux/implicit/mpnc/mpncproperties.hh @@ -20,7 +20,7 @@ #define DUMUX_MPNC_PROPERTIES_HH #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> /*! * \ingroup Properties diff --git a/dumux/implicit/richards/richardsproperties.hh b/dumux/implicit/richards/richardsproperties.hh index 5425891e3f..fec3d44daf 100644 --- a/dumux/implicit/richards/richardsproperties.hh +++ b/dumux/implicit/richards/richardsproperties.hh @@ -29,7 +29,7 @@ #define DUMUX_RICHARDS_PROPERTIES_HH #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> #include <dumux/implicit/nonisothermal/niproperties.hh> namespace Dumux diff --git a/dumux/linear/amgparallelhelpers.hh b/dumux/linear/amgparallelhelpers.hh index dbef9da7b4..3cc263c19e 100644 --- a/dumux/linear/amgparallelhelpers.hh +++ b/dumux/linear/amgparallelhelpers.hh @@ -26,7 +26,7 @@ #define DUMUX_AMGPARALLELHELPERS_HH #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> #include <dumux/decoupled/common/pressureproperties.hh> #include <dumux/linear/amgproperties.hh> diff --git a/dumux/linear/amgproperties.hh b/dumux/linear/amgproperties.hh index cab57bbb9c..86ce952fa1 100644 --- a/dumux/linear/amgproperties.hh +++ b/dumux/linear/amgproperties.hh @@ -36,7 +36,7 @@ #include <dune/common/version.hh> #include <dumux/implicit/box/properties.hh> -#include <dumux/implicit/cellcentered/ccproperties.hh> +#include <dumux/implicit/cellcentered/properties.hh> #include <dumux/decoupled/common/decoupledproperties.hh> #include <dumux/decoupled/common/pressureproperties.hh> #include "linearsolverproperties.hh" diff --git a/test/porousmediumflow/2p/implicit/cc2pcornerpointproblem.hh b/test/porousmediumflow/2p/implicit/cc2pcornerpointproblem.hh index e395d5952a..25cfbc24f2 100644 --- a/test/porousmediumflow/2p/implicit/cc2pcornerpointproblem.hh +++ b/test/porousmediumflow/2p/implicit/cc2pcornerpointproblem.hh @@ -23,7 +23,7 @@ #include <dumux/material/components/dnapl.hh> #include <dumux/implicit/2p/2pmodel.hh> #include <dumux/porousmediumflow/implicit/problem.hh> -#include <dumux/implicit/cellcentered/ccpropertydefaults.hh> +#include <dumux/implicit/cellcentered/propertydefaults.hh> #include <dumux/linear/amgbackend.hh> #include <dune/grid/CpGrid.hpp> diff --git a/test/porousmediumflow/2p/implicit/lensproblem.hh b/test/porousmediumflow/2p/implicit/lensproblem.hh index ebc0d40a5a..e6c2056adc 100644 --- a/test/porousmediumflow/2p/implicit/lensproblem.hh +++ b/test/porousmediumflow/2p/implicit/lensproblem.hh @@ -40,7 +40,7 @@ #include <dumux/material/components/dnapl.hh> #include <dumux/implicit/2p/2pmodel.hh> #include <dumux/porousmediumflow/implicit/problem.hh> -#include <dumux/implicit/cellcentered/ccpropertydefaults.hh> +#include <dumux/implicit/cellcentered/propertydefaults.hh> #include <dumux/implicit/2p/2pgridadaptindicator.hh> #include <dumux/implicit/adaptive/gridadaptinitializationindicator.hh> -- GitLab