ProblemStat.hpp 14.5 KB
Newer Older
1 2 3 4
#pragma once

#include <list>
#include <map>
5 6 7 8 9
#include <memory>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
10

11 12
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
13 14
#include <dune/grid/common/grid.hh>

15 16 17 18
#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/DirichletBC.hpp>
19
//#include <amdis/Estimator.hpp>
20 21 22
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
23
#include <amdis/Marker.cpp>
24 25 26 27 28 29 30 31 32 33 34 35 36 37
#include <amdis/Mesh.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>

#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/TypeDefs.hpp>
#include <amdis/common/Utility.hpp>

#include <amdis/GridFunctions.hpp>
#include <amdis/gridfunctions/DOFVectorView.hpp>

#include <amdis/io/FileWriterInterface.hpp>

38
#include <amdis/utility/TreeData.hpp>
39
#include <amdis/utility/TreePath.hpp>
40

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

50
  public: // typedefs and static constants
51

52
    using GlobalBasis = typename Traits::GlobalBasis;
53 54 55
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
    using Element  = typename GridView::template Codim<0>::Entity;
56
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
57

58
    /// Dimension of the grid
59
    static constexpr int dim = Grid::dimension;
60

61
    /// Dimension of the world
62
    static constexpr int dow = Grid::dimensionworld;
63

64
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
65
    using SystemVector = DOFVector<Traits, double>;
66

67
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
68

69
  public:
70 71
    /**
     * \brief Constructor. Takes the name of the problem that is used to
72
     * access values corresponding to this problem in the parameter file.
73
     **/
74
    explicit ProblemStat(std::string name)
75
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
76
      , name_(std::move(name))
77
    {}
78

79 80
    /// Constructor taking additionally a reference to a grid that is used
    /// instead of the default created grid, \ref ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
81
    ProblemStat(std::string name, Grid& grid)
82
      : ProblemStat(std::move(name))
83
    {
84 85 86 87 88 89 90 91 92 93 94
      grid_ = Dune::stackobject_to_shared_ptr(grid);
      Parameters::get(name_ + "->mesh", gridName_);
    }

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
    ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
      : ProblemStat(std::move(name), grid)
    {
      globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
      initGlobalBasis(*globalBasis_);
95 96
    }

97

98
    /**
99
     * \brief Initialisation of the problem.
100
     *
101 102 103
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
104
    void initialize(Flag initFlag,
105
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
106
                    Flag adoptFlag = INIT_NOTHING);
107

108 109

    /// Adds an operator to \ref A.
110
    /** @{ */
111 112
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath, ColTreePath,
113
                           double* factor = nullptr, double* estFactor = nullptr);
114

115
    // operator evaluated on the boundary
116 117
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath, ColTreePath,
118
                           double* factor = nullptr, double* estFactor = nullptr);
119
    /** @} */
120 121 122


    /// Adds an operator to \ref rhs.
123
    /** @{ */
124 125
    template <class Operator, class TreePath>
    void addVectorOperator(Operator const& op, TreePath,
126
                           double* factor = nullptr, double* estFactor = nullptr);
127

128
    // operator evaluated on the boundary
129 130
    template <class Operator, class TreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath,
131
                           double* factor = nullptr, double* estFactor = nullptr);
132
    /** @} */
133

134

135
    /// Adds a Dirichlet boundary condition
136
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
137
    void addDirichletBC(Predicate const& predicate,
138
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
139
                        Values const& values);
140

141
  public:
142

143
    /// Implementation of \ref ProblemStatBase::solve
144
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
145 146
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
147 148


149 150 151 152 153 154 155 156 157 158 159 160
    /// Implementation of \ref ProblemStatBase::estimate.
    virtual void estimate(AdaptInfo& adaptInfo) override;


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


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


161
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
162
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
163 164 165
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
166

167 168

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

Praetorius, Simon's avatar
Praetorius, Simon committed
171

172
  public: // get-methods
173

174 175 176
    /// Returns a pointer to system-matrix, \ref systemMatrix_
    std::shared_ptr<SystemMatrix> getSystemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
177

178
    /// Returns a pointer to the solution vector, \ref solution_
179
    std::shared_ptr<SystemVector> getSolutionVector() const
180
    {
181
      return solution_;
182 183 184 185 186 187 188
    }

    /// Return a mutable view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {})
    {
      auto&& tp = makeTreePath(path);
189
      return makeDOFVectorView(*solution_, tp);
190
    }
191

192 193 194
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
195
    {
196
      auto&& tp = makeTreePath(path);
197
      return makeDOFVectorView(*solution_, tp);
198
    }
199 200


201
    /// Return a point to the rhs system-vector, \ref rhs
202 203
    std::shared_ptr<SystemVector> getRhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
204 205


206
    /// Return a pointer to the linear solver, \ref linearSolver
207
    std::shared_ptr<LinearSolverType> getSolver() { return linearSolver_; }
208

Praetorius, Simon's avatar
Praetorius, Simon committed
209
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
210
    {
211
      linearSolver_ = solver;
212 213
    }

214
    /// Return a pointer to the grid, \ref grid
215
    std::shared_ptr<Grid> getGrid() { return grid_; }
216

217
    /// Set the grid. Stores pointer to passed reference and initializes feSpaces
218
    /// matrices and vectors, as well as the file-writer.
219
    void setGrid(Grid& grid)
