pdelabboxassembler.hh 13.8 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
// $Id$
2
3
4
5
6
7
8
9
10
11
12
13
14
15
/*****************************************************************************
 *   Copyright (C) 2009-2010 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.                       *
 *****************************************************************************/
Andreas Lauser's avatar
Andreas Lauser committed
16
17
#ifndef DUMUX_PDELAB_BOX_ASSEMBLER_HH
#define DUMUX_PDELAB_BOX_ASSEMBLER_HH
18
19
20
21
22
23
24
25

#include<dune/pdelab/finiteelementmap/p1fem.hh>
#include<dune/pdelab/finiteelementmap/q1fem.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/genericdatahandle.hh>
#include<dune/pdelab/backend/istlvectorbackend.hh>
#include<dune/pdelab/backend/istlmatrixbackend.hh>

Andreas Lauser's avatar
Andreas Lauser committed
26
27
28
29
//#include "pdelabboundarytypes.hh"
#include "pdelabboxlocaloperator.hh"

namespace Dumux {
30

Andreas Lauser's avatar
Andreas Lauser committed
31
namespace PDELab {
32

33
template<class TypeTag>
Andreas Lauser's avatar
Andreas Lauser committed
34
class BoxAssembler
35
{
36
37
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
38
    enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
39
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
40
    enum{dim = GridView::dimension};
41
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
42
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
43
44
    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
45
46
47
48
49
50
51
52

    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;

53
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) LocalJacobian;
Andreas Lauser's avatar
Andreas Lauser committed
54
55
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
56
57
58
59
60
61

    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;

    typedef typename GridView::template Codim<dim>::Entity Vertex;
    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
Andreas Lauser's avatar
Andreas Lauser committed
62
63
64

