ProblemStat.hpp 14.9 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
    using MeshView = typename Mesh::LeafGridView;
47

48
49
50
    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
    template <size_t I>
    using FeSpace = std::tuple_element_t<I, FeSpaces>;
67

68
    using WorldVector = typename Codim0::Geometry::GlobalCoordinate;
69

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
    /**
     * \brief Constructor. Takes the name of the problem that is used to
80
     * access values correpsonding to this püroblem in the parameter file.
81
     *
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
     * \brief Initialisation of the problem.
96
     *
97
98
99
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
100
101
    void initialize(Flag initFlag,
		    Self* adoptProblem = NULL,
102
		    Flag adoptFlag = INIT_NOTHING);
103

104
105

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

112
    void addMatrixOperator(shared_ptr<OperatorType> op,
113
                           int i, int j,
114
                           double* factor = NULL,
115
                           double* estFactor = NULL);
116
117

    void addMatrixOperator(std::map< std::pair<int,int>,
118
                                     shared_ptr<OperatorType> > ops);
119
    void addMatrixOperator(std::map< std::pair<int,int>,
120
121
                                     std::pair<shared_ptr<OperatorType>, bool> > ops);
    /** @} */
122
123
124


    /// Adds an operator to \ref rhs.
125
    /** @{ */
126
    void addVectorOperator(OperatorType& op,
127
                           int i,
128
                           double* factor = NULL,
129
                           double* estFactor = NULL);
130
131

    void addVectorOperator(shared_ptr<OperatorType> op,
132
                           int i,
133
                           double* factor = NULL,
134
                           double* estFactor = NULL);
135
136

    void addVectorOperator(std::map< int,
137
                                     shared_ptr<OperatorType> > ops);
138
    void addVectorOperator(std::map< int,
139
140
                                     std::pair<shared_ptr<OperatorType>, bool> > ops);
    /** @} */
141
142
143

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

148

149
    /// Implementation of \ref ProblemStatBase::solve
150
151
    virtual void solve(AdaptInfo& adaptInfo,
		       bool createMatrixData = true,
152
		       bool storeMatrixData = false) override;
153
154


155
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
156
157
158
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
				   Flag flag,
				   bool asmMatrix = true,
159
				   bool asmVector = true) override;
160

161
162

    /// 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
    {
      return nComponents;
    }
172
173


174
175
176
    /// Return the \ref systemMatrix
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
177
178


179
180
181
    /// Return the \ref solution
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }
182

183
    /// 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
    auto getMeshView() { return meshView; }
204

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
    virtual std::string getName() const override
    {
      return name;
    }
221

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
      meshName = "";
      Parameters::get(name + "->mesh", meshName);
234
235
236
      test_exit(!meshName.empty(),
                "No mesh name specified for '", name, "->mesh'!");

237
      mesh = MeshCreator<Mesh>::create(meshName);
238
      meshView = make_shared<MeshView>(mesh->leafGridView());
239
240
241
242
243
244

      msg("Create mesh:");
      msg("#elements = "   , mesh->size(0));
      msg("#faces/edges = ", mesh->size(1));
      msg("#vertices = "   , mesh->size(dim));
      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
275
      filewriter = make_shared<FileWriter<Traits>>(name + "->output",
                                                        *meshView,
276
                                                        componentNames);
277
    }
278

279
  protected: // sub-methods to assemble DOFMatrix
280

281
282
    template <class LhsData, class RhsData, class Elements>
    void assemble(LhsData lhs, RhsData rhs, Elements const& elements);
283

284
    template <class RowView, class ColView>
285
    bool getElementMatrix(RowView const& rowView,
286
287
			  ColView const& colView,
			  ElementMatrix& elementMatrix,
288
			  std::list<shared_ptr<OperatorType>>& operators,
289
                          std::list<double*> const& factors);
290

291
    template <class RowView>
292
    bool getElementVector(RowView const& rowView,
293
			  ElementVector& elementVector,
294
			  std::list<shared_ptr<OperatorType>>& operators,
295
                          std::list<double*> const& factors);
296
297


298
299
300
301
302
    template <class Matrix, class RowIndexSet, class ColIndexSet>
    void addElementMatrix(Matrix& matrix,
			  RowIndexSet const& rowIndexSet,
			  ColIndexSet const& colIndexSet,
			  ElementMatrix const& elementMatrix);
303

304
305
306
307
    template <class Vector, class IndexSet>
    void addElementVector(Vector& vector,
			  IndexSet const& indexSet,
			  ElementVector const& elementvector);
308

309
  private: // some internal methods
310

