boundarytypes.hh 11.9 KB
Newer Older
1
2
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
3
/*****************************************************************************
4
 *   See the file COPYING for full copying permissions.                      *
5
 *                                                                           *
6
 *   This program is free software: you can redistribute it and/or modify    *
7
 *   it under the terms of the GNU General Public License as published by    *
8
9
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
10
 *                                                                           *
11
12
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
14
15
16
17
 *   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/>.   *
18
19
20
 *****************************************************************************/
/*!
 * \file
21
 * \ingroup Common
22
23
 * \brief Class to specify the type of a boundary.
 */
24
25
#ifndef DUMUX_BOUNDARY_TYPES_HH
#define DUMUX_BOUNDARY_TYPES_HH
26

Timo Koch's avatar
Timo Koch committed
27
28
29
#include <algorithm>
#include <array>

30
31
32
namespace Dumux
{

33
/*!
34
 * \ingroup Common
35
36
 * \brief Class to specify the type of a boundary.
 */
37
38
39
40
41
42
43
template <int numEq>
class BoundaryTypes
{
public:
    BoundaryTypes()
    { reset(); }

44
45
46
47
    //! we have a boundary condition for each equation
    static constexpr int size()
    { return numEq; }

48
    /*!
49
     * \brief Reset the boundary types for all equations.
50
51
     *
     * After this method no equations will be disabled and neither
52
     * Neumann nor Dirichlet conditions will be evaluated. This
53
     * corresponds to a Neumann zero boundary.
54
55
56
     */
    void reset()
    {
57
58
        for (int eqIdx=0; eqIdx < numEq; ++eqIdx)
            resetEq(eqIdx);
59
60
    }

61
62
63
64
65
    /*!
     * \brief Reset the boundary types for one equation.
     */
    void resetEq(int eqIdx)
    {
66
          boundaryInfo_[eqIdx].visited = false;
67

68
69
70
71
72
          boundaryInfo_[eqIdx].isDirichlet = false;
          boundaryInfo_[eqIdx].isNeumann = false;
          boundaryInfo_[eqIdx].isOutflow = false;
          boundaryInfo_[eqIdx].isCouplingDirichlet = false;
          boundaryInfo_[eqIdx].isCouplingNeumann = false;
73
74
75
76
77

          eq2pvIdx_[eqIdx] = eqIdx;
          pv2eqIdx_[eqIdx] = eqIdx;
    }

78
    /*!
79
80
81
82
     * \brief Returns true if the boundary types for a given equation
     *        has been specified.
     *
     * \param eqIdx The index of the equation
83
84
     */
    bool isSet(int eqIdx) const
85
    { return boundaryInfo_[eqIdx].visited; }
86
87
88
89

    /*!
     * \brief Make sure the boundary conditions are well-posed.
     *
90
     * If they are not, an assertion fails and the program aborts!
91
     * (if the NDEBUG macro is not defined)
92
93
94
     */
    void checkWellPosed() const
    {
95
        // if this fails, at least one condition is missing.
96
        for (int i=0; i < numEq; ++i)
97
            assert(boundaryInfo_[i].visited);
98
    }
99
100

    /*!
101
     * \brief Set all boundary conditions to Neumann.
102
103
104
     */
    void setAllNeumann()
    {
105
106
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
            setNeumann(eqIdx);
107
108
109
    }

    /*!
110
     * \brief Set all boundary conditions to Dirichlet.
111
112
113
     */
    void setAllDirichlet()
    {
114
115
        for (int eqIdx = 0; eqIdx < numEq; ++ eqIdx)
            setDirichlet(eqIdx);
116
117
    }

118
    /*!
119
     * \brief Set all boundary conditions to Neumann.
120
121
122
     */
    void setAllOutflow()
    {
123
124
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
            setOutflow(eqIdx);
125
126
    }

127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
    /*!
     * \brief Set all boundary conditions to Dirichlet-like coupling
     */
    void setAllCouplingDirichlet()
    {
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
            setCouplingDirichlet(eqIdx);
    }

    /*!
     * \brief Set all boundary conditions to Neumann-like coupling.
     */
    void setAllCouplingNeumann()
    {
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
            setCouplingNeumann(eqIdx);
    }

145
    /*!
146
     * \brief Set a Neumann boundary condition for a single a single
147
     *        equation.
148
149
     *
     * \param eqIdx The index of the equation
150
151
152
     */
    void setNeumann(int eqIdx)
    {
153
        resetEq(eqIdx);
154
155
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isNeumann = true;
156
157
158
    }