220
    {
221 222
      grid_ = Dune::stackobject_to_shared_ptr(grid);

223
      createGlobalBasis();
224 225 226 227
      createMatricesAndVectors();
      createFileWriter();
    }

228
/*
229
    /// Return the gridView of the leaf-level
230 231 232 233
    auto leafGridView() { return grid_->leafGridView(); }

    /// Return the gridView of levle `level`
    auto levelGridView(int level) { return grid_->levelGridView(level); }
234
*/
235
    /// Return the \ref feSpaces
236
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
237 238


239
    /// Implementation of \ref ProblemStatBase::getName
240 241
    virtual std::string getName() const override
    {
242
      return name_;
243
    }
244

245
  protected: // initialization methods
246

247
    void createGrid()
248
    {
249 250
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
251
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
252

253
      grid_ = MeshCreator<Grid>::create(gridName_);
254

255
      msg("Create grid:");
256 257 258
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
259
      msg("");
260
    }
261

262
    void createGlobalBasis()
263
    {
264 265
      globalBasis_ = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid_->leafGridView()));
      initGlobalBasis(*globalBasis_);
266 267 268 269
    }

    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
270 271 272 273
      localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
      matrixOperators_.init(localView_->tree(), tag::store{});
      rhsOperators_.init(localView_->tree(), tag::store{});
      constraints_.init(localView_->tree(), tag::store{});
274
    }
275

276 277
    void createMatricesAndVectors()
    {
278 279 280
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
281

282 283 284 285 286 287 288 289 290
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
      solution = std::make_shared<SystemVector>(*globalBasis, "solution");

      auto localView = globalBasis->localView();
      forEachNode(localView.tree(), [&,this](auto const& node, auto treePath)
      {
        std::string i = to_string(treePath);
        estimates[i].resize(globalBasis->gridView().indexSet().size(0));
        for (std::size_t j = 0; j < estimates[i].size(); j++)
291
          estimates[i][j] = 0.0; // TODO: Remove when estimate() is implemented
292 293 294
      });

      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
295
    }
296

297
    void createSolver()
298
    {
299
      std::string solverName = "cg";
300
      Parameters::get(name_ + "->solver->name", solverName);
301

302
      auto solverCreator
303
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
304

305
      linearSolver_ = solverCreator->create(name_ + "->solver");
306
    }
307

308
    void createEstimator()
309
    {/* COPIED FROM AMDiS
310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331
      for (std::size_t i = 0, i < nComponents, i++) {
        std::string estName = "";
        Parameters::get(name + "->estimator->name[" + std::to_string(i) + "]", estName);
        auto estimatorCreator
          = named(CreatorMap<Estimator>::getCreator(estName, name + "->estimator->name[" + std::to_string(i) + "]"));
        if (estimatorCreator) {
	        estimatorCreator->setName(estName);
	        estimatorCreator->setRow(i);
	        estimatorCreator->setSolution(solution->getDOFVector(i));
	        estimator[i] = estimatorCreator->create();
        }

        if (estimator[i]) {
          for (std::size_t j = 0; j < nComponents; j++)
	          estimator[i]->addSystem((*systemMatrix)[i][j],
              solution->getDOFVector(j),
				      rhs->getDOFVector(i));
        }
      }*/
    }


332
    void createMarker();
333

334
    void createFileWriter();
335

Praetorius, Simon's avatar
Praetorius, Simon committed
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
  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
    {
      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; }


372
  private:
373
    /// Name of this problem.
374
    std::string name_;
375

376
    /// Grid of this problem.
377
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
378

379 380 381
    /// Number of grids
    int nGrids = 1;

382
    /// Name of the grid
383
    std::string gridName_ = "none";
384

385
    /// FE spaces of this problem.
386 387
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
388

389
    /// A FileWriter object
390
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
391

392
    /// Pointer to the adaptation markers
393
    std::list<std::shared_ptr<Marker<Traits> > > marker;
394 395 396 397

    /// Pointer to the estimators for this problem
//    std::vector<Estimator*> estimator;

398
    /// An object of the linearSolver Interface
399
    std::shared_ptr<LinearSolverType> linearSolver_;
400

401
    /// A block-matrix that is filled during assembling
402
    std::shared_ptr<SystemMatrix> systemMatrix_;
403

404
    /// A block-vector with the solution components
405
    std::shared_ptr<SystemVector> solution_;
406

407
    /// A vector with the local element error estimates
408
    /// reverse indexed by [component index][element index]
409
    std::map<std::string, std::vector<double> > estimates;
410

411
    /// A block-vector (load-vector) corresponding to the right.hand side
412
    /// of the equation, filled during assembling
413
    std::shared_ptr<SystemVector> rhs_;
414 415


416 417
  private: // some internal data-structures

418 419 420
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
421 422

//    VectorData<GlobalBasis, TODO: NodeData> estimates_;
423
  };
424

425 426 427 428 429 430 431 432 433 434 435 436 437 438 439
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
  // Deduction rule
  template <class Grid, class GlobalBasis>
  ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
    -> ProblemStat<DefaultProblemTraits<GlobalBasis>>;
#endif

  // Generator for ProblemStat with given Grid and GlobalBasis
  template <class Grid, class GlobalBasis>
  ProblemStat<DefaultProblemTraits<GlobalBasis>>
  makeProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
  {
    return {std::move(name), grid, globalBasis};
  }

440

441
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
442
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
443
#endif
444

445
} // end namespace AMDiS
446 447 448

#include "ProblemStat.inc.hpp"