boundarytypes.hh 12.6 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
 *   the Free Software Foundation, either version 3 of the License, or       *
9
 *   (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
237
238
    { return boundaryInfo_[eqIdx].isDirichlet ||
             boundaryInfo_[eqIdx].isCouplingDirichlet; }
239

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

253
254
    /*!
     * \brief Returns true if some equation is used to specify a
255
     *        Dirichlet condition.
256
257
258
     */
    bool hasDirichlet() const
    {
259
260
        return std::any_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
261
262
                           [](const BoundaryInfo& b){ return b.isDirichlet ||
                                                             b.isCouplingDirichlet; }
263
                           );
Katherina Baber's avatar
Katherina Baber committed
264
    }
265
266
267

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

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

287
288
    /*!
     * \brief Returns true if some equation is used to specify a
289
     *        Neumann condition.
290
291
292
     */
    bool hasNeumann() const
    {
293
294
295
296
        return std::any_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isNeumann; }
                           );
Katherina Baber's avatar
Katherina Baber committed
297
    }
298

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

    /*!
     * \brief Returns true if some equation is used to specify an
     *        outflow condition.
     */
    bool hasOutflow() const
    {
314
315
316
317
        return std::any_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isOutflow; }
                           );
Katherina Baber's avatar
Katherina Baber committed
318
    }
319

320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
    /*!
     * \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
    {
335
336
337
338
        return std::any_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isCouplingDirichlet; }
                           );
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
    }

    /*!
     * \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
    {
356
357
358
359
        return std::any_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isCouplingNeumann; }
                           );
360
361
362
363
364
    }

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

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

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

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

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

    std::array<BoundaryInfo, numEq> boundaryInfo_;
    std::array<unsigned int, numEq> eq2pvIdx_, pv2eqIdx_;
421
422
};

423
} // end namespace Dumux
424
425

#endif