boxassembler.hh 25.8 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
// $Id$
2
3
/*****************************************************************************
 *   Copyright (C) 2009-2010 by Bernd Flemisch                               *
4
 *   Copyright (C) 2010 by Andreas Lauser                                    *
5
6
7
8
 *   Institute of Hydraulic Engineering                                      *
 *   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
 *****************************************************************************/
Bernd Flemisch's avatar
Bernd Flemisch committed
22
23
24
25
26
/*!
 * \file
 *
 * \brief An assembler for the Jacobian matrix based on PDELab.
 */
27
28
#ifndef DUMUX_BOX_ASSEMBLER_HH
#define DUMUX_BOX_ASSEMBLER_HH
29

30
#if HAVE_DUNE_PDELAB
Andreas Lauser's avatar
Andreas Lauser committed
31
#include "pdelabboxlocaloperator.hh"
32
#endif
Andreas Lauser's avatar
Andreas Lauser committed
33
34

namespace Dumux {
35

36
/*!
37
38
 * \brief An assembler for the Jacobian matrix based on PDELab.
 */
39
template<class TypeTag>
Andreas Lauser's avatar
Andreas Lauser committed
40
class BoxAssembler
41
{
42
43
44
45
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
46
47
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
Andreas Lauser's avatar
Andreas Lauser committed
48

49
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
50
51
52

#if HAVE_DUNE_PDELAB
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
Andreas Lauser's avatar
Andreas Lauser committed
53
54
55
56
57
58
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Constraints)) Constraints;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ScalarGridFunctionSpace)) ScalarGridFunctionSpace;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalOperator)) LocalOperator;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridOperatorSpace)) GridOperatorSpace;
59
#endif
Andreas Lauser's avatar
Andreas Lauser committed
60
61
62

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
Andreas Lauser's avatar
Andreas Lauser committed
63
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
64
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
65

Andreas Lauser's avatar
Andreas Lauser committed
66
    enum{dim = GridView::dimension};
67
68
    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
69
    typedef typename GridView::IntersectionIterator IntersectionIterator;
70
71

    typedef typename GridView::template Codim<dim>::Entity Vertex;
72
    typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
73
    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
Andreas Lauser's avatar
Andreas Lauser committed
74
75
76

    typedef SolutionVector Vector;
    typedef JacobianMatrix Matrix;
77
78
    typedef Matrix RepresentationType;

79
80
    enum {
        enablePartialReassemble = GET_PROP_VALUE(TypeTag, PTAG(EnablePartialReassemble)),
Andreas Lauser's avatar
Andreas Lauser committed
81
82
83
        enableJacobianRecycling = GET_PROP_VALUE(TypeTag, PTAG(EnableJacobianRecycling)),

        numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))
84
85
86
87
88
    };

    // copying the jacobian assembler is not a good idea
    BoxAssembler(const BoxAssembler &);

89
public:
90
91
92
93
    /*!
     * \brief The colors of elements and vertices required for partial
     *        Jacobian reassembly.
     */
94
    enum EntityColor {
95
96
97
98
        /*!
         * Vertex/element that needs to be reassembled because some
         * relative error is above the tolerance
         */
99
100
        Red,

101
102
103
104
105
106
107
108
        /*!
         * Vertex/element that needs to be reassembled because a
         * neighboring element/vertex is red
         */
        Yellow,

        /*!
         * Yellow vertex has only non-green neighbor elements.
109
         *
110
111
112
113
114
115
116
117
118
119
120
         * This means that its relative error is below the tolerance,
         * but its defect can be linearized without any additional
         * cost. This is just an "internal" color which is not used
         * ouside of the jacobian assembler.
         */
        Orange,

        /*!
         * Vertex/element that does not need to be reassembled
         */
        Green
121
122
    };

Andreas Lauser's avatar
Andreas Lauser committed
123
    BoxAssembler()
