ProblemStat.inc.hpp 11.5 KB
Newer Older
1
2
#pragma once

3
4
#include "AdaptInfo.hpp"
#include "Loops.hpp"
5
6
7
8
#include "Timer.hpp"

namespace AMDiS
{
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
  template <class Traits>
  void ProblemStatSeq<Traits>::initialize(Flag initFlag,
                                          Self* adoptProblem,
                                          Flag adoptFlag) 
  {
    // === create meshes ===
    if (mesh) {
      AMDIS_WARNING("mesh already created");
    }
    else {
      if (initFlag.isSet(CREATE_MESH) ||
          (!adoptFlag.isSet(INIT_MESH) &&
          (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) 
      {
        createMesh();
      }
25

26
27
28
29
30
31
32
33
34
35
      if (adoptProblem &&
          (adoptFlag.isSet(INIT_MESH) ||
          adoptFlag.isSet(INIT_SYSTEM) ||
          adoptFlag.isSet(INIT_FE_SPACE))) 
      {
        mesh = adoptProblem->getMesh();
      }
      
      componentMeshes.resize(nComponents, mesh.get());
    }
36

37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
    if (!mesh)
      AMDIS_WARNING("no mesh created");
    
    
    int globalRefinements = 0;
    Parameters::get(meshName + "->global refinements", globalRefinements);
    mesh->globalRefine(globalRefinements);
    
    // === create fespace ===
    if (feSpaces) {
      AMDIS_WARNING("feSpaces already created");
    }
    else {
      if (initFlag.isSet(INIT_FE_SPACE) ||
          (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) 
      {
        createFeSpaces();
      }
55

56
57
58
59
60
      if (adoptProblem &&
          (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) 
      {
        feSpaces = adoptProblem->getFeSpaces();
      }
61
    }
62
63
64
65
66
67
68
69
70
71

    if (!feSpaces)
      AMDIS_WARNING("no feSpaces created\n");
    
    
    // === create system ===
    if (initFlag.isSet(INIT_SYSTEM))
      createMatricesAndVectors();

    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM))
72
    {
73
74
75
      solution = adoptProblem->getSolution();
      rhs = adoptProblem->getRhs();
      systemMatrix = adoptProblem->getSystemMatrix();
76
77
    }

78
79
    // === create solver ===
    if (linearSolver)
80
    {
81
82
83
84
85
86
87
88
89
90
91
92
      AMDIS_WARNING("solver already created\n");
    }
    else
    {
      if (initFlag.isSet(INIT_SOLVER))
        createSolver();

      if (adoptProblem && adoptFlag.isSet(INIT_SOLVER))
      {
        AMDIS_TEST_EXIT(!linearSolver, "solver already created\n");
        linearSolver = adoptProblem->getSolver();
      }
93
94
    }

95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
    if (!linearSolver)
      AMDIS_WARNING("no solver created\n");

    // === create file writer ===
    if (initFlag.isSet(INIT_FILEWRITER))
      createFileWriter();
    
    solution->compress();
  }


  template <class Traits>
  void ProblemStatSeq<Traits>::addMatrixOperator(OperatorType& op, 
                                                 int i, int j,
                                                 double* factor, 
                                                 double* estFactor)
  {
    // TODO: currently the factors are ignored
    assert( factor == NULL && estFactor == NULL );
    
    matrixOperators[std::make_pair(i,j)].push_back(std::make_shared<OperatorType>(op));
  }

  
  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(OperatorType& op, 
                                                 int i,
                                                 double* factor, 
                                                 double* estFactor)
  {
    // TODO: currently the factors are ignored
    assert( factor == NULL && estFactor == NULL );
    
    vectorOperators[i].push_back(std::make_shared<OperatorType>(op));
  }

131
132
133
  
  
  /// Adds a Dirichlet boundary condition
134
  template <class Traits>
135
    template <class Predicate, class Values>
136
137
138
  void ProblemStatSeq<Traits>::addDirichletBC(Predicate const& predicate, 
                                              int row, int col, 
                                              Values const& values)
139
140
141
142
143
144
145
  {
    static_assert(Dune::Functions::Concept::isCallable<Predicate, WorldVector>(), 
      "Function passed to addDirichletBC for predicate does not model the Callable<WorldVector> concept");
    
    static_assert(Dune::Functions::Concept::isCallable<Values, WorldVector>(), 
      "Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");

146
147
    AMDIS_TEST_EXIT(row >= 0 && row < nComponents, "row number out of range: " << row);
    AMDIS_TEST_EXIT(col >= 0 && col < nComponents, "col number out of range: " << col);
148
    
149
//     boundaryConditionSet = true;
150
151
152
153
154
155
156
157

    using BcType = DirichletBC<WorldVector>;
    std::shared_ptr<BcType> dirichlet = std::make_shared<BcType>(predicate, values);
      
    dirichletBc[std::make_pair(row, col)].push_back( dirichlet );
  }
  
  
158
159
160
161
  template <class Traits>
  void ProblemStatSeq<Traits>::solve(AdaptInfo& adaptInfo, 
                                     bool createMatrixData, 
                                     bool storeMatrixData)
162
163
  {
    Timer t;
164
    solution->compress();
165
166
    linearSolver->solve(adaptInfo, systemMatrix->getMatrix(), 
                        solution->getVector(), rhs->getVector());
167
    
168
    AMDIS_MSG("solution of discrete system needed " << t.elapsed() << " seconds");
169

170
171
//     adaptInfo.setSolverIterations(statistics.iterations);
//     adaptInfo.setSolverResidual(statistics.reduction);
172
173
174
  }
  
  
175
176
177
  template <class Traits>
  void ProblemStatSeq<Traits>::buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag,
                                                 bool asmMatrix, bool asmVector)
178
179
  {
    Timer t;
180
181
182
183
184
185
    
    For<0, nComponents>::loop([this](auto const _r) 
    {
	this->getNodeFactory(_r).initializeIndices();
    });
    
186

187
188
    size_t nnz = 0;
    For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector](auto const _r) 
189
    {
190
191
      AMDIS_MSG(this->getFeSpace(_r).size() << " DOFs for " << "FeSpace[" << _r << "]");
      
192
      if (asmVector) {
193
194
	rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
195
196
      }

197
      For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector, _r](auto const _c) 
198
      {
199
200
	auto const& rowFeSpace = this->getFeSpace(_r);
        auto const& colFeSpace = this->getFeSpace(_c);
201
	
202
203
	int r = int(_r), c = int(_c);
	auto row_col = std::make_pair(r, c);
204
205
	
        // The DOFMatrix which should be assembled
206
207
208
        auto& dofmatrix    = (*systemMatrix)(_r, _c);
        auto& solution_vec = solution->getDOFVector(_c);
        auto& rhs_vec      = rhs->getDOFVector(_r);
209
210
	  
	if (asmMatrix) {
211
	  dofmatrix.init();
212
213

	  // init boundary condition
214
215
216
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
217
218
              for (auto bc : bc_list.second)
                bc->init(c == int(_c), dofmatrix, solution_vec, rhs_vec);
219
220
221
	    }
	  }
	  
222
223
	  // assemble the DOFMatrix block and the corresponding rhs vector, of r==c
	  this->assemble(row_col, &dofmatrix, (_r == _c && asmVector ? &rhs_vec : NULL));
224
225

	  // finish boundary condition
226
227
228
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
229
230
              for (auto bc : bc_list.second)
                bc->finish(c == int(_c), dofmatrix, solution_vec, rhs_vec);
231
232
	    }
	  }
