ProblemStat.hpp 14.6 KB
Newer Older
1
2
#pragma once

3
4
5
6
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

7
8
9
10
11
12
#include <tuple>
#include <string>
#include <memory>
#include <list>
#include <map>

13
#include <dune/istl/matrix.hh>
14
15
#include <dune/grid/common/grid.hh>

16
17
18
19
20
21
22
23
24
25
#include <dune/amdis/AdaptInfo.hpp>
#include <dune/amdis/CreatorInterface.hpp>
#include <dune/amdis/DirichletBC.hpp>
#include <dune/amdis/FileWriter.hpp>
#include <dune/amdis/Flag.hpp>
#include <dune/amdis/Initfile.hpp>
#include <dune/amdis/LinearAlgebra.hpp>
#include <dune/amdis/Mesh.hpp>
#include <dune/amdis/Operator.hpp>
#include <dune/amdis/ProblemStatBase.hpp>
26
#include <dune/amdis/ProblemStatTraits.hpp>
27
28
29
30
#include <dune/amdis/StandardProblemIteration.hpp>

#include <dune/amdis/common/Timer.hpp>
#include <dune/amdis/common/TupleUtility.hpp>
31
#include <dune/amdis/common/TypeDefs.hpp>
32
#include <dune/amdis/common/Utility.hpp>
33

34
namespace AMDiS 
35
{
36
37
38
39
40
  
  template <class Traits>
  class ProblemStatSeq : public ProblemStatBase
  {
    using Self = ProblemStatSeq;
41
42
    
  public: // typedefs and static constants
43
      
44
45
    using FeSpaces = typename Traits::FeSpaces;
    using Mesh     = typename Traits::Mesh;
46
47
48
49
50
    using MeshView = typename Mesh::LeafGridView;
    
    using Codim0   = typename MeshView::template Codim<0>;
    using Element  = typename Codim0::Entity;

51
52
    using ElementVector = Impl::ElementVector;
    using ElementMatrix = Impl::ElementMatrix;
53
54
55
    
    
    /// Dimension of the mesh
56
    static constexpr int dim = Traits::dim;
57
58
    
    /// Dimension of the world
59
    static constexpr int dow = Traits::dimworld;
60
61
    
    /// Number of problem components
62
    static constexpr int nComponents = Traits::nComponents;
63
    
64
    
65
66
67
68
69
    template <size_t I>
    using FeSpace = std::tuple_element_t<I, FeSpaces>;
    
    using WorldVector = typename Codim0::Geometry::GlobalCoordinate;
    
70
71
    using SystemVectorType = SystemVector<FeSpaces>;
    using SystemMatrixType = SystemMatrix<FeSpaces>;
72
73
    
    using OperatorType = Operator<MeshView>;
74
75
    using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
                                                   typename SystemVectorType::MultiVector>;
76
    
77
  public:
78
79
80
81
    /** 
     * \brief Constructor. Takes the name of the problem that is used to 
     * access values correpsonding to this püroblem in the parameter file.
     * 
82
     * Parameters read by ProblemStatSeq, with name 'PROB'
83
     *   PROB->names: \ref componentNames
84
     **/
85
    ProblemStatSeq(std::string name)
86
87
      : name(name)
    {
88
      Parameters::get(name + "->names", componentNames);
89
      for (size_t i = componentNames.size(); i < nComponents; ++i)
90
        componentNames.push_back("solution[" + std::to_string(i) + "]");
91
92
93
94
    }
    
    
    /**
95
96
     * \brief Initialisation of the problem.
     * 
97
98
99
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
100
    void initialize(Flag initFlag, 
101
		    Self* adoptProblem = NULL, 
102
		    Flag adoptFlag = INIT_NOTHING);
103
104
    

105
106
    /// Adds an operator to \ref A. 
    /** @{ */
107
108
    void addMatrixOperator(OperatorType& op, 
			   int i, int j,
109
110
                           double* factor = NULL, 
                           double* estFactor = NULL);
111

112
    void addMatrixOperator(shared_ptr<OperatorType> op,
113
114
115
116
                           int i, int j,
                           double* factor = NULL, 
                           double* estFactor = NULL);
    
117
118
119
120
121
    void addMatrixOperator(std::map< std::pair<int,int>, 
                                     shared_ptr<OperatorType> > ops);
    void addMatrixOperator(std::map< std::pair<int,int>, 
                                     std::pair<shared_ptr<OperatorType>, bool> > ops);
    /** @} */
122
    
123
124
125
    
    /// Adds an operator to \ref rhs. 
    /** @{ */
126
127
128
129
    void addVectorOperator(OperatorType& op, 
                           int i,
                           double* factor = NULL, 
                           double* estFactor = NULL);
130
    
131
    void addVectorOperator(shared_ptr<OperatorType> op, 
132
133
134
135
                           int i,
                           double* factor = NULL, 
                           double* estFactor = NULL);
    
136
137
138
139
140
    void addVectorOperator(std::map< int, 
                                     shared_ptr<OperatorType> > ops);
    void addVectorOperator(std::map< int, 
                                     std::pair<shared_ptr<OperatorType>, bool> > ops);
    /** @} */
141
142
143
144
145
146
147
148
149
150

    /// Adds a Dirichlet boundary condition
    template <class Predicate, class Values>
    void addDirichletBC(Predicate const& predicate, 
			int row, int col, 
			Values const& values);

    
    /// Implementation of \ref ProblemStatBase::solve
    virtual void solve(AdaptInfo& adaptInfo, 
151
152
		       bool createMatrixData = true, 
		       bool storeMatrixData = false) override;
153
154
155
156
157
    
		       
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo, 
				   Flag flag, 
158
159
160
161
162
				   bool asmMatrix = true, 
				   bool asmVector = true) override;
    

    /// Writes output files.
163
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
164
    
165
  public: // get-methods
166
167

    /// Returns nr of components \ref nComponents
168
    static constexpr int getNumComponents()
169
170
171
172
    {
      return nComponents;
    }
    
173
    
174
175
176
    /// Return the \ref systemMatrix
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
177
    
178
    
179
180
181
182
183
    /// Return the \ref solution
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }
    
    /// Return the i'th \ref solution component
