boundarytypes.hh 12.1 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
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
     *        Dirichlet condition.
     */
    bool hasOnlyDirichlet() const
    {
        return std::all_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
247
248
                           [](const BoundaryInfo& b){ return b.isDirichlet ||
                                                             b.isCouplingDirichlet; }
249
250
251
                           );
    }

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

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

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

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

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

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

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

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

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

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

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

421
} // end namespace Dumux
422
423

#endif