233
          dofmatrix.finish();
234
	  
235
	  nnz += dofmatrix.getNnz();
236
	}
237
238
      });
    });
239
    
240
241
    AMDIS_MSG("fillin of assembled matrix: " << nnz);
    
242
243
244
245
    AMDIS_MSG("buildAfterCoarsen needed " << t.elapsed() << " seconds");
  }
  

246
247
  template <class Traits>
  void ProblemStatSeq<Traits>::writeFiles(AdaptInfo& adaptInfo, bool force)
248
  {
249
    Timer t;    
250
    filewriter->write(adaptInfo.getTime(), *solution);
251
    AMDIS_MSG("writeFiles needed " << t.elapsed() << " seconds");
252
253
254
  }


255
256
257
258
259
  template <class Traits>
    template <class Matrix, class Vector>
  void ProblemStatSeq<Traits>::assemble(std::pair<int, int> row_col,
                                        Matrix* dofmatrix,
                                        Vector* rhs)
260
  {
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
    auto const& rowFeSpace = dofmatrix->getRowFeSpace();
    auto const& colFeSpace = dofmatrix->getColFeSpace();
    
    for (auto op : matrixOperators[row_col])
      op->init(rowFeSpace, colFeSpace);
    for (auto op : vectorOperators[row_col.first])
      op->init(rowFeSpace, colFeSpace);
            
    auto rowLocalView = rowFeSpace.localView();
    auto rowIndexSet = rowFeSpace.localIndexSet();
    
    auto colLocalView = colFeSpace.localView();
    auto colIndexSet = colFeSpace.localIndexSet();
    
    for (const auto& element : elements(*meshView)) {
      rowLocalView.bind(element);
      rowIndexSet.bind(rowLocalView);
278
      
279
280
      colLocalView.bind(element);
      colIndexSet.bind(colLocalView);
281
      
282
283
284
285
286
287
288
289
290
291
292
293
      if (dofmatrix) {
        ElementMatrix elementMatrix;
        bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix, matrixOperators[row_col]);
        if (add)
          addElementMatrix(*dofmatrix, rowIndexSet, colIndexSet, elementMatrix);
      }
      
      if (rhs) {
        ElementVector elementVector;
        bool add = getElementVector(rowLocalView, elementVector, vectorOperators[row_col.first]);
        if (add)
          addElementVector(*rhs, rowIndexSet, elementVector);
294
      }
295
    }