184
    template <size_t I = 0>
185
    decltype(auto) getSolution(const index_<I> _i = {})
186
    {
187
      return (*solution)[_i];
188
189
    }
    
190
    
191
192
193
    /// Return the \ref rhs
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
194
    
195
    
196
197
    /// Return the \ref linearSolver
    auto getSolver() { return linearSolver; }
198
    
199
200
    /// Return the \ref mesh
    auto getMesh() { return mesh; }
201
    
202
    /// Return the \ref meshView
203
204
    auto getMeshView() { return meshView; }
    
205
206
    /// Return the \ref feSpaces
    auto getFeSpaces() { return feSpaces; }
207
    
208
    
209
    /// Return the I'th \ref feSpaces component
210
    template <size_t I = 0>
211
    FeSpace<I> const& getFeSpace(const index_<I> = {}) const
212
    {
213
      return std::get<I>(*feSpaces);
214
215
    }
    
216
    /// Implementation of \ref ProblemStatBase::getName
217
218
219
220
221
    virtual std::string getName() const override
    {
      return name;
    }
    
222
    /// Return a vector of names of the problem components
223
224
225
226
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
227
    
228
  protected: // initialization methods
229
    
230
231
    void createMesh()
    {
232
233
234
235
236
237
      meshName = "";
      Parameters::get(name + "->mesh", meshName);
      AMDIS_TEST_EXIT(!meshName.empty(),
          "No mesh name specified for '" << name << "->mesh'!");
      
      mesh = MeshCreator<Mesh>::create(meshName);
238
      meshView = make_shared<MeshView>(mesh->leafGridView());
239
240
241
242
243
244
      
      AMDIS_MSG("Create mesh:");
      AMDIS_MSG("#elements = "    << mesh->size(0));
      AMDIS_MSG("#faces/edges = " << mesh->size(1));
      AMDIS_MSG("#vertices = "    << mesh->size(dim));
      AMDIS_MSG("");
245
246
247
248
    }
    
    void createFeSpaces()
    {
249
      feSpaces = make_shared<FeSpaces>(constructTuple<FeSpaces>(*meshView));
250
251
    }
    
