ProblemStat.hpp 15.4 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
#include <dune/amdis/StandardProblemIteration.hpp>

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

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

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

46
  public: // typedefs and static constants
47

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

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

55
56
    using ElementVector = Impl::ElementVector;
    using ElementMatrix = Impl::ElementMatrix;
57
58


59
    /// Dimension of the mesh
60
    static constexpr int dim = Traits::dim;
61

62
    /// Dimension of the world
63
    static constexpr int dow = Traits::dow;
64

65
    /// Number of problem components
66
    static constexpr int nComponents = Traits::nComponents;
67
68


69
    template <std::size_t I>
70
    using FeSpace = std::tuple_element_t<I, FeSpaces>;
71

72
    using WorldVector = typename Codim0::Geometry::GlobalCoordinate;
73

74
75
    using SystemVectorType = SystemVector<FeSpaces>;
    using SystemMatrixType = SystemMatrix<FeSpaces>;
76

77
78
    using ElementOperator = Operator<MeshView, Element>;
    using IntersectionOperator = Operator<MeshView, typename MeshView::Intersection>;
79

80
81
    using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
                                                   typename SystemVectorType::MultiVector>;
82

83
84
    /// Helper class to store an operator with optionally factors and
    /// a boundary indicator
85
86
    template <class OperatorType>
    struct Scaled
87
88
    {
      std::shared_ptr<OperatorType> op;
89
90
91
      double* factor = nullptr;
      double* estFactor = nullptr;
      BoundaryType b = {0};
92
93
    };

94
  public:
95
96
    /**
     * \brief Constructor. Takes the name of the problem that is used to
97
     * access values correpsonding to this püroblem in the parameter file.
98
     *
99
     * Parameters read by ProblemStat, with name 'PROB'
100
     *   PROB->names: \ref componentNames
101
     **/
102
103
104
    ProblemStat(std::string name)
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      , name(name)
105
    {
106
      Parameters::get(name + "->names", componentNames);
107
      for (std::size_t i = componentNames.size(); i < nComponents; ++i)
108
        componentNames.push_back("solution[" + std::to_string(i) + "]");
109
    }
110

111
112
113
114
    /// 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)
115
116
    {
      this->mesh = std::shared_ptr<Mesh>(&mesh, optional_delete(false));
117
118
119
120
      componentMeshes.resize(nComponents, this->mesh.get());

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

123

124
    /**
125
     * \brief Initialisation of the problem.
126
     *
127
128
129
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
130
    void initialize(Flag initFlag,
131
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
132
                    Flag adoptFlag = INIT_NOTHING);
133

134
135

    /// Adds an operator to \ref A.
136
    /** @{ */
137
    template <class OperatorType>
138
139
    void addMatrixOperator(OperatorType const& op, int i, int j,
                           double* factor = nullptr, double* estFactor = nullptr);
140

141
    // operator evaluated on the boundary
142
    void addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j,
143
                           double* factor = nullptr, double* estFactor = nullptr);
144

Praetorius, Simon's avatar
Praetorius, Simon committed
145
146
    void addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops);
    void addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator, bool> > ops);
147
    /** @} */
148
149
150


    /// Adds an operator to \ref rhs.
151
    /** @{ */
152
    template <class OperatorType>
153
154
    void addVectorOperator(OperatorType const& op, int i,
                           double* factor = nullptr, double* estFactor = nullptr);
155

156
    // operator evaluated on the boundary
157
    void addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i,
158
                           double* factor = nullptr, double* estFactor = nullptr);
159

Praetorius, Simon's avatar
Praetorius, Simon committed
160
161
    void addVectorOperator(std::map< int, ElementOperator> ops);
    void addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops);
162
    /** @} */
163
164
165

    /// Adds a Dirichlet boundary condition
    template <class Predicate, class Values>
166
    void addDirichletBC(Predicate const& predicate,
Praetorius, Simon's avatar
Praetorius, Simon committed
167
168
                        int row, int col,
                        Values const& values);
169

170

171
    /// Implementation of \ref ProblemStatBase::solve
172
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
173
174
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
175
176


177
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
178
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
179
180
181
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
182

183
184

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

Praetorius, Simon's avatar
Praetorius, Simon committed
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
221
222
223
224
225
226
227
228
  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; }


229
  public: // get-methods
230
231

    /// Returns nr of components \ref nComponents
232
    static constexpr int getNumComponents()
233
234
235
    {
      return nComponents;
    }
236
237


238
    /// Returns a pointer to system-matrix, \ref systemMatrix
239
240
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
241
242


243
    /// Return a pointer to the solution, \ref solution
244
245
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }
246

247
    /// Return the i'th \ref solution component
248
249
    template <std::size_t I = 0>
    decltype(auto) getSolution(const index_t<I> _i = {})
250
    {
251
      return (*solution)[_i];
252
    }
253
254


255
    /// Return a point to the rhs system-vector, \ref rhs
256
257
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
258
259


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

Praetorius, Simon's avatar
Praetorius, Simon committed
263
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
264
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
265
      linearSolver = solver;
266
267
    }

268
269
    /// Return a pointer to the mesh, \ref mesh
    std::shared_ptr<Mesh> getMesh() { return mesh; }
270

271
272
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
273
    void setMesh(Mesh& mesh_)
