Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind ü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. The accounts of external users can be accessed via the "Standard" tab.
The administrators

ProblemStat.hpp 12.2 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 211 212 213 214 215 216 217 218 219 220 221
    std::shared_ptr<SystemVector> getSolutionVec() const
    {
      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);
    }
222

223 224 225
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
226
    {
227
      auto&& tp = makeTreePath(path);
228
      return makeDOFVectorView(*solution, tp);
229
    }
230 231


232
    /// Return a point to the rhs system-vector, \ref rhs
233 234
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
235 236


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

Praetorius, Simon's avatar
Praetorius, Simon committed
240
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
241
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
242
      linearSolver = solver;
243 244
    }

245 246
    /// Return a pointer to the grid, \ref grid
    std::shared_ptr<Grid> getGrid() { return grid; }
247

248 249
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
250
    void setGrid(Grid& grid_)
251
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
252
      grid = Dune::stackobject_to_shared_ptr(grid_);
253
      std::fill(componentGrids.begin(), componentGrids.end(), this->grid.get());
254

255
      createGlobalBasis();
256 257 258 259 260
      createMatricesAndVectors();
      createFileWriter();
    }


261
    /// Return the gridView of the leaf-level
262
    auto leafGridView() { return grid->leafGridView(); }
263 264

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

267
    /// Return the \ref feSpaces
268
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis; }
269 270


271
    /// Implementation of \ref ProblemStatBase::getName
272 273 274 275
    virtual std::string getName() const override
    {
      return name;
    }
276

277
    /// Return a vector of names of the problem components
278 279 280 281
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
282

283
    std::size_t getNumComponents() const
284 285 286 287
    {
      return nComponents;
    }

288
  protected: // initialization methods
289

290
    void createGrid()
291
    {
292 293 294
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
      test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
295

296
      grid = MeshCreator<Grid>::create(gridName);
297

298 299 300 301
      msg("Create grid:");
      msg("#elements = "   , grid->size(0));
      msg("#faces/edges = ", grid->size(1));
      msg("#vertices = "   , grid->size(dim));
302
      msg("");
303
    }
304

305
    void createGlobalBasis()
306
    {
307
      globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(leafGridView()));
308 309 310 311
      globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis->localView().tree());
      matrixOperators.init(*globalTree);
      rhsOperators.init(*globalTree);
      constraints.init(*globalTree);
312
    }
313

314 315
    void createMatricesAndVectors()
    {
316 317
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
      solution = std::make_shared<SystemVector>(*globalBasis, componentNames[0]);
318

319
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
320
      rhs = std::make_shared<SystemVector>(*globalBasis, rhsNames[0]);
321
    }
322

323
    void createSolver()
324
    {
325 326
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
327

328 329
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
330

331
      linearSolver = solverCreator->create(name + "->solver");
332
    }
333

334
    void createFileWriter();
335

Praetorius, Simon's avatar
Praetorius, Simon committed
336

337
  private:
338 339 340 341 342 343
    /// 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;
344

345 346
    /// Grid of this problem.
    std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
347

348
    /// Name of the mesh
349
    std::string gridName = "none";
350

351
    /// Pointer to the meshes for the different problem components
352
    std::vector<Grid*> componentGrids;
353

354
    /// FE spaces of this problem.
355
    std::shared_ptr<GlobalBasis> globalBasis;
356
    std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
357

358
    /// A FileWriter object
359
    std::list<std::shared_ptr<FileWriterInterface>> filewriter;
360

361
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
362
    std::shared_ptr<LinearSolverType> linearSolver;
363

364
    /// A block-matrix that is filled during assembling
365
    std::shared_ptr<SystemMatrix> systemMatrix;
366

367
    /// A block-vector with the solution components
368
    std::shared_ptr<SystemVector> solution;
369 370

    /// A block-vector (load-vector) corresponding to the right.hand side
371
    /// of the equation, filled during assembling
372
    std::shared_ptr<SystemVector> rhs;
373 374


375 376
  private: // some internal data-structures

377 378
    MatrixOperators<GlobalBasis> matrixOperators;
    VectorOperators<GlobalBasis> rhsOperators;
379

380
    Constraints<GlobalBasis> constraints;
381
  };
382 383


384
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
385
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
386
#endif
387

388
} // end namespace AMDiS
389 390 391

#include "ProblemStat.inc.hpp"