ProblemStat.hpp 13.8 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
#include <dune/amdis/AdaptInfo.hpp>
17
#include <dune/amdis/Assembler.hpp>
18
19
20
21
22
23
24
25
26
#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>
27
#include <dune/amdis/ProblemStatTraits.hpp>
28
29
#include <dune/amdis/StandardProblemIteration.hpp>

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

36
namespace AMDiS
37
{
38
  struct BoundaryType { int b; };
39

40
  template <class Traits>
41
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
42
      : public ProblemStatBase
43
      , public StandardProblemIteration
44
  {
45
    using Self = ProblemStat;
46

47
  public: // typedefs and static constants
48

49
50
    using FeSpaces = typename Traits::FeSpaces;
    using Mesh     = typename Traits::Mesh;
51
    using MeshView = typename Mesh::LeafGridView;
52

53
54
55
    using Codim0   = typename MeshView::template Codim<0>;
    using Element  = typename Codim0::Entity;

56

57
    /// Dimension of the mesh
58
    static constexpr int dim = Traits::dim;
59

60
    /// Dimension of the world
61
    static constexpr int dow = Traits::dow;
62

63
    /// Number of problem components
64
    static constexpr int nComponents = Traits::nComponents;
65
66


67
    template <std::size_t I>
68
    using FeSpace = std::tuple_element_t<I, FeSpaces>;
69

70
    using WorldVector = typename Codim0::Geometry::GlobalCoordinate;
71

72
73
    using SystemVectorType = SystemVector<FeSpaces>;
    using SystemMatrixType = SystemMatrix<FeSpaces>;
74

75
76
    using ElementOperator = Operator<MeshView, Element>;
    using IntersectionOperator = Operator<MeshView, typename MeshView::Intersection>;
77

78
79
    using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
                                                   typename SystemVectorType::MultiVector>;
80

81
  public:
82
83
    /**
     * \brief Constructor. Takes the name of the problem that is used to
84
     * access values correpsonding to this püroblem in the parameter file.
85
     *
86
     * Parameters read by ProblemStat, with name 'PROB'
87
     *   PROB->names: \ref componentNames
88
     **/
89
90
91
    ProblemStat(std::string name)
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      , name(name)
92
    {
93
      Parameters::get(name + "->names", componentNames);
94
      for (std::size_t i = componentNames.size(); i < nComponents; ++i)
95
        componentNames.push_back("solution[" + std::to_string(i) + "]");
96
    }
97

98
99
100
101
    /// Constructor taking additionally a reference to a mesh that is used
    /// instead of the default created mesh, \ref ProblemStat
    ProblemStat(std::string name, Mesh& mesh)
      : ProblemStat(name)
102
103
    {
      this->mesh = std::shared_ptr<Mesh>(&mesh, optional_delete(false));
104
105
106
107
      componentMeshes.resize(nComponents, this->mesh.get());

      meshName = "";
      Parameters::get(name + "->mesh", meshName);
108
109
    }

110

111
    /**
112
     * \brief Initialisation of the problem.
113
     *
114
115
116
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
117
    void initialize(Flag initFlag,
118
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
119
                    Flag adoptFlag = INIT_NOTHING);
120

121
122

    /// Adds an operator to \ref A.
123
    /** @{ */
124
125
126
127
    void addMatrixOperator(ElementOperator const& op, int i, int j,
                           double* factor = nullptr, double* estFactor = nullptr);

    void addMatrixOperator(IntersectionOperator const& op, int i, int j,
128
                           double* factor = nullptr, double* estFactor = nullptr);
129

130
    // operator evaluated on the boundary
131
    void addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j,
132
                           double* factor = nullptr, double* estFactor = nullptr);
133

Praetorius, Simon's avatar
Praetorius, Simon committed
134
135
    void addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops);
    void addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator, bool> > ops);
136
    /** @} */
137
138
139


    /// Adds an operator to \ref rhs.
140
    /** @{ */
141
142
143
144
    void addVectorOperator(ElementOperator const& op, int i,
                           double* factor = nullptr, double* estFactor = nullptr);

    void addVectorOperator(IntersectionOperator const& op, int i,
145
                           double* factor = nullptr, double* estFactor = nullptr);
146

