boundarytypes.hh 13.2 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
// $Id$
2
3
4
5
6
7
/*****************************************************************************
 *   Copyright (C) 2009 by Andreas Lauser                                    *
 *   Institute of Hydraulic Engineering                                      *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
8
 *   This program is free software: you can redistribute it and/or modify    *
9
 *   it under the terms of the GNU General Public License as published by    *
10
11
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
12
 *                                                                           *
13
14
15
16
17
18
19
 *   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/>.   *
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
 *****************************************************************************/
/*!
 * \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
{

35
36
37
/*!
 * \brief Class to specify the type of a boundary.
 */
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
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
58
            boundaryInfo_[i].isNeumann = 0;
59
            boundaryInfo_[i].isOutflow = 0;
60
61
            boundaryInfo_[i].isCouplingInflow = 0;
            boundaryInfo_[i].isCouplingOutflow = 0;
62
63
64
65
66
67
68

            eq2pvIdx_[i] = i;
            pv2eqIdx_[i] = i;
        };
    }

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

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

    /*!
     * \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
100
            boundaryInfo_[eqIdx].isNeumann = 1;
101
            boundaryInfo_[eqIdx].isOutflow = 0;
102
103
            boundaryInfo_[eqIdx].isCouplingInflow = 0;
            boundaryInfo_[eqIdx].isCouplingOutflow = 0;
104
105
106
107
108
109
110
111
112
113
114
115
116

            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
117
            boundaryInfo_[eqIdx].isNeumann = 0;
118
            boundaryInfo_[eqIdx].isOutflow = 0;
119
120
            boundaryInfo_[eqIdx].isCouplingInflow = 0;
            boundaryInfo_[eqIdx].isCouplingOutflow = 0;
121
122
123
124
125
126
127
128

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

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

129
130
131
132
133
134
135
136
137
138
    /*!
     * \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;
139
140
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
            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;
172
            boundaryInfo_[eqIdx].isOutflow = 0;
173
174
            boundaryInfo_[eqIdx].isCouplingInflow = 0;
            boundaryInfo_[eqIdx].isCouplingOutflow = 1;
175
176
177
178
179

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

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

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

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

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

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

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

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

242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
    /*!
     * \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;
        boundaryInfo_[eqIdx].isCouplingInflow = 1;
        boundaryInfo_[eqIdx].isCouplingOutflow = 0;

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

272
273
    /*!
     * \brief Set a dirichlet boundary condition for a single primary
274
275
276
277
278
     *        variable. 
     *
     * 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!
279
280
281
282
283
284
285
286
287
288
289
290
     *
     * \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.
291
292
     *
     * \param eqIdx The index of the equation
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
     */
    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.
312
313
     *
     * \param eqIdx The index of the equation
314
315
     */
    bool isNeumann(unsigned eqIdx) const
Andreas Lauser's avatar
Andreas Lauser committed
316
    { return boundaryInfo_[eqIdx].isNeumann; };
317
318
319
320
321
322
323
324

    /*!
     * \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
325
            if (boundaryInfo_[i].isNeumann)
326
327
328
329
                return true;
        return false;
    };

330
331
332
    /*!
     * \brief Returns true if an equation is used to specify an
     *        outflow condition.
333
334
     *
     * \param eqIdx The index of the equation
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
     */
    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;
    };

351
352
353
    /*!
     * \brief Returns true if an equation is used to specify an
     *        inflow coupling condition.
354
355
     *
     * \param eqIdx The index of the equation
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
     */
    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.
375
376
     *
     * \param eqIdx The index of the equation
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
     */
    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;
    };

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

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

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

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

}

#endif