124
    {
125
#if HAVE_DUNE_PDELAB
126
127
128
129
130
131
132
        fem_ = 0;
        cn_ = 0;
        scalarGridFunctionSpace_ = 0;
        gridFunctionSpace_ = 0;
        constraintsTrafo_ = 0;
        localOperator_ = 0;
        gridOperatorSpace_ = 0;
133
134
135
#endif // HAVE_DUNE_PDELAB

        problemPtr_ = 0;
136
        matrix_ = 0;
137

138
139
140
141
142
        // set reassemble tolerance to 0, so that if partial
        // reassembly of the jacobian matrix is disabled, the
        // reassemble tolerance is always smaller than the current
        // relative tolerance
        reassembleTolerance_ = 0.0;
143
144
    }

Andreas Lauser's avatar
Andreas Lauser committed
145
    ~BoxAssembler()
146
147
    {
        delete matrix_;
148
149

#if HAVE_DUNE_PDELAB
150
151
152
153
154
155
156
        delete gridOperatorSpace_;
        delete localOperator_;
        delete constraintsTrafo_;
        delete gridFunctionSpace_;
        delete scalarGridFunctionSpace_;
        delete cn_;
        delete fem_;
157
#endif
158
159
    }

160
161
162
163
164
165
166
167
168
    /*!
     * \brief Initialize the jacobian assembler.
     *
     * At this point we can assume that all objects in the problem and
     * the model have been allocated. We can not assume that they are
     * fully initialized, though.
     *
     * \param problem The problem object
     */
Andreas Lauser's avatar
Andreas Lauser committed
169
170
    void init(Problem& problem)
    {
171
        problemPtr_ = &problem;
172
173
174
175
176

#if !HAVE_DUNE_PDELAB
        // initialize the BCRS matrix
        createMatrix_();
#else
Andreas Lauser's avatar
Andreas Lauser committed
177
        fem_ = new FEM();
178
        //cn_ = new Constraints(*problemPtr_);
Andreas Lauser's avatar
Andreas Lauser committed
179
        cn_ = new Constraints();
180
        scalarGridFunctionSpace_ = new ScalarGridFunctionSpace(problemPtr_->gridView(), *fem_, *cn_);
Andreas Lauser's avatar
Andreas Lauser committed
181
182
        gridFunctionSpace_ = new GridFunctionSpace(*scalarGridFunctionSpace_);
        //cn_->compute_ghosts(*gridFunctionSpace_);
183

Andreas Lauser's avatar
Andreas Lauser committed
184
185
186
187
188
        //typedef BoundaryIndexHelper<TypeTag> BoundaryFunction;
        //BoundaryFunction *bTypes = new BoundaryFunction();
        constraintsTrafo_ = new ConstraintsTrafo();
        //Dune::PDELab::constraints(*bTypes, *gridFunctionSpace_, *constraintsTrafo_, false);

189
190
        // initialize the grid operator spaces
        localOperator_ = new LocalOperator(problemPtr_->model());
191
        gridOperatorSpace_ =
Andreas Lauser's avatar
Andreas Lauser committed
192
193
            new GridOperatorSpace(*gridFunctionSpace_, *constraintsTrafo_,
                                  *gridFunctionSpace_, *constraintsTrafo_, *localOperator_);
194
195
        matrix_ = new Matrix(*gridOperatorSpace_);
#endif
Andreas Lauser's avatar
Andreas Lauser committed
196

197
198
        // initialize the jacobian matrix and the right hand side
        // vector
Andreas Lauser's avatar
Andreas Lauser committed
199
        *matrix_ = 0;
200
        reuseMatrix_ = false;
201

202
203
204
        int numVerts = gridView_().size(dim);
        int numElems = gridView_().size(0);
        residual_.resize(numVerts);
205

Bernd Flemisch's avatar
Bernd Flemisch committed
206
207
208
209
        if (gridView_().comm().size() > 1)
        	totalElems_ = gridView_().comm().sum(numElems);
        else
        	totalElems_ = numElems;
210

211
212
213
        // initialize data needed for partial reassembly
        if (enablePartialReassemble) {
            vertexColor_.resize(numVerts);
214
            vertexDelta_.resize(numVerts);
215
216
            elementColor_.resize(numElems);
        }
Andreas Lauser's avatar
Andreas Lauser committed
217
        reassembleAll();
Andreas Lauser's avatar
Andreas Lauser committed
218
    }