252
253
    void createMatricesAndVectors()
    {
254
255
      systemMatrix = make_shared<SystemMatrixType>(*feSpaces);
      solution = make_shared<SystemVectorType>(*feSpaces, componentNames);
256
257
      
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
258
      rhs = make_shared<SystemVectorType>(*feSpaces, rhsNames);
259
260
261
    }
    
    void createSolver()
262
    {    
263
264
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
265
      
266
267
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
268
269
        
      linearSolver = solverCreator->create(name + "->solver");
270
271
272
273
    }
    
    void createFileWriter()
    {
274
      filewriter = make_shared<FileWriter<Traits>>(name + "->output", 
275
276
                                                        *meshView, 
                                                        componentNames);
277
278
    }
    
279
  protected: // sub-methods to assemble DOFMatrix
280
281
282
    
    template <class Matrix, class Vector>
    void assemble(std::pair<int, int> row_col,
283
284
                  Matrix& matrix, bool asmMatrix,
                  Vector& rhs,    bool asmVector);
285
286
    
    template <class RowView, class ColView>
287
    bool getElementMatrix(RowView const& rowView,
288
289
			  ColView const& colView,
			  ElementMatrix& elementMatrix,
290
			  std::list<shared_ptr<OperatorType>>& operators,
291
                          std::list<double*> const& factors);
292
293
    
    template <class RowView>
294
    bool getElementVector(RowView const& rowView,
295
			  ElementVector& elementVector,
296
			  std::list<shared_ptr<OperatorType>>& operators,
297
                          std::list<double*> const& factors);
298
    
299
300
301
302
303
304
305
306
307
308
309
310
    
    template <class Matrix, class RowIndexSet, class ColIndexSet>
    void addElementMatrix(Matrix& matrix,
			  RowIndexSet const& rowIndexSet,
			  ColIndexSet const& colIndexSet,
			  ElementMatrix const& elementMatrix);
    
    template <class Vector, class IndexSet>
    void addElementVector(Vector& vector,
			  IndexSet const& indexSet,
			  ElementVector const& elementvector);
    
311
312
  private: // some internal methods
    
313
314
315
316
    template <size_t I = 0>
    typename FeSpace<I>::NodeFactory& 
    getNodeFactory(const index_<I> = index_<I>())
    {
317
      return const_cast<typename FeSpace<I>::NodeFactory&>(std::get<I>(*feSpaces).nodeFactory());
318
319
    }
    
320
321
    
  private:
322
323
324
325
326
327
328
329
    /// Name of this problem.
    std::string name;

    /// Stores the names for all components. Is used for naming the solution
    /// vectors, \ref solution.
    std::vector<std::string> componentNames;
    
    /// Mesh of this problem.
330
    shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
331
332
    
    /// Name of the mesh
333
334
    std::string meshName;
    
335
    /// A gridView object
336
    shared_ptr<MeshView> meshView;
337
338
339
340
341
    
    /// Pointer to the meshes for the different problem components
    std::vector<Mesh*> componentMeshes;
    
    /// FE spaces of this problem.
342
    shared_ptr<FeSpaces> feSpaces; // eventuell const
343
    
344
    /// A FileWriter object
345
    shared_ptr<FileWriter<Traits>> filewriter;
346
347
    
    /// An object of the linearSolver Interface
348
    shared_ptr<LinearSolverType> linearSolver;
349
    
350
    /// A block-matrix that is filled during assembling
351
    shared_ptr<SystemMatrixType> systemMatrix;
352
353
    
    /// A block-vector with the solution components
354
    shared_ptr<SystemVectorType> solution;
355
356
357
    
    /// A block-vector (load-vector) corresponding to the right.hand side 
    /// of the equation, filled during assembling
358
359
360
361
362
363
364
365
366
    shared_ptr<SystemVectorType> rhs;
    
    
    template <class T>
    using MatrixEntries = std::map< std::pair<int,int>, std::list<T> >;
    
    template <class T>
    using VectorEntries = std::map< int, std::list<T> >;
  
