boundarytypes.hh 13.3 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
/*****************************************************************************
 *   Copyright (C) 2009 by Andreas Lauser                                    *
Andreas Lauser's avatar
Andreas Lauser committed
5
 *   Institute for Modelling Hydraulic and Environmental Systems             *
6
7
8
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
9
 *   This program is free software: you can redistribute it and/or modify    *
10
 *   it under the terms of the GNU General Public License as published by    *
11
12
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
13
 *                                                                           *
14
15
16
17
18
19
20
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           *
 *   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/>.   *
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 *****************************************************************************/
/*!
 * \file
 * \brief Class to specify the type of a boundary.
 */
#ifndef BOUNDARY_TYPES_HH
#define BOUNDARY_TYPES_HH



#include <dumux/common/valgrind.hh>

namespace Dumux
{

36
/*!
37
 * \ingroup BC
38
39
 * \brief Class to specify the type of a boundary.
 */
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
template <int numEq>
class BoundaryTypes
{
public:
    BoundaryTypes()
    { reset(); }

    /*!
     * \brief Reset the boundary types.
     *
     * After this method no equations will be disabled and neither
     * neumann nor dirichlet conditions will be evaluated. This
     * corrosponds to a neumann zero boundary.
     */
    void reset()
    {
        for (int i=0; i < numEq; ++i) {
            boundaryInfo_[i].visited = 0;

            boundaryInfo_[i].isDirichlet = 0;
Andreas Lauser's avatar
Andreas Lauser committed
60
            boundaryInfo_[i].isNeumann = 0;
61
            boundaryInfo_[i].isOutflow = 0;
62
63
            boundaryInfo_[i].isCouplingInflow = 0;
            boundaryInfo_[i].isCouplingOutflow = 0;
64
65
66

            eq2pvIdx_[i] = i;
            pv2eqIdx_[i] = i;
67
        }
68
69
70
    }

    /*!
71
72
73
74
     * \brief Returns true if the boundary types for a given equation
     *        has been specified.
     *
     * \param eqIdx The index of the equation
75
76
77
78
79
80
81
     */
    bool isSet(int eqIdx) const
    { return boundaryInfo_[eqIdx].visited; };

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

    /*!
     * \brief Set all boundary conditions to neuman.
     */
    void setAllNeumann()
    {
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            boundaryInfo_[eqIdx].visited = 1;
            boundaryInfo_[eqIdx].isDirichlet = 0;
Andreas Lauser's avatar
Andreas Lauser committed
102
            boundaryInfo_[eqIdx].isNeumann = 1;
103
            boundaryInfo_[eqIdx].isOutflow = 0;
104
105
            boundaryInfo_[eqIdx].isCouplingInflow = 0;
            boundaryInfo_[eqIdx].isCouplingOutflow = 0;
106
107
108
109
110
111
112
113
114
115
116
117
118

            Valgrind::SetDefined(boundaryInfo_[eqIdx]);
        }
    }

    /*!
     * \brief Set all boundary conditions to dirichlet.
     */
    void setAllDirichlet()
    {
        for (int eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
            boundaryInfo_[eqIdx].visited = 1;
            boundaryInfo_[eqIdx].isDirichlet = 1;
Andreas Lauser's avatar
Andreas Lauser committed
119
            boundaryInfo_[eqIdx].isNeumann = 0;
120
            boundaryInfo_[eqIdx].isOutflow = 0;
121
122
            boundaryInfo_[eqIdx].isCouplingInflow = 0;
            boundaryInfo_[eqIdx].isCouplingOutflow = 0;
123
124
125
126
127
128
129
130

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

            Valgrind::SetDefined(boundaryInfo_[eqIdx]);
        }
    }

131
132
133
134
135
136
137
138
139
140
    /*!
     * \brief Set all boundary conditions to neuman.
     */
    void setAllOutflow()
    {
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            boundaryInfo_[eqIdx].visited = 1;
            boundaryInfo_[eqIdx].isDirichlet = 0;
            boundaryInfo_[eqIdx].isNeumann = 0;
            boundaryInfo_[eqIdx].isOutflow = 1;
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
            boundaryInfo_[eqIdx].isCouplingInflow = 0;
            boundaryInfo_[eqIdx].isCouplingOutflow = 0;

            Valgrind::SetDefined(boundaryInfo_[eqIdx]);
        }
    }

    /*!
     * \brief Set all boundary conditions to coupling inflow.
     */
    void setAllCouplingInflow()
    {
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            boundaryInfo_[eqIdx].visited = 1;
            boundaryInfo_[eqIdx].isDirichlet = 0;
            boundaryInfo_[eqIdx].isNeumann = 0;
            boundaryInfo_[eqIdx].isOutflow = 0;
            boundaryInfo_[eqIdx].isCouplingInflow = 1;
            boundaryInfo_[eqIdx].isCouplingOutflow = 0;

            Valgrind::SetDefined(boundaryInfo_[eqIdx]);
        }
    }

    /*!
     * \brief Set all boundary conditions to coupling outflow.
     */
    void setAllCouplingOutflow()
    {
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            boundaryInfo_[eqIdx].visited = 1;
            boundaryInfo_[eqIdx].isDirichlet = 0;
            boundaryInfo_[eqIdx].isNeumann = 0;
174
            boundaryInfo_[eqIdx].isOutflow = 0;
175
176
            boundaryInfo_[eqIdx].isCouplingInflow = 0;
            boundaryInfo_[eqIdx].isCouplingOutflow = 1;
177
178
179
180
181

            Valgrind::SetDefined(boundaryInfo_[eqIdx]);
        }
    }

