Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist über den Reiter "Standard" erreichbar. Die Administratoren

Dear Gitlab user, it is now possible to log in to our service using the ZIH login/LDAP. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

ProblemStat.hpp 11 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 84
      gridName = "";
      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
      globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis->localView().tree());
      matrixOperators.init(*globalTree);
      rhsOperators.init(*globalTree);
      constraints.init(*globalTree);
284
    }
285

286 287
    void createMatricesAndVectors()
    {
288
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
289 290
      solution = std::make_shared<SystemVector>(*globalBasis, "solution");
      rhs = std::make_shared<SystemVector>(*globalBasis, "rhs");
291
    }
292

293
    void createSolver()
294
    {
295 296
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
297

298 299
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
300

301
      linearSolver = solverCreator->create(name + "->solver");
302
    }
303

304
    void createFileWriter();
305

Praetorius, Simon's avatar
Praetorius, Simon committed
306

307
  private:
308 309 310
    /// Name of this problem.
    std::string name;

311 312
    /// Grid of this problem.
    std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
313

314
    /// Name of the mesh
315
    std::string gridName = "none";
316

317
    /// FE spaces of this problem.
318
    std::shared_ptr<GlobalBasis> globalBasis;
319
    std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
320

321
    /// A FileWriter object
322
    std::list<std::shared_ptr<FileWriterInterface>> filewriter;
323

324
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
325
    std::shared_ptr<LinearSolverType> linearSolver;
326

327
    /// A block-matrix that is filled during assembling
328
    std::shared_ptr<SystemMatrix> systemMatrix;
329

330
    /// A block-vector with the solution components
331
    std::shared_ptr<SystemVector> solution;
332 333

    /// A block-vector (load-vector) corresponding to the right.hand side
334
    /// of the equation, filled during assembling
335
    std::shared_ptr<SystemVector> rhs;
336 337


338 339
  private: // some internal data-structures

340 341
    MatrixOperators<GlobalBasis> matrixOperators;
    VectorOperators<GlobalBasis> rhsOperators;
342

343
    Constraints<GlobalBasis> constraints;
344
  };
345 346


347
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
348
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
349
#endif
350

351
} // end namespace AMDiS
352 353 354

#include "ProblemStat.inc.hpp"