ProblemStat.hpp 12.8 KB
Newer Older
1
2
3
4
5
6
7
8
#pragma once

#include <tuple>
#include <string>
#include <memory>
#include <list>
#include <map>

9
#include <dune/istl/matrix.hh>
10
11
#include <dune/grid/common/grid.hh>

12
#include "AdaptInfo.hpp"
13
#include "Basic.hpp"
14
#include "DirichletBC.hpp"
15
#include "FileWriter.hpp"
16
#include "Flag.hpp"
17
#include "Initfile.hpp"
18
19
#include "LinearAlgebra.hpp"
#include "Mesh.hpp"
20
21
#include "Operator.hpp"
#include "ProblemStatBase.hpp"
22
#include "StandardProblemIteration.hpp"
23
24
#include "Timer.hpp"

25
namespace AMDiS 
26
{
27
28
29
30
31
  
  template <class Traits>
  class ProblemStatSeq : public ProblemStatBase
  {
    using Self = ProblemStatSeq;
32
33
    
  public: // typedefs and static constants
34
      
35
36
    using FeSpaces = typename Traits::FeSpaces;
    using Mesh     = typename Traits::Mesh;
37
38
39
40
41
    using MeshView = typename Mesh::LeafGridView;
    
    using Codim0   = typename MeshView::template Codim<0>;
    using Element  = typename Codim0::Entity;

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

    /// Adds an operator to \ref A.
    void addMatrixOperator(OperatorType& op, 
			   int i, int j,
99
100
                           double* factor = NULL, 
                           double* estFactor = NULL);
101
102
103

    
    /// Adds an operator to \ref rhs.
104
105
106
107
    void addVectorOperator(OperatorType& op, 
                           int i,
                           double* factor = NULL, 
                           double* estFactor = NULL);
108
109
110
111
112
113
114
115
116
117
118
    

    /// 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, 
119
120
		       bool createMatrixData = true, 
		       bool storeMatrixData = false) override;
121
122
123
124
125
    
		       
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo, 
				   Flag flag, 
126
127
128
129
130
				   bool asmMatrix = true, 
				   bool asmVector = true) override;
    

    /// Writes output files.
131
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
132
    
133
  public: // get-methods
134
135

    /// Returns nr of components \ref nComponents
136
    static constexpr size_t getNumComponents()
137
138
139
140
    {
      return nComponents;
    }
    
141
    
142
143
144
    /// Return the \ref systemMatrix
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
145
    
146
    
147
148
149
150
151
    /// Return the \ref solution
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }
    
    /// Return the i'th \ref solution component
152
    template <size_t I = 0>
153
    decltype(auto) getSolution(const index_<I> _i = {})
154
    {
155
      return (*solution)[_i];
156
157
    }
    
158
    
159
160
161
    /// Return the \ref rhs
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
162
    
163
    
164
165
    /// Return the \ref linearSolver
    auto getSolver() { return linearSolver; }
166
    
167
168
    /// Return the \ref mesh
    auto getMesh() { return mesh; }
169
    
170
    /// Return the \ref meshView
171
172
    auto getMeshView() { return meshView; }
    
173
174
    /// Return the \ref feSpaces
    auto getFeSpaces() { return feSpaces; }
175
    
176
    
177
    /// Return the I'th \ref feSpaces component
178
    template <size_t I = 0>
179
    FeSpace<I> const& getFeSpace(const index_<I> = {}) const
180
    {
181
      return std::get<I>(*feSpaces);
182
183
    }
    
184
    /// Implementation of \ref ProblemStatBase::getName
185
186
187
188
189
    virtual std::string getName() const override
    {
      return name;
    }
    
190
    /// Return a vector of names of the problem components
191
192
193
194
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
195
    
196
  protected: // initialization methods
197
    
198
199
    void createMesh()
    {
200
201
202
203
204
205
206
      meshName = "";
      Parameters::get(name + "->mesh", meshName);
      AMDIS_TEST_EXIT(!meshName.empty(),
          "No mesh name specified for '" << name << "->mesh'!");
      
      mesh = MeshCreator<Mesh>::create(meshName);
      meshView = std::make_shared<MeshView>(mesh->leafGridView());
207
208
209
210
211
212
      
      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("");
213
214
215
216
    }
    
    void createFeSpaces()
    {
217
      feSpaces = std::make_shared<FeSpaces>(construct_tuple<FeSpaces>(*meshView));
218
219
    }
    
220
221
222
223
    void createMatricesAndVectors()
    {
      systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
      solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);
224
225
226
      
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
      rhs = std::make_shared<SystemVectorType>(*feSpaces, rhsNames);
227
228
229
230
    }
    
    void createSolver()
    {
231
232
233
234
235
236
      using Creator = SolverCreator<typename SystemMatrixType::MultiMatrix,
                                    typename SystemVectorType::MultiVector>;
                                    
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
      linearSolver = Creator::create(solverName, name + "->solver");
237
238
239
240
    }
    
    void createFileWriter()
    {
241
242
243
      filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", 
                                                        *meshView, 
                                                        componentNames);