182
183
184
    /*!
     * \brief Set a neumann boundary condition for a single a single
     *        equation.
185
186
     *
     * \param eqIdx The index of the equation
187
188
189
190
     */
    void setNeumann(int eqIdx)
    {
        boundaryInfo_[eqIdx].visited = 1;
191
192
193
        boundaryInfo_[eqIdx].isDirichlet = 0;
        boundaryInfo_[eqIdx].isNeumann = 1;
        boundaryInfo_[eqIdx].isOutflow = 0;
194
195
        boundaryInfo_[eqIdx].isCouplingInflow = 0;
        boundaryInfo_[eqIdx].isCouplingOutflow = 0;
196
197
198
199
200
201
202
203

        Valgrind::SetDefined(boundaryInfo_[eqIdx]);
    }

    /*!
     * \brief Set a dirichlet boundary condition for a single primary
     *        variable
     *
204
205
206
207
     * \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
208
     */
209
    void setDirichlet(int pvIdx, int eqIdx)
210
211
212
    {
        boundaryInfo_[eqIdx].visited = 1;
        boundaryInfo_[eqIdx].isDirichlet = 1;
Andreas Lauser's avatar
Andreas Lauser committed
213
        boundaryInfo_[eqIdx].isNeumann = 0;
214
        boundaryInfo_[eqIdx].isOutflow = 0;
215
216
        boundaryInfo_[eqIdx].isCouplingInflow = 0;
        boundaryInfo_[eqIdx].isCouplingOutflow = 0;
217
218
219
220
221
222
223
224

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

        Valgrind::SetDefined(boundaryInfo_[eqIdx]);
    }

225
226
227
    /*!
     * \brief Set a neumann boundary condition for a single a single
     *        equation.
228
229
230
     *
     * \param eqIdx The index of the equation on which the outflow
     *              condition applies.
231
232
233
234
235
236
237
     */
    void setOutflow(int eqIdx)
    {
        boundaryInfo_[eqIdx].visited = 1;
        boundaryInfo_[eqIdx].isDirichlet = 0;
        boundaryInfo_[eqIdx].isNeumann = 0;
        boundaryInfo_[eqIdx].isOutflow = 1;
238
239
        boundaryInfo_[eqIdx].isCouplingInflow = 0;
        boundaryInfo_[eqIdx].isCouplingOutflow = 0;
240
241
242
243

        Valgrind::SetDefined(boundaryInfo_[eqIdx]);
    }