    typedef SolutionVector Vector;
    typedef JacobianMatrix Matrix;
65
66
    typedef Matrix RepresentationType;

67
68
69
70
71
72
73
74
    enum {
        enablePartialReassemble = GET_PROP_VALUE(TypeTag, PTAG(EnablePartialReassemble)),
        enableJacobianRecycling = GET_PROP_VALUE(TypeTag, PTAG(EnableJacobianRecycling))
    };

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

75
public:
76
77
78
79
80
81
    enum EntityColor {
        Red, //!< entity needs to be reassembled because it's above the tolerance
        Yellow, //!< entity needs to be reassembled because a neighboring element is red
        Green //!< entity does not need to be reassembled
    };

Andreas Lauser's avatar
Andreas Lauser committed
82
    BoxAssembler()
83
    {
84
        problemPtr_ = 0;
85
86
87
88
89
90
91
92
93
94
        fem_ = 0;
        cn_ = 0;
        scalarGridFunctionSpace_ = 0;
        gridFunctionSpace_ = 0;
        constraintsTrafo_ = 0;
        localOperator_ = 0;
        gridOperatorSpace_ = 0;
        matrix_ = 0;
    }

Andreas Lauser's avatar
Andreas Lauser committed
95
    ~BoxAssembler()
96
97
98
99
100
101
102
103
104
105
106
    {
        delete matrix_;
        delete gridOperatorSpace_;
        delete localOperator_;
        delete constraintsTrafo_;
        delete gridFunctionSpace_;
        delete scalarGridFunctionSpace_;
        delete cn_;
        delete fem_;
    }

Andreas Lauser's avatar
Andreas Lauser committed
107
108
    void init(Problem& problem)
    {
109
        problemPtr_ = &problem;
Andreas Lauser's avatar
Andreas Lauser committed
110
        fem_ = new FEM();
111
        //cn_ = new Constraints(*problemPtr_);
Andreas Lauser's avatar
Andreas Lauser committed
112
        cn_ = new Constraints();
113
        scalarGridFunctionSpace_ = new ScalarGridFunctionSpace(problemPtr_->gridView(), *fem_, *cn_);
Andreas Lauser's avatar
Andreas Lauser committed
114
115
116
        gridFunctionSpace_ = new GridFunctionSpace(*scalarGridFunctionSpace_);

        //cn_->compute_ghosts(*gridFunctionSpace_);
117

Andreas Lauser's avatar
Andreas Lauser committed
118
119
120
121
122
        //typedef BoundaryIndexHelper<TypeTag> BoundaryFunction;
        //BoundaryFunction *bTypes = new BoundaryFunction();
        constraintsTrafo_ = new ConstraintsTrafo();
        //Dune::PDELab::constraints(*bTypes, *gridFunctionSpace_, *constraintsTrafo_, false);

123
124
        // initialize the grid operator spaces
        localOperator_ = new LocalOperator(problemPtr_->model());
125
        gridOperatorSpace_ =
Andreas Lauser's avatar
Andreas Lauser committed
126
127
128
            new GridOperatorSpace(*gridFunctionSpace_, *constraintsTrafo_,
                                  *gridFunctionSpace_, *constraintsTrafo_, *localOperator_);

129
130
        // initialize the jacobian matrix and the right hand side
        // vector
Andreas Lauser's avatar
Andreas Lauser committed
131
132
        matrix_ = new Matrix(*gridOperatorSpace_);
        *matrix_ = 0;
133
        reuseMatrix_ = false;
134
135
136
137
138
139
140
141
142
143

        // calculate the ghost vertices
        VertexIterator vIt = gridView_().template begin<dim>();
        VertexIterator vEndIt = gridView_().template end<dim>();
        for (; vIt != vEndIt; ++vIt) {
            if (vIt->partitionType() == Dune::GhostEntity) {
                int vIdx = vertexMapper_().map(*vIt);
                ghostIndices_.push_back(vIdx);
            }
        };
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159

        int numVerts = gridView_().size(dim);
        int numElems = gridView_().size(0);
        residual_.resize(numVerts);
        
        // initialize data needed for partial reassembly
        if (enablePartialReassemble) {
            vertexColor_.resize(numVerts);
            elementColor_.resize(numElems);
        }
        std::fill(vertexColor_.begin(),
                  vertexColor_.end(),
                  Red);
        std::fill(elementColor_.begin(),
                  elementColor_.end(), 
                  Red);
Andreas Lauser's avatar
Andreas Lauser committed
160
161
    }

162
163
    void assemble(SolutionVector &u)
    {
164
        // assemble the global jacobian matrix
165
        if (!reuseMatrix_) {
166
167
            // we actually need to reassemle!
            resetMatrix_();
168
169
170
            gridOperatorSpace_->jacobian(u, *matrix_);
        }
        reuseMatrix_ = false;
171

Bernd Flemisch's avatar
Bernd Flemisch committed
172
173
174
175
        // calculate the global residual
        residual_ = 0;
        gridOperatorSpace_->residual(u, residual_);

176
        typedef typename Matrix::block_type BlockType;
Bernd Flemisch's avatar
Bernd Flemisch committed
177
        // set the entries for the ghost nodes
178
179
180
181
182
183
        BlockType Id(0.0);
        for (int i=0; i < BlockType::rows; ++i)
            Id[i][i] = 1.0;

        for (int i=0; i < ghostIndices_.size(); ++i) {
            int globI = ghostIndices_[i];
Bernd Flemisch's avatar
Bernd Flemisch committed
184
185

            (*matrix_)[globI] = 0;
186
            (*matrix_)[globI][globI] = Id;
Bernd Flemisch's avatar
Bernd Flemisch committed
187
188
            residual_[globI] = 0;
            u[globI] = 0;
189
        }
190
    }
191

192
193
194
195
196
    void setMatrixReuseable(bool yesno = true)
    {
        if (enableJacobianRecycling)
            reuseMatrix_ = yesno;
    }
197

198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
    void reassembleAll()
    {
        std::fill(vertexColor_.begin(),
                  vertexColor_.end(),
                  Red);
        std::fill(elementColor_.begin(),
                  elementColor_.end(), 
                  Red);
    }
    
    void markVertexRed(int globalVertIdx, bool yesno)
    {
        if (enablePartialReassemble)
            vertexColor_[globalVertIdx] = yesno?Red:Green;
    }
    
    void computeColors()
    { 

        if (!enablePartialReassemble)
            return;

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

        // Mark all red elements
        for (; elemIt != elemEndIt; ++elemIt) {
            bool needReassemble = false;
            int numVerts = elemIt->template count<dim>();
            for (int i=0; i < numVerts; ++i) {
                int globalI = vertexMapper_().map(*elemIt, i, dim);
                if (vertexColor_[globalI] == Red) {
                    needReassemble = true;
                    break;
                }
            };
            
            int globalElemIdx = elementMapper_().map(*elemIt);
            elementColor_[globalElemIdx] = needReassemble?Red:Green;
237
        }
238
        
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
        // Mark all yellow vertices
        elemIt = gridView_().template begin<0>();
        for (; elemIt != elemEndIt; ++elemIt) {
            int elemIdx = this->elementMapper_().map(*elemIt);
            if (elementColor_[elemIdx] == Green)
                continue; // green elements do not tint vertices
                          // yellow!
            
            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)
                    vertexColor_[globalI] = Yellow;
            };
        }

