boundarytypes.hh 12 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
namespace Dumux {
31

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

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

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

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

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

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

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

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

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

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

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

126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
    /*!
     * \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);
    }

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

    /*!
158
     * \brief Set a Dirichlet boundary condition for a single primary
159
160
     *        variable
     *
161
162
163
164
     * \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
165
     */
166
    void setDirichlet(int pvIdx, int eqIdx)
167
    {
168
        resetEq(eqIdx);
169
170
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isDirichlet = true;
171
172
173
174
175
176

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

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

    /*!
     * \brief Set a boundary condition for a single equation to
     *        a Dirichlet-like coupling condition.
     */
    void setCouplingDirichlet(int eqIdx)
    {
        resetEq(eqIdx);
198
199
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isCouplingDirichlet = true;
200
        boundaryInfo_[eqIdx].isDirichlet = 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
        boundaryInfo_[eqIdx].isNeumann = true;
213
214
    }

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

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

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

251
252
    /*!
     * \brief Returns true if some equation is used to specify a
253
     *        Dirichlet condition.
254
255
256
257
258
259
260
     */
    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
261
    }
262
263
264

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

Dennis Gläser's avatar
Dennis Gläser committed
272
273
274
275
276
277
278
279
280
281
282
283
    /*!
     * \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; }
                           );
    }

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

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

    /*!
     * \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
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
361
    /*!
     * \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.
362
     *
363
     * \param eqIdx The index of the equation
364
365
366
     */
    bool isCoupling(unsigned eqIdx) const
    {
367
        return boundaryInfo_[eqIdx].isCouplingDirichlet
368
               || boundaryInfo_[eqIdx].isCouplingNeumann;
369
370
371
372
373
374
375
376
    }

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

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

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

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

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

420
} // end namespace Dumux
421
422

#endif