147
    // operator evaluated on the boundary
148
    void addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i,
149
                           double* factor = nullptr, double* estFactor = nullptr);
150

Praetorius, Simon's avatar
Praetorius, Simon committed
151
152
    void addVectorOperator(std::map< int, ElementOperator> ops);
    void addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops);
153
    /** @} */
154

155

156
157
    /// Adds a Dirichlet boundary condition
    template <class Predicate, class Values>
158
    void addDirichletBC(Predicate const& predicate,
Praetorius, Simon's avatar
Praetorius, Simon committed
159
160
                        int row, int col,
                        Values const& values);
161

162

163
    /// Implementation of \ref ProblemStatBase::solve
164
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
165
166
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
167
168


169
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
170
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
171
172
173
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
174

175
176

    /// Writes output files.
177
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
178

Praetorius, Simon's avatar
Praetorius, Simon committed
179

180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
  public: // implementation of iteration interface methods

    /**
      * \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.
      **/
    virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
    {
      for (std::size_t i = 0; i < 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. */ }

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

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

    /// Implementation of \ref ProblemStatBase::refineMesh.
    virtual Flag refineMesh(AdaptInfo& adaptInfo) override { return 0; }

    /// Implementation of \ref ProblemStatBase::coarsenMesh.
    virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return 0; }

    /// Implementation of \ref ProblemStatBase::markElements.
    virtual Flag markElements(AdaptInfo& adaptInfo) override { return 0; }


221
  public: // get-methods
222
223

    /// Returns nr of components \ref nComponents
224
    static constexpr int getNumComponents()
225
226
227
    {
      return nComponents;
    }
228
229


230
    /// Returns a pointer to system-matrix, \ref systemMatrix
231
232
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
233
234


235
    /// Return a pointer to the solution, \ref solution
236
237
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }
238

239
    /// Return the i'th \ref solution component
240
241
    template <std::size_t I = 0>
    decltype(auto) getSolution(const index_t<I> _i = {})
242
    {
243
      return (*solution)[_i];
244
    }
245
246


247
    /// Return a point to the rhs system-vector, \ref rhs
248
249
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
250
251


252
253
    /// Return a pointer to the linear solver, \ref linearSolver
    std::shared_ptr<LinearSolverType> getSolver() { return linearSolver; }
254

Praetorius, Simon's avatar
Praetorius, Simon committed
255
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
256
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
257
      linearSolver = solver;
258
259
    }

260
261
    /// Return a pointer to the mesh, \ref mesh
    std::shared_ptr<Mesh> getMesh() { return mesh; }
262

263
264
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
265
    void setMesh(Mesh& mesh_)
266
    {
267
      mesh = std::shared_ptr<Mesh>(&mesh_, optional_delete(false));
268
      std::fill(componentMeshes.begin(), componentMeshes.end(), this->mesh.get());
269
270
271
272
273
274
275

      createFeSpaces();
      createMatricesAndVectors();
      createFileWriter();
    }


276
277
278
279
280
    /// Return the gridView of the leaf-level
    auto leafGridView() { return mesh->leafGridView(); }

    /// Return the gridView of levle `level`
    auto levelGridView(int level) { return mesh->levelGridView(level); }
281

282
283
    /// Return the \ref feSpaces
    auto getFeSpaces() { return feSpaces; }
284
285


286
    /// Return the I'th \ref feSpaces component
287
288
289
290
291
292
293
294
    template <std::size_t I = 0>
    FeSpace<I> const& getFeSpace(const index_t<I> = {}) const
    {
      return std::get<I>(*feSpaces);
    }

    template <std::size_t I = 0>
    FeSpace<I>& getFeSpace(const index_t<I> = {})
295
    {
296
      return std::get<I>(*feSpaces);
297
    }
298

299
    /// Implementation of \ref ProblemStatBase::getName
300
301
302
303
    virtual std::string getName() const override
    {
      return name;
    }
304

305
    /// Return a vector of names of the problem components
306
307
308
309
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
310

311
  protected: // initialization methods
312