        // Mark all yellow elements
        int numGreen = 0;
        elemIt = gridView_().template begin<0>();
        for (; elemIt != elemEndIt; ++elemIt) {
            int elemIdx = this->elementMapper_().map(*elemIt);
            if (elementColor_[elemIdx] == Red) {
                continue; // element is red already!
264
265
            }

266
267
268
269
270
271
272
            bool isYellow = false;
            int numVerts = elemIt->template count<dim>();
            for (int i=0; i < numVerts; ++i) {
                int globalI = vertexMapper_().map(*elemIt, i, dim);
                if (vertexColor_[globalI] == Yellow) {
                    isYellow = true;
                    break;
273
                }
274
275
276
277
            };
            
            if (isYellow) {
                elementColor_[elemIdx] = Yellow;
278
            }
279
            else // elementColor_[elemIdx] == Green
280
                ++ numGreen;
281
        }
282
        
283
284
285
286
287
288
        int numTot = gridView_().size(0);
        problem_().newtonController().endIterMsg()
            << ", reassemble " 
            << numTot - numGreen << "/" << numTot
            << " (" << 100*Scalar(numTot - numGreen)/numTot << "%) elems";
    };
289
    
290
    int vertexColor(const Element &element, int vertIdx) const
291
292
293
294
295
296
297
298
    {
        if (!enablePartialReassemble)
            return Red; // reassemble unconditionally!
        
        int globalIdx = vertexMapper_().map(element, vertIdx, dim);
        return vertexColor_[globalIdx];
    }

299
    int vertexColor(int globalVertIdx) const
300
301
302
303
304
305
    {
        if (!enablePartialReassemble)
            return Red; // reassemble unconditionally!
        return vertexColor_[globalVertIdx];
    }

306
    int elementColor(const Element &element) const
307
308
309
310
311
312
313
    {
        if (!enablePartialReassemble)
            return Red; // reassemble unconditionally!
        
        int globalIdx = elementMapper_().map(element);
        return elementColor_[globalIdx];
    }
314

315
    int elementColor(int globalElementIdx) const
316
317
318
319
320
321
    {
        if (!enablePartialReassemble)
            return Red; // reassemble unconditionally!
        return elementColor_[globalElementIdx];
    }
   
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
    const GridFunctionSpace& gridFunctionSpace() const
    {
        return *gridFunctionSpace_;
    }

    const ConstraintsTrafo& constraintsTrafo() const
    {
        return *constraintsTrafo_;
    }

    //! return const reference to matrix
    const Matrix& matrix() const
    { return *matrix_; }

    //! return const reference to residual
    const SolutionVector& residual() const
    { return residual_; }

private:
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
    void resetMatrix_()
    {
        if (!enablePartialReassemble) {
            (*matrix_) = 0;
            return;
        }
      
        // 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
            // reset 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;
            }
        };

        //printSparseMatrix(std::cout, *matrix_, "J", "row");
    }
    
    Problem &problem_()
    { return *problemPtr_; }
367
368
    const Problem &problem_() const
    { return *problemPtr_; }
369
370
371
372
    const Model &model_() const
    { return problem_().model(); }
    Model &model_()
    { return problem_().model(); }
373
374
375
376
377
378
379
380
    const GridView &gridView_() const
    { return problem_().gridView(); }
    const VertexMapper &vertexMapper_() const
    { return problem_().vertexMapper(); }
    const ElementMapper &elementMapper_() const
    { return problem_().elementMapper(); }
    
    Problem *problemPtr_;
381
382
383
384
385
386
387
388
389
    Constraints *cn_;
    FEM *fem_;
    ScalarGridFunctionSpace *scalarGridFunctionSpace_;
    GridFunctionSpace *gridFunctionSpace_;
    ConstraintsTrafo *constraintsTrafo_;
    LocalOperator *localOperator_;
    GridOperatorSpace *gridOperatorSpace_;

    Matrix *matrix_;
390
    bool reuseMatrix_;
391
392
    std::vector<EntityColor> vertexColor_;
    std::vector<EntityColor> elementColor_;
393
394
    std::vector<int> ghostIndices_;

395
396
397
    SolutionVector residual_;
};

Andreas Lauser's avatar
Andreas Lauser committed
398
399
} // namespace PDELab
} // namespace Dumux
400
401

#endif