staggeredfvproblem.hh

Go to the documentation of this file.

// -*- 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 3 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_STAGGERD_FV_PROBLEM_HH
#define DUMUX_STAGGERD_FV_PROBLEM_HH
 
#include <dune/common/rangeutilities.hh>
 
#include <dumux/common/properties.hh>
#include <dumux/common/fvproblem.hh>
#include <dumux/common/numeqvector.hh>
 
namespace Dumux {
 
template<class TypeTag>
class StaggeredFVProblem : public FVProblem<TypeTag>
{
    using ParentType = FVProblem<TypeTag>;
    using Implementation = GetPropType<TypeTag, Properties::Problem>;
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
 
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
    using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
    using GridFaceVariables = typename GridVariables::GridFaceVariables;
    using ElementFaceVariables = typename GridFaceVariables::LocalView;
 
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GridGeometry::LocalView;
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
 
    using CoordScalar = typename GridView::ctype;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
    static constexpr auto cellCenterIdx = GridGeometry::cellCenterIdx();
    static constexpr auto faceIdx = GridGeometry::faceIdx();
 
    static constexpr auto numEqCellCenter = getPropValue<TypeTag, Properties::NumEqCellCenter>();
    static constexpr auto numEqFace = getPropValue<TypeTag, Properties::NumEqFace>();
 
public:
    StaggeredFVProblem(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
    : ParentType(gridGeometry, paramGroup)
    { }
 
    bool isDirichletCell(const Element& element,
                         const FVElementGeometry& fvGeometry,
                         const SubControlVolume& scv,
                         int pvIdx) const
    { return false; }
 
    template<class ElementVolumeVariables, class ElementFaceVariables, class Entity>
    NumEqVector source(const Element &element,
                       const FVElementGeometry& fvGeometry,
                       const ElementVolumeVariables& elemVolVars,
                       const ElementFaceVariables& elementFaceVars,
                       const Entity &e) const
    {
        // forward to solution independent, fully-implicit specific interface
        return asImp_().sourceAtPos(e.center());
    }
 
    NumEqVector neumann(const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
                        const ElementFaceVariables& elemFaceVars,
                        const SubControlVolumeFace& scvf) const
    {
        // forward it to the interface with only the global position
        return asImp_().neumannAtPos(scvf.ipGlobal());
    }
 
    template<class Entity>
    PrimaryVariables initial(const Entity& entity) const
    {
        return asImp_().initialAtPos(entity.center());
    }
 
    template<class SolutionVector>
    void applyInitialSolution(SolutionVector& sol) const
    {
        sol[cellCenterIdx].resize(this->gridGeometry().numCellCenterDofs());
        sol[faceIdx].resize(this->gridGeometry().numFaceDofs());
 
        for (const auto& element : elements(this->gridGeometry().gridView()))
        {
            auto fvGeometry = localView(this->gridGeometry());
            fvGeometry.bindElement(element);
 
            // loop over sub control volumes
            for (auto&& scv : scvs(fvGeometry))
            {
                // let the problem do the dirty work of nailing down
                // the initial solution.
                auto initPriVars = asImp_().initial(scv);
                asImp_().applyInitialCellCenterSolution(sol, scv, initPriVars);
            }
 
            // loop over faces
            for(auto&& scvf : scvfs(fvGeometry))
            {
                auto initPriVars = asImp_().initial(scvf);
                asImp_().applyInitialFaceSolution(sol, scvf, initPriVars);
            }
        }
    }
 
 
    template<class SolutionVector>
    void applyInitialCellCenterSolution(SolutionVector& sol,
                                        const SubControlVolume& scv,
                                        const PrimaryVariables& initSol) const
    {
        // while the container within the actual solution vector holds numEqCellCenter
        // elements, we need to specify an offset to get the correct entry of the initial solution
        static constexpr auto offset = PrimaryVariables::dimension - numEqCellCenter;
 
        for(int pvIdx = 0; pvIdx < numEqCellCenter; ++pvIdx)
            sol[cellCenterIdx][scv.dofIndex()][pvIdx] = initSol[pvIdx + offset];
    }
 
    template<class SolutionVector>
    void applyInitialFaceSolution(SolutionVector& sol,
                                  const SubControlVolumeFace& scvf,
                                  const PrimaryVariables& initSol) const
    {
        for(int pvIdx = 0; pvIdx < numEqFace; ++pvIdx)
            sol[faceIdx][scvf.dofIndex()][pvIdx] = initSol[pvIdx];
    }
 
protected:
    Implementation &asImp_()
    { return *static_cast<Implementation *>(this); }
 
    const Implementation &asImp_() const
    { return *static_cast<const Implementation *>(this); }
};
 
} // end namespace Dumux
 
#endif