boxlocalresidual.hh 20.8 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// $Id$
/*****************************************************************************
 *   Copyright (C) 2007 by Peter Bastian                                     *
 *   Institute of Parallel and Distributed System                            *
 *   Department Simulation of Large Systems                                  *
 *   University of Stuttgart, Germany                                        *
 *                                                                           *
 *   Copyright (C) 2008-2010 by Andreas Lauser                               *
 *   Copyright (C) 2007-2009 by Bernd Flemisch                               *
 *   Institute of Hydraulic Engineering                                      *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
 *   This program is free software; you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation; either version 2 of the License, or       *
 *   (at your option) any later version, as long as this copyright notice    *
 *   is included in its original form.                                       *
 *                                                                           *
 *   This program is distributed WITHOUT ANY WARRANTY.                       *
 *****************************************************************************/
/*!
 * \file
 * \brief Caculates the residual of models based on the box scheme element-wise.
 */
#ifndef DUMUX_BOX_LOCAL_RESIDUAL_HH
#define DUMUX_BOX_LOCAL_RESIDUAL_HH

Andreas Lauser's avatar
Andreas Lauser committed
29
30
#include <dune/istl/matrix.hh>
#include <dune/grid/common/geometry.hh>
Andreas Lauser's avatar
Andreas Lauser committed
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

#include <dumux/common/valgrind.hh>

#include "boxproperties.hh"

namespace Dumux
{
/*!
 * \ingroup BoxModel
 * \brief Element-wise caculation of the residual matrix for models
 *        based on the box scheme .
 *
 * \todo Please doc me more!
 */
template<class TypeTag>
class BoxLocalResidual
{
private:
Andreas Lauser's avatar
Andreas Lauser committed
49
    typedef BoxLocalResidual<TypeTag> ThisType;
Andreas Lauser's avatar
Andreas Lauser committed
50
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) Implementation;
51
52
53
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
Andreas Lauser's avatar
Andreas Lauser committed
54
55

    enum {
56
        numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
Andreas Lauser's avatar
Andreas Lauser committed
57

58
59
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
Andreas Lauser's avatar
Andreas Lauser committed
60
61
62
    };

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
63
    typedef typename GridView::Grid::ctype CoordScalar;
Andreas Lauser's avatar
Andreas Lauser committed
64

Andreas Lauser's avatar
Andreas Lauser committed
65
66
    typedef Dune::FieldVector<Scalar, dim>  LocalPosition;
    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
Andreas Lauser's avatar
Andreas Lauser committed
67

68
69
    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
70
71
    typedef typename GridView::template Codim<dim>::Entity Vertex;
    typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
Andreas Lauser's avatar
Andreas Lauser committed
72

73
74
    typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
    typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
Andreas Lauser's avatar
Andreas Lauser committed
75

76
77
78
    typedef typename GridView::IntersectionIterator IntersectionIterator;
    typedef typename Element::Geometry Geometry;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
Andreas Lauser's avatar
Andreas Lauser committed
79
80
81

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
82
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
Andreas Lauser's avatar
Andreas Lauser committed
83
84
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
85
86
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
Andreas Lauser's avatar
Andreas Lauser committed
87
88
89
90

    typedef Dune::FieldMatrix<Scalar, numEq, numEq>  MatrixBlock;
    typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;

91
92
93
    // copying the local residual class is not a good idea
    BoxLocalResidual(const BoxLocalResidual &);

Andreas Lauser's avatar
Andreas Lauser committed
94
95
96
97
98
99
public:
    BoxLocalResidual()
    { }

    ~BoxLocalResidual()
    { }
100
101
102
103
104
105
106
107
108
109
    