367
    
368
369
    /// A map (i,j) -> list<DirichleBC> string a boundary conditions for 
    /// each matrix block
370
    MatrixEntries<shared_ptr<DirichletBC<WorldVector>>> dirichletBc;
371
372
373
    
    /// A map (i,j) -> list<OperatorType> string the differential operators for 
    /// each matrix block
374
375
    MatrixEntries<shared_ptr<OperatorType>> matrixOperators;
    MatrixEntries<double*>                  matrixFactors;
376
377
    std::map< std::pair<int,int>, bool >    matrixAssembled; // if false, do reassemble
    std::map< std::pair<int,int>, bool >    matrixChanging;  // if true, or vectorAssembled false, do reassemble
378
379
380
    
    /// A map (i) -> list<OperatorType> string the differential operators for 
    /// each rhs block
381
382
    VectorEntries<shared_ptr<OperatorType>> vectorOperators;
    VectorEntries<double*>                  vectorFactors;
383
384
    std::map< int, bool >                   vectorAssembled; // if false, do reassemble
    std::map< int, bool >                   vectorChanging;  // if true, or vectorAssembled false, do reassemble
385
  };
386
  
387
388
  
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
389
390
391
392
//   extern template class ProblemStatSeq<ProblemStatTraits<1>>;
  extern template class ProblemStatSeq<ProblemStatTraits<2,2,1>>;
  extern template class ProblemStatSeq<ProblemStatTraits<2,2,2>>;
//   extern template class ProblemStatSeq<ProblemStatTraits<3>>;
393
394
#endif
  
395
  namespace Impl
396
397
398
399
400
401
402
403
404
  {
    template <class ProblemStatType>
    struct ProblemStat : public ProblemStatType,
                         public StandardProblemIteration
    {
      using ProblemStatType::getName;

      /// Constructor
      explicit ProblemStat(std::string nameStr)
405
406
        : ProblemStatType(nameStr)
        , StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
407
408
      {}

409
410
411
412
413
414
415
416
417
      /**
       * \brief Determines the execution order of the single adaption steps. 
       * 
       * If adapt is true, mesh adaption will be performed. This allows to avoid 
       * mesh adaption, e.g. in timestep adaption loops of timestep adaptive 
       * strategies. 
       * 
       * Implementation of \ref StandardProblemIteration::oneIteration.
       **/
418
419
      virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
      {
420
        for (size_t i = 0; i < ProblemStatType::getNumComponents(); ++i)
421
422
423
424
425
426
427
428
          if (adaptInfo.spaceToleranceReached(i))
            adaptInfo.allowRefinement(false, i);
          else
            adaptInfo.allowRefinement(true, i);

        return StandardProblemIteration::oneIteration(adaptInfo, toDo);
      }
      
429
430
      /// Implementation of \ref ProblemStatBase::buildBeforeRefine.
      virtual void buildBeforeRefine(AdaptInfo&, Flag) override { /* Does nothing. */ }
431

432
433
      /// Implementation of \ref ProblemStatBase::buildBeforeCoarsen.
      virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing. */ }
434
      
435
436
      /// Implementation of \ref ProblemStatBase::estimate.
      virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
437
      
438
      /// Implementation of \ref ProblemStatBase::refineMesh.
439
      virtual Flag refineMesh(AdaptInfo& adaptInfo) override { return {0}; }
440

441
      /// Implementation of \ref ProblemStatBase::coarsenMesh.
442
      virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return {0}; }
443
      
444
      /// Implementation of \ref ProblemStatBase::markElements.
445
      virtual Flag markElements(AdaptInfo& adaptInfo) override { return {0}; }
446
    };
447
448
    
  } // end namespace Impl
449
450
  
  template <class Traits>
451
  using ProblemStat = Impl::ProblemStat<ProblemStatSeq<Traits>>;
452
  
453
} // end namespace AMDiS
454
455
456

#include "ProblemStat.inc.hpp"