ProblemStat.hpp 11.7 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
32
  
  template <class Traits>
  class ProblemStatSeq : public ProblemStatBase
  {
    using Self = ProblemStatSeq;
  public:
33
      
34
35
    using FeSpaces = typename Traits::FeSpaces;
    using Mesh     = typename Traits::Mesh;
36
37
38
39
40
    using MeshView = typename Mesh::LeafGridView;
    
    using Codim0   = typename MeshView::template Codim<0>;
    using Element  = typename Codim0::Entity;

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

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

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

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

    /// Writes output files.
128
    void writeFiles(AdaptInfo& adaptInfo, bool force);
129
    
130
131
    // get-methods
    // -------------------------------------------------------------------------
132
133

    /// Returns nr of components \ref nComponents
134
    static constexpr size_t getNumComponents()
135
136
137
138
    {
      return nComponents;
    }
    
139
140
141
    /// Return the \ref systemMatrix
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
142
    
143
144
145
146
147
    /// Return the \ref solution
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }
    
    /// Return the i'th \ref solution component
148
    template <size_t I = 0>
149
    decltype(auto) getSolution(const index_<I> _i = index_<I>())
150
    {
151
      return (*solution)[_i];
152
153
    }
    
154
155
156
    /// Return the \ref rhs
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
157
    
158
159
    /// Return the \ref linearSolver
    auto getSolver() { return linearSolver; }
160
    
161
162
    /// Return the \ref mesh
    auto getMesh() { return mesh; }
163
    
164
165
    /// Return the \ref feSpaces
    auto getFeSpaces() { return feSpaces; }
166
    
167
    /// Return the I'th \ref feSpaces component
168
    template <size_t I = 0>
169
    FeSpace<I> const& getFeSpace(const index_<I> = index_<I>()) const
170
    {
171
      return std::get<I>(*feSpaces);
172
173
    }
    
174
    /// Return the \ref name of the problem. Implementation of \ref ProblemStatBase::getName
175
176
177
178
179
    virtual std::string getName() const override
    {
      return name;
    }
    
180
181
182
183
184
  protected:
    
    // initialization methods
    // -------------------------------------------------------------------------
    
185
186
187
188
    // Reads a filename from the init file and creates the mesh object with 
    // respect to this file.
    void createMesh()
    {
189
190
191
192
193
194
195
196
      // TODO: Creation from meshname must be generalized to other meshes than AlbertaGrid.
      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());
197
198
199
200
201
    }
    
    //
    void createFeSpaces()
    {
202
      feSpaces = std::make_shared<FeSpaces>(construct_tuple<FeSpaces>(*meshView));
203
204
205
    }
    
    //
206
207
208
209
210
211
212
213
214
215
    void createMatricesAndVectors()
    {
      systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
      solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);
      rhs = std::make_shared<SystemVectorType>(*feSpaces, std::vector<std::string>(nComponents, "rhs"));
    }
    
    //
    void createSolver()
    {
216
217
218
219
220
221
      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");
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
    }
    
    //
    void createFileWriter()
    {
      filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", *meshView, componentNames);
    }
    
    // sub-methods to assemble DOFMatrix
    // -------------------------------------------------------------------------
    
    //
    template <class Matrix, class Vector>
    void assemble(std::pair<int, int> row_col,
                  Matrix* matrix,
                  Vector* rhs);
238
239
240
    
    //
    template <class RowView, class ColView>
241
    bool getElementMatrix(RowView const& rowView,
242
243
244
245
246
247
			  ColView const& colView,
			  ElementMatrix& elementMatrix,
			  std::list<std::shared_ptr<OperatorType>>& operators);
    
    //
    template <class RowView>
248
    bool getElementVector(RowView const& rowView,
249
250
251
			  ElementVector& elementVector,
			  std::list<std::shared_ptr<OperatorType>>& operators);
    
252
253
254
255
256
257
258
259
260
261
262
263
    
    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);
    
264
  private:
265
266
267
268
    template <size_t I = 0>
    typename FeSpace<I>::NodeFactory& 
    getNodeFactory(const index_<I> = index_<I>())
    {
269
      return const_cast<typename FeSpace<I>::NodeFactory&>(std::get<I>(*feSpaces).nodeFactory());
270
271
    }
    
272
273
    
  private:
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
    /// 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.
    // TODO: generalize to multi-mesh problems
    std::shared_ptr<Mesh> mesh;
    std::string meshName;
    
    ///
    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
    
295
    std::shared_ptr<FileWriter<Traits>> filewriter;
296
    std::shared_ptr<LinearSolverType> linearSolver;
297
    
298
    std::shared_ptr<SystemMatrixType> systemMatrix;
299
300
    std::shared_ptr<SystemVectorType> solution;
    std::shared_ptr<SystemVectorType> rhs;
301
302
303
304
    
    IdxPairList<DirichletBC<WorldVector>>	dirichletBc;
    IdxPairList<OperatorType>			matrixOperators;
    IdxList<OperatorType>			vectorOperators;
305
  };
306
  
307
308
  
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
309
310
311
312
//   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>>;
313
314
#endif
  
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
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
367
368
369
370
371
  
  

  namespace detail
  {
    template <class ProblemStatType>
    struct ProblemStat : public ProblemStatType,
                         public StandardProblemIteration
    {
      using ProblemStatType::getName;

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

      /// 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.
      // implements StandardProblemIteration::oneIteration(AdaptInfo&, Flag)
      virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
      {
        for (int i = 0; i < ProblemStatType::getNumComponents(); i++)
          if (adaptInfo.spaceToleranceReached(i))
            adaptInfo.allowRefinement(false, i);
          else
            adaptInfo.allowRefinement(true, i);

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

      /// Implementation of \ref ProblemStatBase::buildBeforeCoarsen().
      virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing here. */ }
      
      /// Implementation of \ref ProblemStatBase::estimate().
      virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing here. */ }
      
      /// Implementation of \ref ProblemStatBase::refineMesh().
      virtual Flag refineMesh(AdaptInfo& adaptInfo) override { /* Does nothing here. */ }

      /// Implementation of \ref ProblemStatBase::coarsenMesh().
      virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { /* Does nothing here. */ }
      
      /// Implementation of \ref ProblemStatBase::markElements().
      virtual Flag markElements(AdaptInfo& adaptInfo) override { /* Does nothing here. */ }
    };
  }
  
  template <class Traits>
  using ProblemStat = detail::ProblemStat<ProblemStatSeq<Traits>>;
  
  
  
372
373
374
375
}

#include "ProblemStat.inc.hpp"