    /*!
     * \brief Initialize the local residual.
     *
     * This assumes that all objects of the simulation have been fully
     * allocated but not necessarrily initialized completely.
     *
     * \param prob The representation of the physical problem to be
     *             solved.
     */
Andreas Lauser's avatar
Andreas Lauser committed
110
111
    void init(Problem &prob)
    { problemPtr_ = &prob; }
112

Andreas Lauser's avatar
Andreas Lauser committed
113
114
115
    /*!
     * \brief Compute the local residual, i.e. the deviation of the
     *        equations from zero.
116
117
118
     *
     * \param element The DUNE Codim<0> entity for which the residual
     *                ought to be calculated
Andreas Lauser's avatar
Andreas Lauser committed
119
120
121
122
123
     */
    void eval(const Element &element)
    {
        FVElementGeometry fvGeom;
        fvGeom.update(gridView_(), element);
124
125
126
        ElementVolumeVariables volVarsPrev, volVarsCur;
        volVarsPrev.update(problem_(),
                           element,
Andreas Lauser's avatar
Andreas Lauser committed
127
128
                           fvGeom,
                           true /* oldSol? */);
129
130
        volVarsCur.update(problem_(),
                          element,
Andreas Lauser's avatar
Andreas Lauser committed
131
132
133
134
                          fvGeom,
                          false /* oldSol? */);
        ElementBoundaryTypes bcTypes;
        bcTypes.update(problem_(), element, fvGeom);
135

Andreas Lauser's avatar
Andreas Lauser committed
136
137
138
139
140
141
        // this is pretty much a HACK because the internal state of
        // the problem is not supposed to be changed during the
        // evaluation of the residual. (Reasons: It is a violation of
        // abstraction, makes everything more prone to errors and is
        // not thread save.) The real solution are context objects!
        problem_().updateCouplingParams(element);
142
143

        asImp_().eval(element, fvGeom, volVarsPrev, volVarsCur, bcTypes);
Andreas Lauser's avatar
Andreas Lauser committed
144
145
    }

Andreas Lauser's avatar
Andreas Lauser committed
146
147
148
149
150
    /*!
     * \brief Compute the storage term for the current solution.
     *
     * This can be used to figure out how much of each conservation
     * quantity is inside the element.
151
152
153
     *
     * \param element The DUNE Codim<0> entity for which the storage
     *                term ought to be calculated
Andreas Lauser's avatar
Andreas Lauser committed
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
     */
    void evalStorage(const Element &element)
    {
        FVElementGeometry fvGeom;
        fvGeom.update(gridView_(), element);
        ElementBoundaryTypes bcTypes;
        bcTypes.update(problem_(), element, fvGeom);
        ElementVolumeVariables volVars;
        volVars.update(problem_(), element, fvGeom, false);
        
        residual_.resize(fvGeom.numVertices);
        residual_ = 0;

        elemPtr_ = &element;
        fvElemGeomPtr_ = &fvGeom;
        bcTypesPtr_ = &bcTypes;
        prevVolVarsPtr_ = 0;
        curVolVarsPtr_ = &volVars;
        asImp_().evalStorage_();
    }

Andreas Lauser's avatar
Andreas Lauser committed
175
176
    /*!
     * \brief Compute the flux term for the current solution.
177
178
179
180
181
     *
     * \param element The DUNE Codim<0> entity for which the residual
     *                ought to be calculated
     * \param curVolVars The volume averaged variables for all
     *                   sub-contol volumes of the element
Andreas Lauser's avatar
Andreas Lauser committed
182
     */
183
184
    void evalFluxes(const Element &element,
                    const ElementVolumeVariables &curVolVars)
Andreas Lauser's avatar
Andreas Lauser committed
185
186
187
188
189
    {
        FVElementGeometry fvGeom;
        fvGeom.update(gridView_(), element);
        ElementBoundaryTypes bcTypes;
        bcTypes.update(problem_(), element, fvGeom);
190

Andreas Lauser's avatar
Andreas Lauser committed
191
192
193
194
195
196
        residual_.resize(fvGeom.numVertices);
        residual_ = 0;

        elemPtr_ = &element;
        fvElemGeomPtr_ = &fvGeom;
        bcTypesPtr_ = &bcTypes;
197
198
        prevVolVarsPtr_ = 0;
        curVolVarsPtr_ = &curVolVars;
Andreas Lauser's avatar
Andreas Lauser committed
199
200
201
202
203
204
205
        asImp_().evalFluxes_();
    }


