boundarytypes.hh 13.2 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
73
          boundaryInfo_[eqIdx].isDirichlet = false;
          boundaryInfo_[eqIdx].isNeumann = false;
          boundaryInfo_[eqIdx].isOutflow = false;
          boundaryInfo_[eqIdx].isCouplingDirichlet = false;
          boundaryInfo_[eqIdx].isCouplingNeumann = false;
          boundaryInfo_[eqIdx].isCouplingMortar = false;
74
75
76
77
78

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

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

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

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

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

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

128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
    /*!
     * \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);
    }
    /*!
     * \brief Set all boundary conditions to mortar coupling.
     */
    void setAllCouplingMortar()
    {
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
            setCouplingDirichlet(eqIdx);
    }

154
    /*!
155
     * \brief Set a Neumann boundary condition for a single a single
156
     *        equation.
157
158
     *
     * \param eqIdx The index of the equation
159
160
161
     */
    void setNeumann(int eqIdx)
    {
162
        resetEq(eqIdx);
163
164
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isNeumann = true;
165
166
167
    }

    /*!
168
     * \brief Set a Dirichlet boundary condition for a single primary
169
170
     *        variable
     *
171
172
173
174
     * \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
175
     */
176
    void setDirichlet(int pvIdx, int eqIdx)
177
    {
178
        resetEq(eqIdx);
179
180
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isDirichlet = true;
181
182
183
184
185
186

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

187
    /*!
188
     * \brief Set a Neumann boundary condition for a single a single
189
     *        equation.
190
191
192
     *
     * \param eqIdx The index of the equation on which the outflow
     *              condition applies.
193
194
195
     */
    void setOutflow(int eqIdx)
    {
196
        resetEq(eqIdx);
197
198
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isOutflow = true;
199
200
201
202
203
204
205
206
207
    }

    /*!
     * \brief Set a boundary condition for a single equation to
     *        a Dirichlet-like coupling condition.
     */
    void setCouplingDirichlet(int eqIdx)
    {
        resetEq(eqIdx);
208
209
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isCouplingDirichlet = true;
210
211
212
213
214
215
216
217
218
    }

    /*!
     * \brief Set a boundary condition for a single equation to
     *        a Neumann-like coupling condition.
     */
    void setCouplingNeumann(int eqIdx)
    {
        resetEq(eqIdx);
219
220
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isCouplingNeumann = true;
221
222
223
224
225
226
227
228
229
    }

    /*!
     * \brief Set a boundary condition for a single equation to
     *        a mortar coupling condition.
     */
    void setCouplingMortar(int eqIdx)
    {
        resetEq(eqIdx);
230
231
        boundaryInfo_[eqIdx].visited = true;
        boundaryInfo_[eqIdx].isCouplingMortar = true;
232
233
    }

234
    /*!
235
     * \brief Set a Dirichlet boundary condition for a single primary
236
     *        variable.
237
     *
238
239
     * Depending on the discretization, setting the Dirichlet condition
     * will replace the balance equation with index equal to pvIdx.
240
     *
241
242
     * \param pvIdx The index of the primary variable inside a
     *              PrimaryVariables object.
243
     */
244
    void setDirichlet(int pvIdx)
245
    {
246
        setDirichlet(pvIdx, pvIdx);
247
248
249
250
    }

    /*!
     * \brief Returns true if an equation is used to specify a
251
     *        Dirichlet condition.
252
253
     *
     * \param eqIdx The index of the equation
254
255
     */
    bool isDirichlet(unsigned eqIdx) const
Katherina Baber's avatar
Katherina Baber committed
256
    { return boundaryInfo_[eqIdx].isDirichlet; }
257

258
    /*!
Dennis Gläser's avatar
Dennis Gläser committed
259
     * \brief Returns true if all equations are used to specify a
260
261
262
263
264
265
266
267
268
269
     *        Dirichlet condition.
     */
    bool hasOnlyDirichlet() const
    {
        return std::all_of(boundaryInfo_.begin(),
                           boundaryInfo_.end(),
                           [](const BoundaryInfo& b){ return b.isDirichlet; }
                           );
    }

270
271
    /*!
     * \brief Returns true if some equation is used to specify a
272
     *        Dirichlet condition.
273
274
275
276
277
278
279
     */
    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
280
    }
