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
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
        boundaryInfo_[eqIdx].isDirichlet = true;
202
203
204
205
206
207
208
209
210
    }

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

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

    /*!
     * \brief Returns true if an equation is used to specify a
233
     *        Dirichlet condition.
234
235
     *
     * \param eqIdx The index of the equation
236
237
     */
    bool isDirichlet(unsigned eqIdx) const
Katherina Baber's avatar
Katherina Baber committed
238
    { return boundaryInfo_[eqIdx].isDirichlet; }
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
248
249
250
251
     *        Dirichlet condition.
     */
    bool hasOnlyDirichlet() const
    {
        return std::all_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isDirichlet; }
                           );
    }

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