    /*!
     * \brief Compute the local residual, i.e. the deviation of the
     *        equations from zero.
206
207
208
209
     *
     * \param element The DUNE Codim<0> entity for which the residual
     *                ought to be calculated
     * \param fvGeom The finite-volume geometry of the element
210
211
212
     * \param prevVolVars The volume averaged variables for all
     *                   sub-contol volumes of the element at the previous 
     *                   time level
213
     * \param curVolVars The volume averaged variables for all
214
215
     *                   sub-contol volumes of the element at the current
     *                   time level
216
217
     * \param bcTypes The types of the boundary conditions for all 
     *                vertices of the element
Andreas Lauser's avatar
Andreas Lauser committed
218
     */
219
    void eval(const Element &element,
Andreas Lauser's avatar
Andreas Lauser committed
220
              const FVElementGeometry &fvGeom,
221
222
              const ElementVolumeVariables &prevVolVars,
              const ElementVolumeVariables &curVolVars,
Andreas Lauser's avatar
Andreas Lauser committed
223
224
225
              const ElementBoundaryTypes &bcTypes)
    {
        //Valgrind::CheckDefined(fvGeom);
226
227
        Valgrind::CheckDefined(prevVolVars);
        Valgrind::CheckDefined(curVolVars);
Andreas Lauser's avatar
Andreas Lauser committed
228
229
230

#if HAVE_VALGRIND
        for (int i=0; i < fvGeom.numVertices; i++) {
231
232
            Valgrind::CheckDefined(prevVolVars[i]);
            Valgrind::CheckDefined(curVolVars[i]);
Andreas Lauser's avatar
Andreas Lauser committed
233
234
235
236
237
238
        }
#endif // HAVE_VALGRIND

        elemPtr_ = &element;
        fvElemGeomPtr_ = &fvGeom;
        bcTypesPtr_ = &bcTypes;
239
240
241
        prevVolVarsPtr_ = &prevVolVars;
        curVolVarsPtr_ = &curVolVars;

Andreas Lauser's avatar
Andreas Lauser committed
242
243
244
245
246
247
248
249
250
        // reset residual
        int numVerts = fvElemGeom_().numVertices;
        residual_.resize(numVerts);
        residual_ = 0;

        asImp_().evalFluxes_();
        asImp_().evalVolumeTerms_();

        // evaluate the boundary
251
        asImp_().evalBoundary_();
Andreas Lauser's avatar
Andreas Lauser committed
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268

#if HAVE_VALGRIND
        for (int i=0; i < fvElemGeom_().numVertices; i++)
            Valgrind::CheckDefined(residual_[i]);
#endif // HAVE_VALGRIND
    }

    /*!
     * \brief Returns the local residual for a given sub-control
     *        volume of the element.
     */
    const ElementSolutionVector &residual() const
    { return residual_; }

    /*!
     * \brief Returns the local residual for a given sub-control
     *        volume of the element.
269
270
271
     *
     * \param scvIdx The local index of the sub-control volume
     *               (i.e. the element's local vertex index)
Andreas Lauser's avatar
Andreas Lauser committed
272
     */
273
    const PrimaryVariables &residual(int scvIdx) const
Andreas Lauser's avatar
Andreas Lauser committed
274
275
    { return residual_[scvIdx]; }

276
277
278
    /*!
     * \brief Returns the local residual for a given sub-control
     *        volume of the element - as a reference!
279
280
281
     *
     * \param scvIdx The local index of the sub-control volume
     *               (i.e. the element's local vertex index)
282
283
284
285
     */
    PrimaryVariables& residualReference(int scvIdx)
    { return residual_[scvIdx]; }

286
protected:
Andreas Lauser's avatar
Andreas Lauser committed
287
288
    Implementation &asImp_()
    { return *static_cast<Implementation*>(this); }
289

Andreas Lauser's avatar
Andreas Lauser committed
290
291
292
    const Implementation &asImp_() const
    { return *static_cast<const Implementation*>(this); }

293
294
295
    /*!
     * \brief Evaluate all boundary conditions on the current element.
     */
296
297
298
299
300
301
302
    void evalBoundary_()
    {
        if (bcTypes_().hasNeumann())
            asImp_().evalNeumann_();
        if (bcTypes_().hasDirichlet())
            asImp_().evalDirichlet_();
    }
303
304
305
306
307