219

220
221
222
223
224
225
    /*!
     * \brief Assemble the local jacobian of the problem.
     *
     * The current state of affairs (esp. the previous and the current
     * solutions) is represented by the model object.
     */
226
227
228
    void assemble()
    {
        resetSystem_();
229

230
231
        greenElems_ = 0;

232
233
234
235
236
237
238
239
240
        ElementIterator elemIt = gridView_().template begin<0>();
        ElementIterator elemEndIt = gridView_().template end<0>();
        for (; elemIt != elemEndIt; ++elemIt) {
            const Element &elem = *elemIt;
            if (elem.partitionType() == Dune::GhostEntity)
                assembleGhostElement_(elem);
            else
                assembleElement_(elem);
        };
241

242
        if (enablePartialReassemble) {
Bernd Flemisch's avatar
Bernd Flemisch committed
243
244
        	if (gridView_().comm().size() > 1)
        		greenElems_ = gridView_().comm().sum(greenElems_);
245

246
            reassembleTolerance_ = nextReassembleTolerance_;
247
248
            // print some information at the end of the iteration
            problem_().newtonController().endIterMsg()
249
                << ", reassembled "
250
                << totalElems_ - greenElems_ << "/" << totalElems_
251
                << " (" << 100*Scalar(totalElems_ - greenElems_)/totalElems_ << "%) elems";
252
        }
Bernd Flemisch's avatar
Bernd Flemisch committed
253

254
        return;
255
    }
256

257
258
259
260
261
262
263
    /*!
     * \brief If Jacobian matrix recycling is enabled, this method
     *        specifies whether the next call to assemble() just
     *        rescales the storage term or does a full reassembly
     *
     * \param yesno If true, only rescale; else do full Jacobian assembly.
     */
264
265
266
267
268
    void setMatrixReuseable(bool yesno = true)
    {
        if (enableJacobianRecycling)
            reuseMatrix_ = yesno;
    }
269

270
271
272
273
274
    /*!
     * \brief If partial Jacobian matrix reassembly is enabled, this
     *        method causes all elements to be reassembled in the next
     *        assemble() call.
     */
275
276
    void reassembleAll()
    {
277
        nextReassembleTolerance_ = 0.0;
278

Andreas Lauser's avatar
Andreas Lauser committed
279
280
281
282
283
        if (enablePartialReassemble) {
            std::fill(vertexColor_.begin(),
                      vertexColor_.end(),
                      Red);
            std::fill(elementColor_.begin(),
284
                      elementColor_.end(),
Andreas Lauser's avatar
Andreas Lauser committed
285
                      Red);
286
            std::fill(vertexDelta_.begin(),
287
                      vertexDelta_.end(),
288
                      0.0);
Andreas Lauser's avatar
Andreas Lauser committed
289
        }
290
    }
291

292
    /*!
293
294
295
     * \brief Returns the relative error below which a vertex is
     *        considered to be "green" if partial Jacobian reassembly
     *        is enabled.
296
297
298
299
     *
     * This returns the _actual_ relative computed seen by
     * computeColors(), not the tolerance which it was given.
     */
