fluxstencil.hh 7.68 KB
Newer Older
Timo Koch's avatar
Timo Koch committed
1
2
3
4
5
6
7
// -*- 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    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
Timo Koch's avatar
Timo Koch committed
9
10
11
12
13
14
15
16
17
18
19
20
 *   (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
21
22
 * \ingroup Discretization
 * \brief The flux stencil specialized for different discretization schemes
Timo Koch's avatar
Timo Koch committed
23
24
25
26
 */
#ifndef DUMUX_DISCRETIZATION_FLUXSTENCIL_HH
#define DUMUX_DISCRETIZATION_FLUXSTENCIL_HH

27
28
#include <vector>

29
#include <dune/common/reservedvector.hh>
30
#include <dumux/common/indextraits.hh>
31
#include <dumux/discretization/method.hh>
Timo Koch's avatar
Timo Koch committed
32

33
namespace Dumux {
Timo Koch's avatar
Timo Koch committed
34
35
36

/*!
 * \ingroup Discretization
37
 * \brief The flux stencil specialized for different discretization schemes
Timo Koch's avatar
Timo Koch committed
38
39
40
41
42
43
 * \note There might be different stencils used for e.g. advection and diffusion for schemes
 *       where the stencil depends on variables. Also schemes might even have solution dependent
 *       stencil. However, we always reserve the stencil or all DOFs that are possibly involved
 *       since we use the flux stencil for matrix and assembly. This might lead to some zeros stored
 *       in the matrix.
 */
44
template<class FVElementGeometry, DiscretizationMethod discMethod = FVElementGeometry::GridGeometry::discMethod>
45
class FluxStencil;
Timo Koch's avatar
Timo Koch committed
46

47
48
49
/*
 * \ingroup Discretization
 * \brief Flux stencil specialization for the cell-centered tpfa scheme
50
 * \tparam FVElementGeometry The local view on the finite volume grid geometry
51
 */
52
template<class FVElementGeometry>
53
class FluxStencil<FVElementGeometry, DiscretizationMethod::cctpfa>
Timo Koch's avatar
Timo Koch committed
54
{
55
56
57
    using GridGeometry = typename FVElementGeometry::GridGeometry;
    using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
    using GridView = typename GridGeometry::GridView;
Timo Koch's avatar
Timo Koch committed
58
    using Element = typename GridView::template Codim<0>::Entity;
59
    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
Timo Koch's avatar
Timo Koch committed
60
61

public:
62
    //! Each cell I couples to a cell J always only via one face
63
    using ScvfStencilIForJ = Dune::ReservedVector<GridIndexType, 1>;
64
65

    //! The flux stencil type
66
    using Stencil = typename SubControlVolumeFace::Traits::GridIndexStorage;
67

68
    //! Returns the flux stencil
69
    static Stencil stencil(const Element& element,
Timo Koch's avatar
Timo Koch committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
                           const FVElementGeometry& fvGeometry,
                           const SubControlVolumeFace& scvf)
    {
        if (scvf.boundary())
            return Stencil({scvf.insideScvIdx()});
        else if (scvf.numOutsideScvs() > 1)
        {
            Stencil stencil({scvf.insideScvIdx()});
            for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
                stencil.push_back(scvf.outsideScvIdx(i));
            return stencil;
        }
        else
            return Stencil({scvf.insideScvIdx(), scvf.outsideScvIdx()});
    }
};

87
88
89
/*
 * \ingroup Discretization
 * \brief Flux stencil specialization for the cell-centered mpfa scheme
90
 * \tparam FVElementGeometry The local view on the finite volume grid geometry
91
 */
92
template<class FVElementGeometry>
93
class FluxStencil<FVElementGeometry, DiscretizationMethod::ccmpfa>
Timo Koch's avatar
Timo Koch committed
94
{
95
96
97
    using GridGeometry = typename FVElementGeometry::GridGeometry;
    using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
    using GridView = typename GridGeometry::GridView;
Timo Koch's avatar
Timo Koch committed
98
    using Element = typename GridView::template Codim<0>::Entity;
99
    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
100
101

    // Use the stencil type of the primary interaction volume
102
    using NodalIndexSet = typename GridGeometry::GridIVIndexSets::DualGridIndexSet::NodalIndexSet;
Timo Koch's avatar
Timo Koch committed
103
104

public:
105
    //! We don't know yet how many faces couple to a neighboring element
106
    using ScvfStencilIForJ = std::vector<GridIndexType>;
107
108

    //! The flux stencil type
109
    using Stencil = typename NodalIndexSet::NodalGridStencilType;
110

111
    //! Returns the indices of the elements required for flux calculation on an scvf.
112
113
114
    static const Stencil& stencil(const Element& element,
                                  const FVElementGeometry& fvGeometry,
                                  const SubControlVolumeFace& scvf)
Timo Koch's avatar
Timo Koch committed
115
    {
116
        const auto& gridGeometry = fvGeometry.gridGeometry();
Timo Koch's avatar
Timo Koch committed
117
118

        // return the scv (element) indices in the interaction region
119
120
        if (gridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
            return gridGeometry.gridInteractionVolumeIndexSets().secondaryIndexSet(scvf).gridScvIndices();
Timo Koch's avatar
Timo Koch committed
121
        else
122
            return gridGeometry.gridInteractionVolumeIndexSets().primaryIndexSet(scvf).gridScvIndices();
Timo Koch's avatar
Timo Koch committed
123
124
125
    }
};

126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
/*
 * \ingroup Discretization
 * \brief Flux stencil specialization for the cell-centered wmpfa scheme
 * \tparam FVElementGeometry The local view on the finite volume grid geometry
 */
template<class FVElementGeometry>
class FluxStencil<FVElementGeometry, DiscretizationMethod::ccwmpfa>
{
    using GridGeometry = typename FVElementGeometry::GridGeometry;
    using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
    using GridView = typename GridGeometry::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
    using GridIndexType = typename IndexTraits<GridView>::GridIndex;

public:
    //! We don't know yet how many faces couple to a neighboring element
    //! Each cell I couples to a cell J always only via one face
    using ScvfStencilIForJ = std::vector<GridIndexType>;

    //! The flux stencil type
    using Stencil = typename SubControlVolumeFace::Traits::GridIndexStorage;

    //! Returns the indices of the elements required for flux calculation on an scvf.
    static Stencil stencil(const Element& element,
                                  const FVElementGeometry& fvGeometry,
                                  const SubControlVolumeFace& scvf)
    {
        const auto& gridGeometry = fvGeometry.gridGeometry();

        //ToDo correct stencil! At the moment simply all neighors of neighbors are added
        Stencil stencil({scvf.insideScvIdx()});
        for (const auto& scvf : scvfs(fvGeometry))
        {
            if (scvf.boundary())
                continue;

            stencil.push_back(scvf.outsideScvIdx());

            const auto outsideScvIdx = scvf.outsideScvIdx();
            const auto outsideElement = gridGeometry.element(outsideScvIdx);
            auto fvGeometryJ = localView(gridGeometry);
            fvGeometryJ.bind(outsideElement);

            for (const auto& scvfJ : scvfs(fvGeometryJ))
            {
                if (scvfJ.boundary() || scvf.insideScvIdx() == scvfJ.outsideScvIdx())
                    continue;

                stencil.push_back(scvfJ.outsideScvIdx());
            }
        }

        return stencil;
    }
};

Timo Koch's avatar
Timo Koch committed
182
183
184
} // end namespace Dumux

#endif