    /*!
     * \brief Set the values of the Dirichlet boundary control volumes
     *        of the current element.
     */
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
    void evalDirichlet_()
    {
        PrimaryVariables tmp;
        for (int i = 0; i < fvElemGeom_().numVertices; ++i) {
            const BoundaryTypes &bcTypes = bcTypes_(i);
            if (! bcTypes.hasDirichlet())
                continue;

            // ask the problem for the dirichlet values
            const VertexPointer vPtr = elem_().template subEntity<dim>(i);
            asImp_().problem_().dirichlet(tmp, *vPtr);
            
            // set the dirichlet conditions
            for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                if (!bcTypes.isDirichlet(eqIdx))
                    continue;
                int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
                this->residual_[i][eqIdx] =
                    curPrimaryVars_(i)[pvIdx] - tmp[pvIdx];
            };
        };
    }

331
332
333
334
    /*!
     * \brief Add all Neumann boundary conditions to the local
     *        residual.
     */
335
    void evalNeumann_()
Andreas Lauser's avatar
Andreas Lauser committed
336
    {
337
        Dune::GeometryType geoType = elem_().geometry().type();
Andreas Lauser's avatar
Andreas Lauser committed
338
339
        const ReferenceElement &refElem = ReferenceElements::general(geoType);

340
341
        IntersectionIterator isIt = gridView_().ibegin(elem_());
        const IntersectionIterator &endIt = gridView_().iend(elem_());
Andreas Lauser's avatar
Andreas Lauser committed
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
        for (; isIt != endIt; ++isIt)
        {
            // handle only faces on the boundary
            if (!isIt->boundary())
                continue;

            // Assemble the boundary for all vertices of the current
            // face
            int faceIdx = isIt->indexInInside();
            int numFaceVerts = refElem.size(faceIdx, 1, dim);
            for (int faceVertIdx = 0;
                 faceVertIdx < numFaceVerts;
                 ++faceVertIdx)
            {
                int elemVertIdx = refElem.subEntity(faceIdx,
                                                    1,
                                                    faceVertIdx,
                                                    dim);
360
                
Andreas Lauser's avatar
Andreas Lauser committed
361
362
363
                int boundaryFaceIdx =
                    fvElemGeom_().boundaryFaceIndex(faceIdx, faceVertIdx);

364
365
366
367
368
                // add the residual of all vertices of the boundary
                // segment
                evalNeumannSegment_(isIt,
                                    elemVertIdx,
                                    boundaryFaceIdx);
Andreas Lauser's avatar
Andreas Lauser committed
369
370
371
372
            }
        }
    }

373
374
375
376
    /*!
     * \brief Add Neumann boundary conditions for a single sub-control
     *        volume face to the local residual.
     */
377
378
379
    void evalNeumannSegment_(const IntersectionIterator &isIt,
                             int scvIdx,
                             int boundaryFaceIdx)
Andreas Lauser's avatar
Andreas Lauser committed
380
381
    {
        // temporary vector to store the neumann boundary fluxes
382
        PrimaryVariables values(0.0);
Andreas Lauser's avatar
Andreas Lauser committed
383
        const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
384

Andreas Lauser's avatar
Andreas Lauser committed
385
386
387
388
389
390
391
392
393
394
395
396
397
        // deal with neumann boundaries
        if (bcTypes.hasNeumann()) {
            problem_().neumann(values,
                               elem_(),
                               fvElemGeom_(),
                               *isIt,
                               scvIdx,
                               boundaryFaceIdx);
            values *= fvElemGeom_().boundaryFace[boundaryFaceIdx].area;
            Valgrind::CheckDefined(values);
            residual_[scvIdx] += values;
        }
    }
398

399
400
401
402
    /*!
     * \brief Add the flux terms to the local residual of all
     *        sub-control volumes of the current element.
     */
