ProblemStat.hpp 11.1 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 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/DirichletBC.hpp>
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
#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>

#include <amdis/utility/TreePath.hpp>
37

38
namespace AMDiS
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
    using GlobalBasis = typename Traits::GlobalBasis;
50 51
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
52

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

    /// Dimension of the mesh
56
    static constexpr int dim = Grid::dimension;
57

58
    /// Dimension of the world
59
    static constexpr int dow = Grid::dimensionworld;
60

61
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
62

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

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

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

78 79
    /// Constructor taking additionally a reference to a mesh that is used
    /// instead of the default created mesh, \ref ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
80
    ProblemStat(std::string name, Grid& grid)
81
      : ProblemStat(name)
82
    {
83
      this->grid = Dune::stackobject_to_shared_ptr(grid);
84
      Parameters::get(name + "->mesh", gridName);
85 86
    }

87

88
    /**
89
     * \brief Initialisation of the problem.
90
     *
91 92 93
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
94
    void initialize(Flag initFlag,
95
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
96
                    Flag adoptFlag = INIT_NOTHING);
97

98 99

    /// Adds an operator to \ref A.
100
    /** @{ */
101 102
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath, ColTreePath,
103
                           double* factor = nullptr, double* estFactor = nullptr);
104

105
    // operator evaluated on the boundary
106 107
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath, ColTreePath,
108
                           double* factor = nullptr, double* estFactor = nullptr);
109
    /** @} */
110 111 112


    /// Adds an operator to \ref rhs.
113
    /** @{ */
114 115
    template <class Operator, class TreePath>
    void addVectorOperator(Operator const& op, TreePath,
116
                           double* factor = nullptr, double* estFactor = nullptr);
117

118
    // operator evaluated on the boundary
119 120
    template <class Operator, class TreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath,
121
                           double* factor = nullptr, double* estFactor = nullptr);
122
    /** @} */
123

124

125
    /// Adds a Dirichlet boundary condition
126
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
127
    void addDirichletBC(Predicate const& predicate,
128
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
129
                        Values const& values);
130

131
  public:
132

133
    /// Implementation of \ref ProblemStatBase::solve
134
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
135 136
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
137 138


139
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
140
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
141 142 143
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
144

145 146

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

Praetorius, Simon's avatar
Praetorius, Simon committed
149

150 151 152 153 154 155 156 157 158 159 160 161 162
  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
    {
163 164 165 166 167
      // for (std::size_t i = 0; i < getNumComponents(); ++i)
      //   if (adaptInfo.spaceToleranceReached(i))
      //     adaptInfo.allowRefinement(false, i);
      //   else
      //     adaptInfo.allowRefinement(true, i);
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190

      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; }


191
  public: // get-methods
192

193
    /// Returns a pointer to system-matrix, \ref systemMatrix
194 195
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
196 197


198
    std::shared_ptr<SystemVector> getSolutionVector() const
199 200 201 202 203 204 205 206 207 208 209
    {
      return solution;
    }

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

211 212 213
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
214
    {
215
      auto&& tp = makeTreePath(path);
216
      return makeDOFVectorView(*solution, tp);
217
    }
218 219


220
    /// Return a point to the rhs system-vector, \ref rhs
221 222
    auto getRhsVector()       { return rhs; }
    auto getRhsVector() const { return rhs; }
223 224


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

Praetorius, Simon's avatar
Praetorius, Simon committed
228
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
229
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
230
      linearSolver = solver;
231 232
    }

233 234
    /// Return a pointer to the grid, \ref grid
    std::shared_ptr<Grid> getGrid() { return grid; }
235

236 237
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
238
    void setGrid(Grid& grid_)
239
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
240
      grid = Dune::stackobject_to_shared_ptr(grid_);
241
      createGlobalBasis();
242 243 244 245 246
      createMatricesAndVectors();
      createFileWriter();
    }


247
    /// Return the gridView of the leaf-level
248
    auto const& gridView() { return globalBasis->gridView(); }
249

250
    /// Return the \ref feSpaces
251
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis; }
252 253


254
    /// Implementation of \ref ProblemStatBase::getName
255 256 257 258
    virtual std::string getName() const override
    {
      return name;
    }
259

260
  protected: // initialization methods
261

262
    void createGrid()
263
    {
264 265 266
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
      test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
267

268
      grid = MeshCreator<Grid>::create(gridName);
269

270 271 272 273
      msg("Create grid:");
      msg("#elements = "   , grid->size(0));
      msg("#faces/edges = ", grid->size(1));
      msg("#vertices = "   , grid->size(dim));
274
      msg("");
275
    }
276

277
    void createGlobalBasis()
278
    {
279
      globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid->leafGridView()));
280 281 282 283 284 285
      initGlobalBasis(*globalBasis);
    }

    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
      globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis.localView().tree());
286 287 288
      matrixOperators.init(*globalTree);
      rhsOperators.init(*globalTree);
      constraints.init(*globalTree);
289
    }
290

291 292
    void createMatricesAndVectors()
    {
293
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
294 295
      solution = std::make_shared<SystemVector>(*globalBasis, "solution");
      rhs = std::make_shared<SystemVector>(*globalBasis, "rhs");
296
    }
297

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

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

306
      linearSolver = solverCreator->create(name + "->solver");
307
    }
308

309
    void createFileWriter();
310

Praetorius, Simon's avatar
Praetorius, Simon committed
311

312
  private:
313 314 315
    /// Name of this problem.
    std::string name;

316 317
    /// Grid of this problem.
    std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
318

319
    /// Name of the mesh
320
    std::string gridName = "none";
321

322
    /// FE spaces of this problem.
323
    std::shared_ptr<GlobalBasis> globalBasis;
324
    std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
325

326
    /// A FileWriter object
327
    std::list<std::shared_ptr<FileWriterInterface>> filewriter;
328

329
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
330
    std::shared_ptr<LinearSolverType> linearSolver;
331

332
    /// A block-matrix that is filled during assembling
333
    std::shared_ptr<SystemMatrix> systemMatrix;
334

335
    /// A block-vector with the solution components
336
    std::shared_ptr<SystemVector> solution;
337 338

    /// A block-vector (load-vector) corresponding to the right.hand side
339
    /// of the equation, filled during assembling
340
    std::shared_ptr<SystemVector> rhs;
341 342


343 344
  private: // some internal data-structures

345 346
    MatrixOperators<GlobalBasis> matrixOperators;
    VectorOperators<GlobalBasis> rhsOperators;
347

348
    Constraints<GlobalBasis> constraints;
349
  };
350 351


352
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
353
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
354
#endif
355

356
} // end namespace AMDiS
357 358 359

#include "ProblemStat.inc.hpp"