300
301
    Scalar reassembleTolerance() const
    { return reassembleTolerance_; }
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325

    /*!
     * \brief Update the distance where the non-linear system was
     *        originally insistently linearized and the point where it
     *        will be linerized the next time.
     *
     * This only has an effect if partial reassemble is enabled.
     */
    void updateDiscrepancy(const SolutionVector &u,
                           const SolutionVector &uDelta)
    {
        if (!enablePartialReassemble)
            return;

        // update the vector with the distances of the current
        // evaluation point used for linearization from the original
        // evaluation point
        for (int i = 0; i < vertexDelta_.size(); ++i) {
            PrimaryVariables uCurrent(u[i]);
            PrimaryVariables uNext(uCurrent);
            uNext -= uDelta[i];

            // we need to add the distance the solution was moved for
            // this vertex
326
            Scalar dist = model_().relativeErrorVertex(i,
327
328
                                                       uCurrent,
                                                       uNext);
329
            vertexDelta_[i] += std::abs(dist);
330
331
332
        }

    }
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350

    /*!
     * \brief Determine the colors of vertices and elements for partial
     *        reassembly given a relative tolerance.
     *
     * The following approach is used:
     *
     * - Set all vertices and elements to 'green'
     * - Mark all vertices as 'red' which exhibit an relative error above
     *   the tolerance
     * - Mark all elements which feature a 'red' vetex as 'red'
     * - Mark all vertices which are not 'red' and are part of a
     *   'red' element as 'yellow'
     * - Mark all elements which are not 'red' and contain a
     *   'yellow' vertex as 'yellow'
     *
     * \param relTol The relative error below which a vertex won't be
     *               reassembled. Note that this specifies the
351
352
353
     *               worst-case relative error between the last
     *               linearization point and the current solution and
     *               _not_ the delta vector of the Newton iteration!
354
     */
Andreas Lauser's avatar
Andreas Lauser committed
355
    void computeColors(Scalar relTol)
356
    {
357
358
359
360
361
362
        if (!enablePartialReassemble)
            return;

        ElementIterator elemIt = gridView_().template begin<0>();
        ElementIterator elemEndIt = gridView_().template end<0>();

363
364
        // mark the red vertices and update the tolerance of the
        // linearization which actually will get achieved
365
        nextReassembleTolerance_ = 0;
Andreas Lauser's avatar
Andreas Lauser committed
366
        for (int i = 0; i < vertexColor_.size(); ++i) {
367
            vertexColor_[i] = Green;
368
            if (vertexDelta_[i] > relTol) {
369
370
371
                // mark vertex as red if discrepancy is larger than
                // the relative tolerance
                vertexColor_[i] = Red;
372
            }
373
            nextReassembleTolerance_ =
374
                std::max(nextReassembleTolerance_, vertexDelta_[i]);
Andreas Lauser's avatar
Andreas Lauser committed
375
376
        };

377
378
        // Mark all red elements
        for (; elemIt != elemEndIt; ++elemIt) {
379
380
381
            // find out whether the current element features a red
            // vertex
            bool isRed = false;
382
383
384
385
            int numVerts = elemIt->template count<dim>();
            for (int i=0; i < numVerts; ++i) {
                int globalI = vertexMapper_().map(*elemIt, i, dim);
                if (vertexColor_[globalI] == Red) {
386
                    isRed = true;
387
388
389
                    break;
                }
            };
390

391
392
            // if yes, the element color is also red, else it is not
            // red, i.e. green for the mean time
393
            int globalElemIdx = elementMapper_().map(*elemIt);
394
395
396
397
            if (isRed)
                elementColor_[globalElemIdx] = Red;
            else
                elementColor_[globalElemIdx] = Green;
398
        }
399

400
        // Mark yellow vertices (as orange for the mean time)
401
402
403
        elemIt = gridView_().template begin<0>();
        for (; elemIt != elemEndIt; ++elemIt) {
            int elemIdx = this->elementMapper_().map(*elemIt);
404
405
            if (elementColor_[elemIdx] != Red)
                continue; // non-red elements do not tint vertices
406
                          // yellow!
407

408
409
410
411
412
413
            int numVerts = elemIt->template count<dim>();
            for (int i=0; i < numVerts; ++i) {
                int globalI = vertexMapper_().map(*elemIt, i, dim);
                // if a vertex is already red, don't recolor it to
                // yellow!
                if (vertexColor_[globalI] != Red)
414
                    vertexColor_[globalI] = Orange;
415
416
417
            };
        }

418
        // Mark yellow elements
419
420
421
422
423
        elemIt = gridView_().template begin<0>();
        for (; elemIt != elemEndIt; ++elemIt) {
            int elemIdx = this->elementMapper_().map(*elemIt);
            if (elementColor_[elemIdx] == Red) {
                continue; // element is red already!
424
425
            }

426
427
            // check whether the element features a yellow
            // (resp. orange at this point) vertex
428
429
430
431
            bool isYellow = false;
            int numVerts = elemIt->template count<dim>();
            for (int i=0; i < numVerts; ++i) {
                int globalI = vertexMapper_().map(*elemIt, i, dim);
432
                if (vertexColor_[globalI] == Orange) {
433
434
                    isYellow = true;
                    break;
435
                }
436
            };
437

438
            if (isYellow)
439
                elementColor_[elemIdx] = Yellow;
440
        }
441
442

        // Demote orange vertices to yellow ones if it has at least
443
        // one green element as a neighbor.
444
445
446
447
448
449
        elemIt = gridView_().template begin<0>();
        for (; elemIt != elemEndIt; ++elemIt) {
            int elemIdx = this->elementMapper_().map(*elemIt);
            if (elementColor_[elemIdx] != Green)
                continue; // yellow and red elements do not make
                          // orange vertices yellow!
450

451
452
453
454
455
456
457
458
459
            int numVerts = elemIt->template count<dim>();
            for (int i=0; i < numVerts; ++i) {
                int globalI = vertexMapper_().map(*elemIt, i, dim);
                // if a vertex is orange, recolor it to yellow!
                if (vertexColor_[globalI] == Orange)
                    vertexColor_[globalI] = Yellow;
            };
        }

460
        // promote the remaining orange vertices to red
461
        for (int i=0; i < vertexColor_.size(); ++i) {
462
            // if a vertex is green or yellow don't do anything!
463
464
            if (vertexColor_[i] == Green || vertexColor_[i] == Yellow)
                continue;
465
466
467
468

            // make sure the vertex is red (this is a no-op vertices
            // which are already red!)
            vertexColor_[i] = Red;
469

470
471
            // set the error of this vertex to 0 because the system
            // will be consistently linearized at this vertex
472
473
            vertexDelta_[i] = 0.0;
        };
474
    };