Andreas Lauser's avatar
Andreas Lauser committed
403
404
405
406
407
408
409
410
411
    void evalFluxes_()
    {
        // calculate the mass flux over the faces and subtract
        // it from the local rates
        for (int k = 0; k < fvElemGeom_().numEdges; k++)
        {
            int i = fvElemGeom_().subContVolFace[k].i;
            int j = fvElemGeom_().subContVolFace[k].j;

412
            PrimaryVariables flux;
Andreas Lauser's avatar
Andreas Lauser committed
413
414
415
416
            Valgrind::SetUndefined(flux);
            this->asImp_().computeFlux(flux, k);
            Valgrind::CheckDefined(flux);

417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
            // The balance equation for a finite volume is:
            // 
            // dStorage/dt = Flux + Source
            //
            // Re-arranging this, we get 
            //
            // dStorage/dt - Source - Flux = 0
            //
            // The residual already contains the "dStorage/dt -
            // Source" term, so we have to subtract the flux
            // term. Since the calculated flux goes from sub-control
            // volume i to sub-control volume j, we need to add the
            // flux to j and subtract it from i
            residual_[i] += flux;
            residual_[j] -= flux;
Andreas Lauser's avatar
Andreas Lauser committed
432
433
434
        }
    }

435
436
437
438
    /*!
     * \brief Set the local residual to the storage terms of all
     *        sub-control volumes of the current element.
     */
Andreas Lauser's avatar
Andreas Lauser committed
439
440
441
442
443
444
445
446
447
448
449
    void evalStorage_()
    {
        // calculate the amount of conservation each quantity inside
        // all sub control volumes
        for (int i=0; i < fvElemGeom_().numVertices; i++) {
            residual_[i] = 0.0;
            this->asImp_().computeStorage(residual_[i], i, false /*isOldSol*/);
            residual_[i] *= fvElemGeom_().subContVol[i].volume;
        }
    }

450
451
452
453
454
    /*!
     * \brief Add the change in the storage terms and the source term
     *        to the local residual of all sub-control volumes of the
     *        current element.
     */
Andreas Lauser's avatar
Andreas Lauser committed
455
456
457
458
459
    void evalVolumeTerms_()
    {
        // evaluate the volume terms (storage + source terms)
        for (int i=0; i < fvElemGeom_().numVertices; i++)
        {
460
            PrimaryVariables dStorage_dt(0), tmp(0);
Andreas Lauser's avatar
Andreas Lauser committed
461
462
463
464
465
466
467

            // mass balance within the element. this is the
            // $\frac{m}{\partial t}$ term if using implicit
            // euler as time discretization.
            //
            // TODO (?): we might need a more explicit way for
            // doing the time discretization...
468
            this->asImp_().computeStorage(dStorage_dt, i, false);
Andreas Lauser's avatar
Andreas Lauser committed
469
470
            this->asImp_().computeStorage(tmp, i, true);

471
472
            dStorage_dt -= tmp;
            dStorage_dt *=
Andreas Lauser's avatar
Andreas Lauser committed
473
474
475
                fvElemGeom_().subContVol[i].volume
                /
                problem_().timeManager().timeStepSize();
476
            residual_[i] += dStorage_dt;
Andreas Lauser's avatar
Andreas Lauser committed
477
478

            // subtract the source term from the local rate
479
            PrimaryVariables source;
Andreas Lauser's avatar
Andreas Lauser committed
480
481
            this->asImp_().computeSource(source, i);
            source *= fvElemGeom_().subContVol[i].volume;
482
483
484
485
486
            residual_[i] -= source;
            
            // make sure that only defined quantities where used
            // to calculate the residual.
            Valgrind::CheckDefined(residual_[i]);
Andreas Lauser's avatar
Andreas Lauser committed
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
        }
    }

    /*!
     * \brief Returns a reference to the problem.
     */
    const Problem &problem_() const
    { return *problemPtr_; };

    /*!
     * \brief Returns a reference to the model.
     */
    const Model &model_() const
    { return problem_().model(); };

    /*!
     * \brief Returns a reference to the vertex mapper.
     */
    const VertexMapper &vertexMapper_() const
506
    { return problem_().vertexMapper(); };
507
508
509
510

    /*!
     * \brief Returns a reference to the grid view.
     */
Andreas Lauser's avatar
Andreas Lauser committed
511
512
513
    const GridView &gridView_() const
    { return problem_().gridView(); }

514
515
516
    /*!
     * \brief Returns a reference to the current element.
     */