296
297
298
  }
  
  
299
  template <class Traits>
300
    template <class RowView, class ColView>
301
302
303
304
  bool ProblemStatSeq<Traits>::getElementMatrix(RowView const& rowLocalView,
                                                ColView const& colLocalView,
                                                ElementMatrix& elementMatrix,
                                                std::list<std::shared_ptr<OperatorType>>& operators)
305
  {
306
307
    const auto nRows = rowLocalView.tree().finiteElement().size();
    const auto nCols = colLocalView.tree().finiteElement().size();
308

Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
309
    elementMatrix.setSize(nRows, nCols);
310
311
312
    elementMatrix = 0.0; // fills the entire matrix with zeroes

    for (auto op : operators)
313
      op->getElementMatrix(rowLocalView, colLocalView, elementMatrix);
314
315
316
317
318
319
    
    return !operators.empty();
  }
  
  
  
320
  template <class Traits>
321
    template <class RowView>
322
323
324
  bool ProblemStatSeq<Traits>::getElementVector(RowView const& rowLocalView,
                                                ElementVector& elementVector,
                                                std::list<std::shared_ptr<OperatorType>>& operators)
325
  {
326
    const auto nRows = rowLocalView.tree().finiteElement().size();
327
328

    // Set all entries to zero
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
329
    elementVector.resize(nRows);
330
331
332
    elementVector = 0.0;

    for (auto op : operators)
333
      op->getElementVector(rowLocalView, elementVector);
334
335
336
    
    return !operators.empty();
  }
337
  
338
  template <class Traits>
339
    template <class Matrix, class RowIndexSet, class ColIndexSet>
340
341
342
343
  void ProblemStatSeq<Traits>::addElementMatrix(Matrix& dofmatrix,
                                                RowIndexSet const& rowIndexSet,
                                                ColIndexSet const& colIndexSet,
                                                ElementMatrix const& elementMatrix)
344
345
  {
    for (size_t i = 0; i < elementMatrix.N(); ++i) {
346
347
348
349
350
351
352
      // The global index of the i−th vertex of the element
      const auto row = rowIndexSet.index(i);
      for (size_t j = 0; j < elementMatrix.M(); ++j) {
        // The global index of the j−th vertex of the element
        const auto col = colIndexSet.index(j);
        dofmatrix(row, col) += elementMatrix[i][j];
      }
353
354
355
    }
  }
  
356
  template <class Traits>
357
    template <class Vector, class IndexSet>
358
359
360
  void ProblemStatSeq<Traits>::addElementVector(Vector& dofvector,
                                                IndexSet const& indexSet,
                                                ElementVector const& elementVector)
361
362
  {
    for (size_t i = 0; i < elementVector.size(); ++i) {
363
364
365
      // The global index of the i-th vertex
      const auto row = indexSet.index(i);
      dofvector[row] += elementVector[i];
366
367
    }
  }
368
369

} // end namespace AMDiS