274
    {
275
      mesh = std::shared_ptr<Mesh>(&mesh_, optional_delete(false));
276
      std::fill(componentMeshes.begin(), componentMeshes.end(), this->mesh.get());
277
278
279
280
281
282
283

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


284
285
286
287
288
    /// 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); }
289

290
291
    /// Return the \ref feSpaces
    auto getFeSpaces() { return feSpaces; }
292
293


294
    /// Return the I'th \ref feSpaces component
295
296
297
298
299
300
301
302
    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> = {})
303
    {
304
      return std::get<I>(*feSpaces);
305
    }
306

307
    /// Implementation of \ref ProblemStatBase::getName
308
309
310
311
    virtual std::string getName() const override
    {
      return name;
    }
312

313
    /// Return a vector of names of the problem components
314
315
316
317
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
318

319
  protected: // initialization methods
320

321
322
    void createMesh()
    {
323
324
      meshName = "";
      Parameters::get(name + "->mesh", meshName);
325
      test_exit(!meshName.empty(), "No mesh name specified for '", name, "->mesh'!");
326

327
      mesh = MeshCreator<Mesh>::create(meshName);
328
329
330
331
332
333

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

336
337
    void createFeSpaces()
    {
338
      feSpaces = std::make_shared<FeSpaces>(constructTuple<FeSpaces>(leafGridView()));
339
    }
340

341
342
    void createMatricesAndVectors()
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
343
344
      systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
      solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);
345

346
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
Praetorius, Simon's avatar
Praetorius, Simon committed
347
      rhs = std::make_shared<SystemVectorType>(*feSpaces, rhsNames);
348
    }
349

350
    void createSolver()
351
    {
352
353
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
354

355
356
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
357

358
      linearSolver = solverCreator->create(name + "->solver");
359
    }
360

361
362
    void createFileWriter()
    {
363
      filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", leafGridView(), componentNames);
364
    }
365

Praetorius, Simon's avatar
Praetorius, Simon committed
366

367
  protected: // sub-methods to assemble DOFMatrix
368

369
370
    template <class LhsData, class RhsData, class GV>
    void assemble(LhsData lhs, RhsData rhs, GV const& gridView);
371

372
    template <class RowView, class ColView>
373
    bool getElementMatrix(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
374
375
                          ColView const& colView,
                          ElementMatrix& elementMatrix,
376
377
378
                          std::list<Scaled<ElementOperator>>& operators,
                          std::list<Scaled<IntersectionOperator>>& boundary_operators,
                          std::list<Scaled<IntersectionOperator>>& intersection_operators);
379

380
    template <class RowView>
381
    bool getElementVector(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
382
                          ElementVector& elementVector,
383
384
385
                          std::list<Scaled<ElementOperator>>& operators,
                          std::list<Scaled<IntersectionOperator>>& boundary_operators,
                          std::list<Scaled<IntersectionOperator>>& intersection_operators);
386
387


388
389
    template <class Matrix, class RowIndexSet, class ColIndexSet>
    void addElementMatrix(Matrix& matrix,
Praetorius, Simon's avatar
Praetorius, Simon committed
390
391
392
                          RowIndexSet const& rowIndexSet,
                          ColIndexSet const& colIndexSet,
                          ElementMatrix const& elementMatrix);
393

394
395
    template <class Vector, class IndexSet>
    void addElementVector(Vector& vector,
Praetorius, Simon's avatar
Praetorius, Simon committed
396
397
                          IndexSet const& indexSet,
                          ElementVector const& elementvector);
398
399


400
  private:
401
402
403
404
405
406
    /// 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;
407

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

411
    /// Name of the mesh
412
    std::string meshName = "none";
413

414
415
    /// Pointer to the meshes for the different problem components
    std::vector<Mesh*> componentMeshes;
416

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

420
    /// A FileWriter object
Praetorius, Simon's avatar
Praetorius, Simon committed
421
    std::shared_ptr<FileWriter<Traits>> filewriter;
422

423
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
424
    std::shared_ptr<LinearSolverType> linearSolver;
425

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

429
    /// A block-vector with the solution components
Praetorius, Simon's avatar
Praetorius, Simon committed
430
    std::shared_ptr<SystemVectorType> solution;
431
432

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


437
438
  private: // some internal data-structures

439
    template <class T>
440
    using MatrixEntries = Dune::FieldMatrix<T, nComponents, nComponents>;
441

442
    template <class T>
443
    using VectorEntries = Dune::FieldVector<T, nComponents>;
444

445
446
447
448
449

    /// Helper class to store an operator with optionally factors and
    /// a boundary indicator
    template <class OperatorType>
    struct Scaled
450
    {
451
452
453
454
      std::shared_ptr<OperatorType> op;
      double* factor = nullptr;
      double* estFactor = nullptr;
      BoundaryType b = {0};
455
    };
456

457
458
    // lists of operators on elements, boundary and intersection
    struct Operators
459
    {
460
461
462
      std::list<Scaled<ElementOperator>> element;
      std::list<Scaled<IntersectionOperator>> boundary;
      std::list<Scaled<IntersectionOperator>> intersection;
463

464
465
466
467
      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
468
    };
469
470
471

    MatrixEntries<Operators> matrixOperators;
    VectorEntries<Operators> rhsOperators;
472
  };
473
474


475
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
476
477
  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
478
#endif
479

480
} // end namespace AMDiS
481
482
483

#include "ProblemStat.inc.hpp"