Andreas Lauser's avatar
Andreas Lauser committed
517
    const Element &elem_() const
518
    {
Andreas Lauser's avatar
Andreas Lauser committed
519
520
521
522
        Valgrind::CheckDefined(elemPtr_);
        return *elemPtr_;
    }

523
524
525
526
    /*!
     * \brief Returns a reference to the current element's finite
     *        volume geometry.
     */
Andreas Lauser's avatar
Andreas Lauser committed
527
    const FVElementGeometry &fvElemGeom_() const
528
    {
Andreas Lauser's avatar
Andreas Lauser committed
529
        Valgrind::CheckDefined(fvElemGeomPtr_);
530
        return *fvElemGeomPtr_;
Andreas Lauser's avatar
Andreas Lauser committed
531
532
    }

533
534
535
536
    /*!
     * \brief Returns a reference to the primary variables of the i'th
     *        sub-control volume of the current element.
     */
537
    const PrimaryVariables &curPrimaryVars_(int i) const
Andreas Lauser's avatar
Andreas Lauser committed
538
    {
539
        return curVolVars_(i).primaryVars();
Andreas Lauser's avatar
Andreas Lauser committed
540
541
    }

542
543
544
545
    /*!
     * \brief Returns a reference to the current volume variables of
     *        all sub-control volumes of the current element.
     */
546
547
548
549
    const ElementVolumeVariables &curVolVars_() const
    {
        Valgrind::CheckDefined(curVolVarsPtr_);
        return *curVolVarsPtr_;
Andreas Lauser's avatar
Andreas Lauser committed
550
    }
551
552
553
554
555

    /*!
     * \brief Returns a reference to the volume variables of the i-th
     *        sub-control volume of the current element.
     */
556
557
558
    const VolumeVariables &curVolVars_(int i) const
    {
        return curVolVars_()[i];
Andreas Lauser's avatar
Andreas Lauser committed
559
560
    }

561
562
563
564
565
    /*!
     * \brief Returns a reference to the previous time step's volume
     *        variables of all sub-control volumes of the current
     *        element.
     */
566
567
568
569
    const ElementVolumeVariables &prevVolVars_() const
    {
        Valgrind::CheckDefined(prevVolVarsPtr_);
        return *prevVolVarsPtr_;
Andreas Lauser's avatar
Andreas Lauser committed
570
    }
571
572
573
574
575
576

    /*!
     * \brief Returns a reference to the previous time step's volume
     *        variables of the i-th sub-control volume of the current
     *        element.
     */
577
578
579
    const VolumeVariables &prevVolVars_(int i) const
    {
        return prevVolVars_()[i];
Andreas Lauser's avatar
Andreas Lauser committed
580
581
    }

582
583
584
585
    /*!
     * \brief Returns a reference to the boundary types of all
     *        sub-control volumes of the current element.
     */
Andreas Lauser's avatar
Andreas Lauser committed
586
    const ElementBoundaryTypes &bcTypes_() const
587
    {
Andreas Lauser's avatar
Andreas Lauser committed
588
589
590
        Valgrind::CheckDefined(bcTypesPtr_);
        return *bcTypesPtr_;
    }
591
592
593
594
595

    /*!
     * \brief Returns a reference to the boundary types of the i-th
     *        sub-control volume of the current element.
     */
Andreas Lauser's avatar
Andreas Lauser committed
596
    const BoundaryTypes &bcTypes_(int i) const
597
    {
Andreas Lauser's avatar
Andreas Lauser committed
598
599
600
601
602
603
604
605
606
607
608
609
610
        return bcTypes_()[i];
    }

protected:
    ElementSolutionVector residual_;

    // The problem we would like to solve
    Problem *problemPtr_;

    const Element *elemPtr_;
    const FVElementGeometry *fvElemGeomPtr_;

    // current and previous secondary variables for the element
611
612
    const ElementVolumeVariables *prevVolVarsPtr_;
    const ElementVolumeVariables *curVolVarsPtr_;
Andreas Lauser's avatar
Andreas Lauser committed
613
614
615
616
617
618

    const ElementBoundaryTypes *bcTypesPtr_;
};
}

#endif