281
282
283

    /*!
     * \brief Returns true if an equation is used to specify a
284
     *        Neumann condition.
285
286
     *
     * \param eqIdx The index of the equation
287
288
     */
    bool isNeumann(unsigned eqIdx) const
Katherina Baber's avatar
Katherina Baber committed
289
    { return boundaryInfo_[eqIdx].isNeumann; }
290

Dennis Gläser's avatar
Dennis Gläser committed
291
292
293
294
295
296
297
298
299
300
301
302
    /*!
     * \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; }
                           );
    }

303
304
    /*!
     * \brief Returns true if some equation is used to specify a
305
     *        Neumann condition.
306
307
308
309
     */
    bool hasNeumann() const
    {
        for (int i = 0; i < numEq; ++i)
Andreas Lauser's avatar
Andreas Lauser committed
310
            if (boundaryInfo_[i].isNeumann)
311
312
                return true;
        return false;
Katherina Baber's avatar
Katherina Baber committed
313
    }
314

315
316
317
    /*!
     * \brief Returns true if an equation is used to specify an
     *        outflow condition.
318
319
     *
     * \param eqIdx The index of the equation
320
321
     */
    bool isOutflow(unsigned eqIdx) const
Katherina Baber's avatar
Katherina Baber committed
322
    { return boundaryInfo_[eqIdx].isOutflow; }
323
324
325
326
327
328
329
330
331
332
333

    /*!
     * \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
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
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
    /*!
     * \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;
    }

378
379
    /*!
     * \brief Returns true if an equation is used to specify a
380
     *        Mortar coupling condition.
381
382
     *
     * \param eqIdx The index of the equation
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
     */
    bool isCouplingMortar(unsigned eqIdx) const
    {
        return boundaryInfo_[eqIdx].isCouplingMortar;
    }

    /*!
     * \brief Returns true if some equation is used to specify an
     *        Mortar coupling condition.
     */
    bool hasCouplingMortar() const
    {
        for (int i = 0; i < numEq; ++i)
            if (boundaryInfo_[i].isCouplingMortar)
                return true;
        return false;
    }

    /*!
     * \brief Returns true if an equation is used to specify a
     *        coupling condition.
404
     *
405
     * \param eqIdx The index of the equation
406
407
408
     */
    bool isCoupling(unsigned eqIdx) const
    {
409
410
411
        return boundaryInfo_[eqIdx].isCouplingDirichlet
               || boundaryInfo_[eqIdx].isCouplingNeumann
               || boundaryInfo_[eqIdx].isCouplingMortar;
412
413
414
415
416
417
418
419
    }

    /*!
     * \brief Returns true if some equation is used to specify a
     *        coupling condition.
     */
    bool hasCoupling() const
    {
420
421
422
        for (int i = 0; i < numEq; ++i)
            if (isCoupling(i))
                return true;
423
424
425
        return false;
    }

426
427
    /*!
     * \brief Returns the index of the equation which should be used
428
     *        for the Dirichlet condition of the pvIdx's primary
429
     *        variable.
430
431
432
     *
     * \param pvIdx The index of the primary variable which is be set
     *              by the Dirichlet condition.
433
434
     */
    unsigned dirichletToEqIndex(unsigned pvIdx) const
Katherina Baber's avatar
Katherina Baber committed
435
    { return pv2eqIdx_[pvIdx]; }
436
437
438

    /*!
     * \brief Returns the index of the primary variable which should
439
     *        be used for the Dirichlet condition given an equation
440
     *        index.
441
442
443
     *
     * \param eqIdx The index of the equation which is used to set
     *              the Dirichlet condition.
444
445
     */
    unsigned eqToDirichletIndex(unsigned eqIdx) const
Katherina Baber's avatar
Katherina Baber committed
446
    { return eq2pvIdx_[eqIdx]; }
447

448
protected:
449
    //! use bitfields to minimize the size
450
    struct BoundaryInfo {
451
452
453
454
455
456
457
        bool visited : 1;
        bool isDirichlet : 1;
        bool isNeumann : 1;
        bool isOutflow : 1;
        bool isCouplingDirichlet : 1;
        bool isCouplingNeumann : 1;
        bool isCouplingMortar : 1;
458
459
460
461
    };

    std::array<BoundaryInfo, numEq> boundaryInfo_;
    std::array<unsigned int, numEq> eq2pvIdx_, pv2eqIdx_;
462
463
};

464
} // end namespace Dumux
465
466

#endif