244
245
    }
    
246
  protected: // sub-methods to assemble DOFMatrix
247
248
249
250
251
    
    template <class Matrix, class Vector>
    void assemble(std::pair<int, int> row_col,
                  Matrix* matrix,
                  Vector* rhs);
252
253
    
    template <class RowView, class ColView>
254
    bool getElementMatrix(RowView const& rowView,
255
256
			  ColView const& colView,
			  ElementMatrix& elementMatrix,
257
258
			  std::list<std::shared_ptr<OperatorType>>& operators,
                          std::list<double*> const& factors);
259
260
    
    template <class RowView>
261
    bool getElementVector(RowView const& rowView,
262
			  ElementVector& elementVector,
263
264
			  std::list<std::shared_ptr<OperatorType>>& operators,
                          std::list<double*> const& factors);
265
    
266
267
268
269
270
271
272
273
274
275
276
277
    
    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);
    
278
279
  private: // some internal methods
    
280
281
282
283
    template <size_t I = 0>
    typename FeSpace<I>::NodeFactory& 
    getNodeFactory(const index_<I> = index_<I>())
    {
284
      return const_cast<typename FeSpace<I>::NodeFactory&>(std::get<I>(*feSpaces).nodeFactory());
285
286
    }
    
287
288
    
  private:
289
290
291
292
293
294
295
296
    /// 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.
297
298
299
    std::shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
    
    /// Name of the mesh
300
301
    std::string meshName;
    
302
    /// A gridView object
303
304
305
306
307
308
309
310
    std::shared_ptr<MeshView> meshView;
    
    /// Pointer to the meshes for the different problem components
    std::vector<Mesh*> componentMeshes;
    
    /// FE spaces of this problem.
    std::shared_ptr<FeSpaces> feSpaces; // eventuell const
    
311
    /// A FileWriter object
312
    std::shared_ptr<FileWriter<Traits>> filewriter;
313
314
    
    /// An object of the linearSolver Interface
315
    std::shared_ptr<LinearSolverType> linearSolver;
316
    
317
    /// A block-matrix that is filled during assembling
318
    std::shared_ptr<SystemMatrixType> systemMatrix;
319
320
    
    /// A block-vector with the solution components
321
    std::shared_ptr<SystemVectorType> solution;
322
323
324
    
    /// A block-vector (load-vector) corresponding to the right.hand side 
    /// of the equation, filled during assembling
325
    std::shared_ptr<SystemVectorType> rhs;
326
    
327
328
329
330
331
332
333
334
335
336
337
338
339
    /// A map (i,j) -> list<DirichleBC> string a boundary conditions for 
    /// each matrix block
    IdxPairList<DirichletBC<WorldVector>> dirichletBc;
    
    /// A map (i,j) -> list<OperatorType> string the differential operators for 
    /// each matrix block
    IdxPairList<OperatorType> matrixOperators;
    std::map<std::pair<int,int>, std::list<double*>> matrixFactors;
    
    /// A map (i) -> list<OperatorType> string the differential operators for 
    /// each rhs block
    IdxList<OperatorType> vectorOperators;
    std::map<int, std::list<double*>> vectorFactors;
340
  };
341
  
342
343
  
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
344
345
346
347
//   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>>;
348
349
#endif
  
350
  namespace Impl
351
352
353
354
355
356
357
358
359
  {
    template <class ProblemStatType>
    struct ProblemStat : public ProblemStatType,
                         public StandardProblemIteration
    {
      using ProblemStatType::getName;

      /// Constructor
      explicit ProblemStat(std::string nameStr)
360
361
        : ProblemStatType(nameStr)
        , StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
362
363
      {}

364
365
366
367
368
369
370
371
372
      /**
       * \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.
       **/
373
374
      virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
      {
375
        for (size_t i = 0; i < ProblemStatType::getNumComponents(); ++i)
376
377
378
379
380
381
382
383
          if (adaptInfo.spaceToleranceReached(i))
            adaptInfo.allowRefinement(false, i);
          else
            adaptInfo.allowRefinement(true, i);

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

387
388
      /// Implementation of \ref ProblemStatBase::buildBeforeCoarsen.
      virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing. */ }
389
      
390
391
      /// Implementation of \ref ProblemStatBase::estimate.
      virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
392
      
393
      /// Implementation of \ref ProblemStatBase::refineMesh.
394
      virtual Flag refineMesh(AdaptInfo& adaptInfo) override { return {0}; }
395

396
      /// Implementation of \ref ProblemStatBase::coarsenMesh.
397
      virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return {0}; }
398
      
399
      /// Implementation of \ref ProblemStatBase::markElements.
400
      virtual Flag markElements(AdaptInfo& adaptInfo) override { return {0}; }
401
    };
402
403
    
  } // end namespace Impl
404
405
  
  template <class Traits>
406
  using ProblemStat = Impl::ProblemStat<ProblemStatSeq<Traits>>;
407
  
408
} // end namespace AMDiS
409
410
411

#include "ProblemStat.inc.hpp"