475

476
477
478
479
480
481
    /*!
     * \brief Returns the reassemble color of a vertex
     *
     * \param element An element which contains the vertex
     * \param vertIdx The local index of the vertex in the element.
     */
482
    int vertexColor(const Element &element, int vertIdx) const
483
484
485
    {
        if (!enablePartialReassemble)
            return Red; // reassemble unconditionally!
486

487
488
489
490
        int globalIdx = vertexMapper_().map(element, vertIdx, dim);
        return vertexColor_[globalIdx];
    }

491
492
493
494
495
    /*!
     * \brief Returns the reassemble color of a vertex
     *
     * \param globalVertIdx The global index of the vertex.
     */
496
    int vertexColor(int globalVertIdx) const
497
498
499
500
501
502
    {
        if (!enablePartialReassemble)
            return Red; // reassemble unconditionally!
        return vertexColor_[globalVertIdx];
    }

503
504
505
506
507
    /*!
     * \brief Returns the Jacobian reassemble color of an element
     *
     * \param element The Codim-0 DUNE entity
     */
508
    int elementColor(const Element &element) const
509
510
511
    {
        if (!enablePartialReassemble)
            return Red; // reassemble unconditionally!
512

513
514
515
        int globalIdx = elementMapper_().map(element);
        return elementColor_[globalIdx];
    }
516

517
518
519
    /*!
     * \brief Returns the Jacobian reassemble color of an element
     *
Felix Bode's avatar
Felix Bode committed
520
     * \param globalElementIdx The global index of the element.
521
     */
522
    int elementColor(int globalElementIdx) const
523
524
525
526
527
    {
        if (!enablePartialReassemble)
            return Red; // reassemble unconditionally!
        return elementColor_[globalElementIdx];
    }