311
    template <size_t I = 0>
312
    typename FeSpace<I>::NodeFactory&
313
314
    getNodeFactory(const index_<I> = index_<I>())
    {
315
      return const_cast<typename FeSpace<I>::NodeFactory&>(std::get<I>(*feSpaces).nodeFactory());
316
    }
317
318


319
  private:
320
321
322
323
324
325
    /// 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;
326

327
    /// Mesh of this problem.
328
    shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
329

330
    /// Name of the mesh
331
    std::string meshName;
332

333
    /// A gridView object
334
    shared_ptr<MeshView> meshView;
335

336
337
    /// Pointer to the meshes for the different problem components
    std::vector<Mesh*> componentMeshes;
338

339
    /// FE spaces of this problem.
340
    shared_ptr<FeSpaces> feSpaces; // eventuell const
341

342
    /// A FileWriter object
343
    shared_ptr<FileWriter<Traits>> filewriter;
344

345
    /// An object of the linearSolver Interface
346
    shared_ptr<LinearSolverType> linearSolver;
347

348
    /// A block-matrix that is filled during assembling
349
    shared_ptr<SystemMatrixType> systemMatrix;
350

351
    /// A block-vector with the solution components
352
    shared_ptr<SystemVectorType> solution;
353
354

    /// A block-vector (load-vector) corresponding to the right.hand side
355
    /// of the equation, filled during assembling
356
    shared_ptr<SystemVectorType> rhs;
357
358


359
360
    template <class T>
    using MatrixEntries = std::map< std::pair<int,int>, std::list<T> >;
361

362
363
    template <class T>
    using VectorEntries = std::map< int, std::list<T> >;
364
365
366


    /// A map (i,j) -> list<DirichleBC> string a boundary conditions for
367
    /// each matrix block
368
    MatrixEntries<shared_ptr<DirichletBC<WorldVector>>> dirichletBc;
369
370

    /// A map (i,j) -> list<OperatorType> string the differential operators for
371
    /// each matrix block
372
373
    MatrixEntries<shared_ptr<OperatorType>> matrixOperators;
    MatrixEntries<double*>                  matrixFactors;
374
375
    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
376
377

    /// A map (i) -> list<OperatorType> string the differential operators for
378
    /// each rhs block
379
380
    VectorEntries<shared_ptr<OperatorType>> vectorOperators;
    VectorEntries<double*>                  vectorFactors;
381
382
    std::map< int, bool >                   vectorAssembled; // if false, do reassemble
    std::map< int, bool >                   vectorChanging;  // if true, or vectorAssembled false, do reassemble
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
    
    template <int r, int c>
    struct MatrixData 
    {
      using DOFMatrixType = 
        std::tuple_element_t<c, std::tuple_element_t<r, typename SystemMatrixType::DOFMatrices>>;
        
      DOFMatrixType&                       matrix;
      std::list<shared_ptr<OperatorType>>& operators;
      std::list<double*> const&            factors;
      bool                                 assemble;
      
      std::pair<int,int> row_col = {r, c};
    };
    
    template <int r>
    struct VectorData
    {
      using DOFVectorType = 
        std::tuple_element_t<r, typename SystemVectorType::DOFVectors>;
      
      DOFVectorType&                       vector;
      std::list<shared_ptr<OperatorType>>& operators;
      std::list<double*> const&            factors;
      bool                                 assemble;
      
      int row = r;
    };
411
  };
412
413


414
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
415
416
417
418
//   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>>;
419
#endif
420

421
  namespace Impl
422
423
424
425
426
427
428
429
430
  {
    template <class ProblemStatType>
    struct ProblemStat : public ProblemStatType,
                         public StandardProblemIteration
    {
      using ProblemStatType::getName;

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

435
      /**
436
437
438
439
440
441
       * \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.
       *
442
443
       * Implementation of \ref StandardProblemIteration::oneIteration.
       **/
444
445
      virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
      {
446
        for (size_t i = 0; i < ProblemStatType::getNumComponents(); ++i)
447
448
449
450
451
452
453
          if (adaptInfo.spaceToleranceReached(i))
            adaptInfo.allowRefinement(false, i);
          else
            adaptInfo.allowRefinement(true, i);

        return StandardProblemIteration::oneIteration(adaptInfo, toDo);
      }
454

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

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

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

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

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

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

474
  } // end namespace Impl
475

476
  template <class Traits>
477
  using ProblemStat = Impl::ProblemStat<ProblemStatSeq<Traits>>;
478

479
} // end namespace AMDiS
480
481
482

#include "ProblemStat.inc.hpp"