    /*!
159
     * \brief Set a Dirichlet boundary condition for a single primary
160
161
     *        variable
     *
162
163
164
165
     * \param pvIdx The index of the primary variable for which the
     *              Dirichlet condition should apply.
     * \param eqIdx The index of the equation which should used to set
     *              the Dirichlet condition
166
     */
167
    void setDirichlet(int pvIdx, int eqIdx)
168
    {
169
        resetEq(eqIdx);
170
171
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isDirichlet = true;
172
173
174
175
176
177

        // update the equation <-> primary variable mapping
        eq2pvIdx_[eqIdx] = pvIdx;
        pv2eqIdx_[pvIdx] = eqIdx;
    }

178
    /*!
179
     * \brief Set a Neumann boundary condition for a single a single
180
     *        equation.
181
182
183
     *
     * \param eqIdx The index of the equation on which the outflow
     *              condition applies.
184
185
186
     */
    void setOutflow(int eqIdx)
    {
187
        resetEq(eqIdx);
188
189
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isOutflow = true;
190
191
192
193
194
195
196
197
198
    }

    /*!
     * \brief Set a boundary condition for a single equation to
     *        a Dirichlet-like coupling condition.
     */
    void setCouplingDirichlet(int eqIdx)
    {
        resetEq(eqIdx);
199
200
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isCouplingDirichlet = true;
201
202
203
204
205
206
207
208
209
    }

    /*!
     * \brief Set a boundary condition for a single equation to
     *        a Neumann-like coupling condition.
     */
    void setCouplingNeumann(int eqIdx)
    {
        resetEq(eqIdx);
210
211
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isCouplingNeumann = true;
212
213
    }

214
    /*!
215
     * \brief Set a Dirichlet boundary condition for a single primary
216
     *        variable.
217
     *
218
219
     * Depending on the discretization, setting the Dirichlet condition
     * will replace the balance equation with index equal to pvIdx.
220
     *
221
222
     * \param pvIdx The index of the primary variable inside a
     *              PrimaryVariables object.
223
     */
224
    void setDirichlet(int pvIdx)
225
    {
226
        setDirichlet(pvIdx, pvIdx);
227
228
229
230
    }

    /*!
     * \brief Returns true if an equation is used to specify a
231
     *        Dirichlet condition.
232
233
     *
     * \param eqIdx The index of the equation
234
235
     */
    bool isDirichlet(unsigned eqIdx) const
Katherina Baber's avatar
Katherina Baber committed
236
    { return boundaryInfo_[eqIdx].isDirichlet; }
237

238
    /*!
Dennis Gläser's avatar
Dennis Gläser committed
239
     * \brief Returns true if all equations are used to specify a
240
241
242
243
244
245
246
247
248
249
     *        Dirichlet condition.
     */
    bool hasOnlyDirichlet() const
    {
        return std::all_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isDirichlet; }
                           );
    }

250
251
    /*!
     * \brief Returns true if some equation is used to specify a
252
     *        Dirichlet condition.
253
254
255
256
257
258
259
     */
    bool hasDirichlet() const
    {
        for (int i = 0; i < numEq; ++i)
            if (boundaryInfo_[i].isDirichlet)
                return true;
        return false;
Katherina Baber's avatar
Katherina Baber committed
260
    }
261
262
263

    /*!
     * \brief Returns true if an equation is used to specify a
264
     *        Neumann condition.
265
266
     *
     * \param eqIdx The index of the equation
267
268
     */
    bool isNeumann(unsigned eqIdx) const
Katherina Baber's avatar
Katherina Baber committed
269
    { return boundaryInfo_[eqIdx].isNeumann; }
270

Dennis Gläser's avatar
Dennis Gläser committed
271
272
273
274
275
276
277
278
279
280
281
282
    /*!
     * \brief Returns true if all equations are used to specify a
     *        Neumann condition.
     */
    bool hasOnlyNeumann() const
    {
        return std::all_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isNeumann; }
                           );
    }

283
284
    /*!
     * \brief Returns true if some equation is used to specify a
285
     *        Neumann condition.
286
287
288
289
     */
    bool hasNeumann() const
    {
        for (int i = 0; i < numEq; ++i)
Andreas Lauser's avatar
Andreas Lauser committed
290
            if (boundaryInfo_[i].isNeumann)
291
292
                return true;
        return false;
Katherina Baber's avatar
Katherina Baber committed
293
    }