528

529
#if HAVE_DUNE_PDELAB
530
531
532
    /*!
     * \brief Returns a pointer to the PDELab's grid function space.
     */
533
534
535
536
537
    const GridFunctionSpace& gridFunctionSpace() const
    {
        return *gridFunctionSpace_;
    }

538
539
540
541
    /*!
     * \brief Returns a pointer to the PDELab's constraints
     *        transformation.
     */
542
543
544
545
    const ConstraintsTrafo& constraintsTrafo() const
    {
        return *constraintsTrafo_;
    }
546
#endif // HAVE_DUNE_PDELAB
547

548
549
550
    /*!
     * \brief Return constant reference to global Jacobian matrix.
     */
551
552
553
    const Matrix& matrix() const
    { return *matrix_; }

554
555
556
    /*!
     * \brief Return constant reference to global residual vector.
     */
557
558
559
    const SolutionVector& residual() const
    { return residual_; }

560

561
private:
562
563
564
565
566
567
568
569
#if !HAVE_DUNE_PDELAB
    // Construct the BCRS matrix for the global jacobian
    void createMatrix_()
    {
        int nVerts = gridView_().size(dim);

        // allocate raw matrix
        matrix_ = new Matrix(nVerts, nVerts, Matrix::random);
570

571
572
        // find out the global indices of the neighboring vertices of
        // each vertex
573
        typedef std::set<int> NeighborSet;
574
        std::vector<NeighborSet> neighbors(nVerts);
575
576
577
578
579
580
581
        ElementIterator eIt = gridView_().template begin<0>();
        const ElementIterator eEndIt = gridView_().template end<0>();
        for (; eIt != eEndIt; ++eIt) {
            const Element &elem = *eIt;

            // loop over all element vertices
            int n = elem.template count<dim>();
582
            for (int i = 0; i < n - 1; ++i) {
583
                int globalI = vertexMapper_().map(*eIt, i, dim);
584
                for (int j = i + 1; j < n; ++j) {
585
                    int globalJ = vertexMapper_().map(*eIt, j, dim);
586
587
                    // make sure that vertex j is in the neighbor set
                    // of vertex i and vice-versa
588
                    neighbors[globalI].insert(globalJ);
589
                    neighbors[globalJ].insert(globalI);
590
591
592
                }
            }
        };
593

594
595
596
        // make vertices neighbors to themselfs
        for (int i = 0; i < nVerts; ++i)
            neighbors[i].insert(i);
597

598
        // allocate space for the rows of the matrix
599
600
601
602
603
        for (int i = 0; i < nVerts; ++i) {
            matrix_->setrowsize(i, neighbors[i].size());
        }
        matrix_->endrowsizes();

604
605
606
        // fill the rows with indices. each vertex talks to all of its
        // neighbors. (it also talks to itself since vertices are
        // sometimes quite egocentric.)
607
608
609
610
611
612
613
614
615
616
617
        for (int i = 0; i < nVerts; ++i) {
            typename NeighborSet::iterator nIt = neighbors[i].begin();
            typename NeighborSet::iterator nEndIt = neighbors[i].end();
            for (; nIt != nEndIt; ++nIt) {
                matrix_->addindex(i, *nIt);
            }
        }
        matrix_->endindices();
    };
#endif

618
619
620
621
622
623
624
625
626
627
628
629
630
631
    // reset the global linear system of equations. if partial
    // reassemble is enabled, this means that the jacobian matrix must
    // only be erased partially!
    void resetSystem_()
    {
        // always reset the right hand side.
        residual_ = 0.0;

        if (!enablePartialReassemble) {
            // If partial reassembly of the jacobian is not enabled,
            // we can just reset everything!
            (*matrix_) = 0;
            return;
        }
632

633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
        // reset all entries corrosponding to a red vertex
        for (int rowIdx = 0; rowIdx < matrix_->N(); ++rowIdx) {
            if (vertexColor_[rowIdx] == Green)
                continue; // the equations for this control volume are
                          // already below the treshold

            // set all entries in the row to 0
            typedef typename JacobianMatrix::ColIterator ColIterator;
            ColIterator colIt = (*matrix_)[rowIdx].begin();
            const ColIterator &colEndIt = (*matrix_)[rowIdx].end();
            for (; colIt != colEndIt; ++colIt) {
                (*colIt) = 0.0;
            }
        };
    }