244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
    /*!
     * \brief Set a boundary condition for a single equation to coupling inflow.
     */
    void setCouplingInflow(int eqIdx)
    {
        boundaryInfo_[eqIdx].visited = 1;
        boundaryInfo_[eqIdx].isDirichlet = 0;
        boundaryInfo_[eqIdx].isNeumann = 0;
        boundaryInfo_[eqIdx].isOutflow = 0;
        boundaryInfo_[eqIdx].isCouplingInflow = 1;
        boundaryInfo_[eqIdx].isCouplingOutflow = 0;

        Valgrind::SetDefined(boundaryInfo_[eqIdx]);
    }

    /*!
     * \brief Set a boundary condition for a single equation to coupling outflow.
     */
    void setCouplingOutflow(int eqIdx)
    {
        boundaryInfo_[eqIdx].visited = 1;
        boundaryInfo_[eqIdx].isDirichlet = 0;
        boundaryInfo_[eqIdx].isNeumann = 0;
        boundaryInfo_[eqIdx].isOutflow = 0;
Klaus Mosthaf's avatar
Klaus Mosthaf committed
268
269
        boundaryInfo_[eqIdx].isCouplingInflow = 0;
        boundaryInfo_[eqIdx].isCouplingOutflow = 1;
270
271
272
273

        Valgrind::SetDefined(boundaryInfo_[eqIdx]);
    }

274
275
    /*!
     * \brief Set a dirichlet boundary condition for a single primary
276
     *        variable.
277
278
279
280
     *
     * WARNING: This method assumes that the equation with the same
     * index as the primary variable to be set is used to specify the
     * Dirichlet condition. USE WITH _GREAT_ CARE!
281
282
283
284
285
286
287
288
289
290
291
292
     *
     * \param eqIdx The index of the equation which is assumed to be
     *              equal to the index of the primary variable
     */
    void setDirichlet(int eqIdx)
    {
        setDirichlet(eqIdx, eqIdx);
    }

    /*!
     * \brief Returns true if an equation is used to specify a
     *        dirichlet condition.
293
294
     *
     * \param eqIdx The index of the equation
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
     */
    bool isDirichlet(unsigned eqIdx) const
    { return boundaryInfo_[eqIdx].isDirichlet; };

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

    /*!
     * \brief Returns true if an equation is used to specify a
     *        neumann condition.
314
315
     *
     * \param eqIdx The index of the equation
316
317
     */
    bool isNeumann(unsigned eqIdx) const
Andreas Lauser's avatar
Andreas Lauser committed
318
    { return boundaryInfo_[eqIdx].isNeumann; };
319
320
321
322
323
324
325
326

    /*!
     * \brief Returns true if some equation is used to specify a
     *        neumann condition.
     */
    bool hasNeumann() const
    {
        for (int i = 0; i < numEq; ++i)
Andreas Lauser's avatar
Andreas Lauser committed
327
            if (boundaryInfo_[i].isNeumann)
328
329
330
331
                return true;
        return false;
    };

332
333
334
    /*!
     * \brief Returns true if an equation is used to specify an
     *        outflow condition.
335
336
     *
     * \param eqIdx The index of the equation
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
     */
    bool isOutflow(unsigned eqIdx) const
    { return boundaryInfo_[eqIdx].isOutflow; };

    /*!
     * \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;
    };

353
354
355
    /*!
     * \brief Returns true if an equation is used to specify an
     *        inflow coupling condition.
356
357
     *
     * \param eqIdx The index of the equation
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
     */
    bool isCouplingInflow(unsigned eqIdx) const
    { return boundaryInfo_[eqIdx].isCouplingInflow; };

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

    /*!
     * \brief Returns true if an equation is used to specify an
     *        outflow coupling condition.
377
378
     *
     * \param eqIdx The index of the equation
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
     */
    bool isCouplingOutflow(unsigned eqIdx) const
    { return boundaryInfo_[eqIdx].isCouplingOutflow; };

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

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

    /*!
     * \brief Returns the index of the primary variable which should
408
     *        be used for the Dirichlet condition given an equation
409
     *        index.
410
411
412
     *
     * \param eqIdx The index of the equation which is used to set
     *              the Dirichlet condition.
413
414
415
416
417
418
     */
    unsigned eqToDirichletIndex(unsigned eqIdx) const
    { return eq2pvIdx_[eqIdx]; };

private:
    // this is a bitfield structure!
419
    struct __attribute__((__packed__)) {
420
421
        unsigned char visited : 1;
        unsigned char isDirichlet : 1;
Andreas Lauser's avatar
Andreas Lauser committed
422
        unsigned char isNeumann : 1;
423
        unsigned char isOutflow : 1;
424
425
        unsigned char isCouplingInflow : 1;
        unsigned char isCouplingOutflow : 1;
426
427
428
429
430
431
432
433
434
    } boundaryInfo_[numEq];

    unsigned char eq2pvIdx_[numEq];
    unsigned char pv2eqIdx_[numEq];
};

}

#endif