294

295
296
297
    /*!
     * \brief Returns true if an equation is used to specify an
     *        outflow condition.
298
299
     *
     * \param eqIdx The index of the equation
300
301
     */
    bool isOutflow(unsigned eqIdx) const
Katherina Baber's avatar
Katherina Baber committed
302
    { return boundaryInfo_[eqIdx].isOutflow; }
303
304
305
306
307
308
309
310
311
312
313

    /*!
     * \brief Returns true if some equation is used to specify an
     *        outflow condition.
     */
    bool hasOutflow() const
    {
        for (int i = 0; i < numEq; ++i)
            if (boundaryInfo_[i].isOutflow)
                return true;
        return false;
Katherina Baber's avatar
Katherina Baber committed
314
    }
315

316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
    /*!
     * \brief Returns true if an equation is used to specify an
     *        Dirichlet coupling condition.
     *
     * \param eqIdx The index of the equation
     */
    bool isCouplingDirichlet(unsigned eqIdx) const
    { return boundaryInfo_[eqIdx].isCouplingDirichlet; }

    /*!
     * \brief Returns true if some equation is used to specify an
     *        Dirichlet coupling condition.
     */
    bool hasCouplingDirichlet() const
    {
        for (int i = 0; i < numEq; ++i)
            if (boundaryInfo_[i].isCouplingDirichlet)
                return true;
        return false;
    }

    /*!
     * \brief Returns true if an equation is used to specify an
     *        Neumann coupling condition.
     *
     * \param eqIdx The index of the equation
     */
    bool isCouplingNeumann(unsigned eqIdx) const
    { return boundaryInfo_[eqIdx].isCouplingNeumann; }

    /*!
     * \brief Returns true if some equation is used to specify an
     *        Neumann coupling condition.
     */
    bool hasCouplingNeumann() const
    {
        for (int i = 0; i < numEq; ++i)
            if (boundaryInfo_[i].isCouplingNeumann)
                return true;
        return false;
    }

    /*!
     * \brief Returns true if an equation is used to specify a
     *        coupling condition.
361
     *
362
     * \param eqIdx The index of the equation
363
364
365
     */
    bool isCoupling(unsigned eqIdx) const
    {
366
        return boundaryInfo_[eqIdx].isCouplingDirichlet
367
               || boundaryInfo_[eqIdx].isCouplingNeumann;
368
369
370
371
372
373
374
375
    }

    /*!
     * \brief Returns true if some equation is used to specify a
     *        coupling condition.
     */
    bool hasCoupling() const
    {
376
377
378
        for (int i = 0; i < numEq; ++i)
            if (isCoupling(i))
                return true;
379
380
381
        return false;
    }

382
383
    /*!
     * \brief Returns the index of the equation which should be used
384
     *        for the Dirichlet condition of the pvIdx's primary
385
     *        variable.
386
387
388
     *
     * \param pvIdx The index of the primary variable which is be set
     *              by the Dirichlet condition.
389
390
     */
    unsigned dirichletToEqIndex(unsigned pvIdx) const
Katherina Baber's avatar
Katherina Baber committed
391
    { return pv2eqIdx_[pvIdx]; }
392
393
394

    /*!
     * \brief Returns the index of the primary variable which should
395
     *        be used for the Dirichlet condition given an equation
396
     *        index.
397
398
399
     *
     * \param eqIdx The index of the equation which is used to set
     *              the Dirichlet condition.
400
401
     */
    unsigned eqToDirichletIndex(unsigned eqIdx) const
Katherina Baber's avatar
Katherina Baber committed
402
    { return eq2pvIdx_[eqIdx]; }
403

404
protected:
405
    //! use bitfields to minimize the size
406
    struct BoundaryInfo {
407
408
409
410
411
412
        bool visited : 1;
        bool isDirichlet : 1;
        bool isNeumann : 1;
        bool isOutflow : 1;
        bool isCouplingDirichlet : 1;
        bool isCouplingNeumann : 1;
413
414
415
416
    };

    std::array<BoundaryInfo, numEq> boundaryInfo_;
    std::array<unsigned int, numEq> eq2pvIdx_, pv2eqIdx_;
417
418
};

419
} // end namespace Dumux
420
421

#endif