649
650
651
    // assemble a non-ghost element
    void assembleElement_(const Element &elem)
    {
652
653
        if (enablePartialReassemble) {
            int globalElemIdx = model_().elementMapper().map(elem);
654
655
            if (elementColor_[globalElemIdx] == Green) {
                ++greenElems_;
656
                return;
657
            }
658
659
        }

660
        model_().localJacobian().assemble(elem);
661

662
663
664
        int numVertices = elem.template count<dim>();
        for (int i=0; i < numVertices; ++ i) {
            int globI = vertexMapper_().map(elem, i, dim);
665

666
            // update the right hand side
667
            if (vertexColor(globI) == Green) {
668
                continue;
669
            }
670

671
            residual_[globI] += model_().localJacobian().residual(i);
672

673
            // update the jacobian matrix
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
            for (int j=0; j < numVertices; ++ j) {
                int globJ = vertexMapper_().map(elem, j, dim);
                (*matrix_)[globI][globJ] +=
                    model_().localJacobian().mat(i,j);
            }
        }
    }

    // "assemble" a ghost element
    void assembleGhostElement_(const Element &elem)
    {
        int n = elem.template count<dim>();
        for (int i=0; i < n; ++i) {
            const VertexPointer vp = elem.template subEntity<dim>(i);
            if (vp->partitionType() != Dune::GhostEntity)
                continue; // ignore a ghost cell's non-ghost vertices

            // set main diagonal entries for the vertex
            int vIdx = vertexMapper_().map(*vp);
            typedef typename Matrix::block_type BlockType;
            BlockType &J = (*matrix_)[vIdx][vIdx];
            for (int j = 0; j < BlockType::rows; ++j)
                J[j][j] = 1.0;

            // set residual for the vertex
            residual_[vIdx] = 0;
        }
    }

703

704
705
    Problem &problem_()
    { return *problemPtr_; }
706
707
    const Problem &problem_() const
    { return *problemPtr_; }
708
709
710
711
    const Model &model_() const
    { return problem_().model(); }
    Model &model_()
    { return problem_().model(); }
712
713
714
715
716
717
    const GridView &gridView_() const
    { return problem_().gridView(); }
    const VertexMapper &vertexMapper_() const
    { return problem_().vertexMapper(); }
    const ElementMapper &elementMapper_() const
    { return problem_().elementMapper(); }
718

719
    Problem *problemPtr_;
720

721
    // the jacobian matrix
722
    Matrix *matrix_;
723
724
    // the right-hand side
    SolutionVector residual_;
725

726
    // attributes required for jacobian matrix recycling
727
    bool reuseMatrix_;
728

729
    // attributes required for partial jacobian reassembly
730
731
    std::vector<EntityColor> vertexColor_;
    std::vector<EntityColor> elementColor_;
732
733
734
735
    std::vector<Scalar> vertexDelta_;

    int totalElems_;
    int greenElems_;
736

737
738
    Scalar nextReassembleTolerance_;
    Scalar reassembleTolerance_;
739

740
#if HAVE_DUNE_PDELAB
741
742
743
744
745
746
747
748
    // PDELab stuff
    Constraints *cn_;
    FEM *fem_;
    ScalarGridFunctionSpace *scalarGridFunctionSpace_;
    GridFunctionSpace *gridFunctionSpace_;
    ConstraintsTrafo *constraintsTrafo_;
    LocalOperator *localOperator_;
    GridOperatorSpace *gridOperatorSpace_;
749
#endif
750
751
};

Andreas Lauser's avatar
Andreas Lauser committed
752
} // namespace Dumux
753
754

#endif