313
314
    void createMesh()
    {
315
316
      meshName = "";
      Parameters::get(name + "->mesh", meshName);
317
      test_exit(!meshName.empty(), "No mesh name specified for '", name, "->mesh'!");
318

319
      mesh = MeshCreator<Mesh>::create(meshName);
320
321
322
323
324
325

      msg("Create mesh:");
      msg("#elements = "   , mesh->size(0));
      msg("#faces/edges = ", mesh->size(1));
      msg("#vertices = "   , mesh->size(dim));
      msg("");
326
    }
327

328
329
    void createFeSpaces()
    {
330
      feSpaces = std::make_shared<FeSpaces>(constructTuple<FeSpaces>(leafGridView()));
331
    }
332

333
334
    void createMatricesAndVectors()
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
335
336
      systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
      solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);
337

338
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
Praetorius, Simon's avatar
Praetorius, Simon committed
339
      rhs = std::make_shared<SystemVectorType>(*feSpaces, rhsNames);
340
    }
341

342
    void createSolver()
343
    {
344
345
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
346

347
348
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
349

350
      linearSolver = solverCreator->create(name + "->solver");
351
    }
352

353
354
    void createFileWriter()
    {
355
      filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", leafGridView(), componentNames);
356
    }
357

Praetorius, Simon's avatar
Praetorius, Simon committed
358

359
  private:
360
361
362
363
364
365
    /// 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;
366

367
    /// Mesh of this problem.
Praetorius, Simon's avatar
Praetorius, Simon committed
368
    std::shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
369

370
    /// Name of the mesh
371
    std::string meshName = "none";
372

373
374
    /// Pointer to the meshes for the different problem components
    std::vector<Mesh*> componentMeshes;
375

376
    /// FE spaces of this problem.
Praetorius, Simon's avatar
Praetorius, Simon committed
377
    std::shared_ptr<FeSpaces> feSpaces; // eventuell const
378

379
    /// A FileWriter object
Praetorius, Simon's avatar
Praetorius, Simon committed
380
    std::shared_ptr<FileWriter<Traits>> filewriter;
381

382
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
383
    std::shared_ptr<LinearSolverType> linearSolver;
384

385
    /// A block-matrix that is filled during assembling
Praetorius, Simon's avatar
Praetorius, Simon committed
386
    std::shared_ptr<SystemMatrixType> systemMatrix;
387

388
    /// A block-vector with the solution components
Praetorius, Simon's avatar
Praetorius, Simon committed
389
    std::shared_ptr<SystemVectorType> solution;
390
391

    /// A block-vector (load-vector) corresponding to the right.hand side
392
    /// of the equation, filled during assembling
Praetorius, Simon's avatar
Praetorius, Simon committed
393
    std::shared_ptr<SystemVectorType> rhs;
394
395


396
397
  private: // some internal data-structures

398
    template <class T>
399
    using MatrixEntries = Dune::FieldMatrix<T, nComponents, nComponents>;
400

401
    template <class T>
402
    using VectorEntries = Dune::FieldVector<T, nComponents>;
403

404
405
406
407
408

    /// Helper class to store an operator with optionally factors and
    /// a boundary indicator
    template <class OperatorType>
    struct Scaled
409
    {
410
411
412
413
      std::shared_ptr<OperatorType> op;
      double* factor = nullptr;
      double* estFactor = nullptr;
      BoundaryType b = {0};
414
    };
415

416
417
    // lists of operators on elements, boundary and intersection
    struct Operators
418
    {
419
420
421
      std::list<Scaled<ElementOperator>> element;
      std::list<Scaled<IntersectionOperator>> boundary;
      std::list<Scaled<IntersectionOperator>> intersection;
422

423
424
425
426
      std::list<std::shared_ptr<DirichletBC<WorldVector>>> dirichlet; // boundary conditions

      bool assembled = false; // if false, do reassemble
      bool changing = false; // if true, or vectorAssembled false, do reassemble
427
    };
428
429
430

    MatrixEntries<Operators> matrixOperators;
    VectorEntries<Operators> rhsOperators;
431
  };
432
433


434
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
435
436
  extern template class ProblemStat<TestTraits<2,1>>; // 1 component with lagrange degree 1
  extern template class ProblemStat<TestTraits<2,1,2>>; // 2 components with different polynomial degree
437
#endif
438

439
} // end namespace AMDiS
440
441
442

#include "ProblemStat.inc.hpp"