boundarytypes.hh 12.8 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
    { return boundaryInfo_[eqIdx].isDirichlet ||
238
239
             boundaryInfo_[eqIdx].isCouplingDirichlet;
    }
240

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

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

    /*!
     * \brief Returns true if an equation is used to specify a
269
     *        Neumann condition.
270
271
     *
     * \param eqIdx The index of the equation
272
273
     */
    bool isNeumann(unsigned eqIdx) const
274
275
276
277
    {
        return boundaryInfo_[eqIdx].isNeumann ||
               boundaryInfo_[eqIdx].isCouplingNeumann;
    }
278

Dennis Gläser's avatar
Dennis Gläser committed
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(),
287
288
                           [](const BoundaryInfo& b){ return b.isNeumann ||
                                                             b.isCouplingNeumann; }
Dennis Gläser's avatar
Dennis Gläser committed
289
290
291
                           );
    }

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

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

    /*!
     * \brief Returns true if some equation is used to specify an
     *        outflow condition.
     */
    bool hasOutflow() const
    {
320
321
322
323
        return std::any_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isOutflow; }
                           );
Katherina Baber's avatar
Katherina Baber committed
324
    }
325

326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
    /*!
     * \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
    {
341
342
343
344
        return std::any_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isCouplingDirichlet; }
                           );
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
     *        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
    {
362
363
364
365
        return std::any_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isCouplingNeumann; }
                           );
366
367
368
369
370
    }

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

    /*!
     * \brief Returns true if some equation is used to specify a
     *        coupling condition.
     */
    bool hasCoupling() const
    {
386
387
388
        for (int i = 0; i < numEq; ++i)
            if (isCoupling(i))
                return true;
389
390
391
        return false;
    }

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

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

414
protected:
415
    //! use bitfields to minimize the size
416
    struct BoundaryInfo {
417
418
419
420
421
422
        bool visited : 1;
        bool isDirichlet : 1;
        bool isNeumann : 1;
        bool isOutflow : 1;
        bool isCouplingDirichlet : 1;
        bool isCouplingNeumann : 1;
423
424
425
426
    };

    std::array<BoundaryInfo, numEq> boundaryInfo_;
    std::array<unsigned int, numEq> eq2pvIdx_, pv2eqIdx_;
427
428
};

429
} // end namespace Dumux
430
431

#endif