ProblemStat.hpp 12.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
    static const std::size_t nComponents = 1; // TODO: count leaf nodes in GlobalBasis
56

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

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

63
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
64

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

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

70
  public:
71 72
    /**
     * \brief Constructor. Takes the name of the problem that is used to
73
     * access values correpsonding to this püroblem in the parameter file.
74
     *
75
     * Parameters read by ProblemStat, with name 'PROB'
76
     *   PROB->names: \ref componentNames
77
     **/
78
    explicit ProblemStat(std::string name)
79 80
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      , name(name)
81
    {
82
      Parameters::get(name + "->names", componentNames);
83
      for (std::size_t i = componentNames.size(); i < nComponents; ++i)
84
        componentNames.push_back("solution[" + std::to_string(i) + "]");
85
    }
86

87 88
    /// 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
89
    ProblemStat(std::string name, Grid& grid)
90
      : ProblemStat(name)
91
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
92
      this->grid = Dune::stackobject_to_shared_ptr(grid);
93
      componentGrids.resize(nComponents, this->grid.get());
94

95 96
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
97 98
    }

99

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

110 111

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

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


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

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

136

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

143
  public:
144

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


151
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
152
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
153 154 155
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
156

157 158

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

Praetorius, Simon's avatar
Praetorius, Simon committed
161

162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
  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; }


203
  public: // get-methods
204

205
    /// Returns a pointer to system-matrix, \ref systemMatrix
206 207
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
208 209


210
    /// Return a pointer to the solution, \ref solution
211 212
    std::shared_ptr<SystemVector>       getSolution()       { return solution; }
    std::shared_ptr<SystemVector const> getSolution() const { return solution; }
213

214
    /// Return a view to a solution component (as shared_ptr)
215
    template <class TreePath>
216
    auto getSolution(TreePath path) const
217
    {
218
      auto tp = makeTreePath(path);
219
      return makeDOFVectorView(*solution, tp);
220
    }
221 222


223
    /// Return a point to the rhs system-vector, \ref rhs
224 225
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
226 227


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

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

236 237
    /// Return a pointer to the grid, \ref grid
    std::shared_ptr<Grid> getGrid() { return grid; }
238

239 240
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
241
    void setGrid(Grid& grid_)
242
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
243
      grid = Dune::stackobject_to_shared_ptr(grid_);
244
      std::fill(componentGrids.begin(), componentGrids.end(), this->grid.get());
245

246
      createGlobalBasis();
247 248 249 250 251
      createMatricesAndVectors();
      createFileWriter();
    }


252
    /// Return the gridView of the leaf-level
253
    auto leafGridView() { return grid->leafGridView(); }
254 255

    /// Return the gridView of levle `level`
256
    auto levelGridView(int level) { return grid->levelGridView(level); }
257

258
    /// Return the \ref feSpaces
259
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis; }
260 261


262
    /// Implementation of \ref ProblemStatBase::getName
263 264 265 266
    virtual std::string getName() const override
    {
      return name;
    }
267

268
    /// Return a vector of names of the problem components
269 270 271 272
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
273

274
    std::size_t getNumComponents() const
275 276 277 278
    {
      return nComponents;
    }

279
  protected: // initialization methods
280

281
    void createGrid()
282
    {
283 284 285
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
      test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
286

287
      grid = MeshCreator<Grid>::create(gridName);
288

289 290 291 292
      msg("Create grid:");
      msg("#elements = "   , grid->size(0));
      msg("#faces/edges = ", grid->size(1));
      msg("#vertices = "   , grid->size(dim));
293
      msg("");
294
    }
295

296
    void createGlobalBasis()
297
    {
298
      globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(leafGridView()));
299 300 301 302
      globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis->localView().tree());
      matrixOperators.init(*globalTree);
      rhsOperators.init(*globalTree);
      constraints.init(*globalTree);
303
    }
304

305 306
    void createMatricesAndVectors()
    {
307 308
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
      solution = std::make_shared<SystemVector>(*globalBasis, componentNames[0]);
309

310
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
311
      rhs = std::make_shared<SystemVector>(*globalBasis, rhsNames[0]);
312
    }
313

314
    void createSolver()
315
    {
316 317
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
318

319 320
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
321

322
      linearSolver = solverCreator->create(name + "->solver");
323
    }
324

325
    void createFileWriter();
326

Praetorius, Simon's avatar
Praetorius, Simon committed
327

328
  private:
329 330 331 332 333 334
    /// 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;
335

336 337
    /// Grid of this problem.
    std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
338

339
    /// Name of the mesh
340
    std::string gridName = "none";
341

342
    /// Pointer to the meshes for the different problem components
343
    std::vector<Grid*> componentGrids;
344

345
    /// FE spaces of this problem.
346
    std::shared_ptr<GlobalBasis> globalBasis; // eventuell const
347
    std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
348

349
    /// A FileWriter object
350
    std::list<std::shared_ptr<FileWriterInterface>> filewriter;
351

352
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
353
    std::shared_ptr<LinearSolverType> linearSolver;
354

355
    /// A block-matrix that is filled during assembling
356
    std::shared_ptr<SystemMatrix> systemMatrix;
357

358
    /// A block-vector with the solution components
359
    std::shared_ptr<SystemVector> solution;
360 361

    /// A block-vector (load-vector) corresponding to the right.hand side
362
    /// of the equation, filled during assembling
363
    std::shared_ptr<SystemVector> rhs;
364 365


366 367
  private: // some internal data-structures

368 369
    MatrixOperators<GlobalBasis> matrixOperators;
    VectorOperators<GlobalBasis> rhsOperators;
370

371
    Constraints<GlobalBasis> constraints;
372
  };
373 374


375
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
376
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
377
#endif
378

379
} // end namespace AMDiS
380 381 382

#include "ProblemStat.inc.hpp"