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 13.6 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
#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>
22
#include <amdis/LocalAssemblerList.hpp>
23
#include <amdis/Marker.hpp>
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))
Praetorius, Simon's avatar
Praetorius, Simon committed
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
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
96 97 98 99 100 101 102 103 104 105 106 107 108

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
    ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      , name(std::move(name))
    {
      gridName = "";
      Parameters::get(name + "->mesh", gridName);

      globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
      initGlobalBasis(*globalBasis_);
    }
109

110

111
    /**
112
     * \brief Initialisation of the problem.
113
     *
114 115 116
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
117
    void initialize(Flag initFlag,
118
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
119
                    Flag adoptFlag = INIT_NOTHING);
120

121 122

    /// Adds an operator to \ref A.
123
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
124 125
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath = {}, ColTreePath = {});
126

127
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
128 129
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath = {}, ColTreePath = {});
130
    /** @} */
131 132 133


    /// Adds an operator to \ref rhs.
134
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
135 136
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(Operator const& op, TreePath = {});
137

138
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
139 140
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath = {});
141
    /** @} */
142

143

144
    /// Adds a Dirichlet boundary condition
145
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
146
    void addDirichletBC(Predicate const& predicate,
147
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
148
                        Values const& values);
149

150
  public:
151

152
    /// Implementation of \ref ProblemStatBase::solve
153
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
154 155
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
156

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

163
    /// Writes output files.
164
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
165

Praetorius, Simon's avatar
Praetorius, Simon committed
166

167
  public: // get-methods
168

169 170 171
    /// Returns a pointer to system-matrix, \ref systemMatrix_
    std::shared_ptr<SystemMatrix> getSystemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
172

173
    /// Returns a pointer to the solution vector, \ref solution_
174
    std::shared_ptr<SystemVector> getSolutionVector() const
175
    {
176
      return solution_;
177 178 179 180 181 182 183
    }

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

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


196
    /// Return a point to the rhs system-vector, \ref rhs
197 198
    std::shared_ptr<SystemVector> getRhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
199 200


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

Praetorius, Simon's avatar
Praetorius, Simon committed
204
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
205
    {
206
      linearSolver_ = solver;
207 208
    }

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

Praetorius, Simon's avatar
Praetorius, Simon committed
212
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
213
    /// matrices and vectors, as well as the file-writer.
214
    void setGrid(Grid& grid)
215
    {
216 217
      grid_ = Dune::stackobject_to_shared_ptr(grid);

218
      createGlobalBasis();
219 220 221 222
      createMatricesAndVectors();
      createFileWriter();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
223 224
    /// Return the gridView of the leaf-level
    auto const& gridView() { return globalBasis_->gridView(); }
225

226
    /// Return the \ref feSpaces
227
    std::shared_ptr<GlobalBasis> const& globalBasis() { return globalBasis_; }
228 229


230
    /// Implementation of \ref ProblemStatBase::getName
231 232
    virtual std::string getName() const override
    {
233
      return name_;
234
    }
235

236
  protected: // initialization methods
237

238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258
    template <class T, class GV>
    using HasCreate = decltype(T::create(std::declval<GV>()));

    void createGlobalBasis()
    {
      createGlobalBasis(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
      initGlobalBasis(*globalBasis_);
    }

    void createGlobalBasis(std::true_type)
    {
      assert( bool(grid) );
      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid->leafGridView()));
    }

    void createGlobalBasis(std::false_type)
    {
      error_exit("Can not create GlobalBasis from type. Pass a BasisCreator instead!");
    }

259
    void createGrid()
260
    {
261 262
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
263
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
264

265
      grid_ = MeshCreator<Grid>::create(gridName_);
266

267
      msg("Create grid:");
268 269 270
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
271
      msg("");
272
    }
273

274 275
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
276
      localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
Praetorius, Simon's avatar
Praetorius, Simon committed
277 278 279
      matrixOperators.init(localView_->tree(), tag::store{});
      rhsOperators.init(localView_->tree(), tag::store{});
      constraints.init(localView_->tree(), tag::store{});
280
    }
281

282 283
    void createMatricesAndVectors()
    {
284 285 286
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
287

288
      AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
289 290
      {
        std::string i = to_string(treePath);
291 292 293
        estimates_[i].resize(globalBasis_->gridView().indexSet().size(0));
        for (std::size_t j = 0; j < estimates_[i].size(); j++)
          estimates_[i][j] = 0.0; // TODO: Remove when estimate() is implemented
294 295 296
      });

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

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

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

307
      linearSolver_ = solverCreator->create(name_ + "->solver");
308
    }
309

310
    void createMarker();
311

312
    void createFileWriter();
313

Praetorius, Simon's avatar
Praetorius, Simon committed
314

315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
  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.
347
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
348 349


350
  private:
351
    /// Name of this problem.
352
    std::string name_;
353

354
    /// Grid of this problem.
355
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
356

357 358 359
    /// Number of grids
    int nGrids = 1;

360
    /// Name of the grid
361
    std::string gridName_ = "none";
362

363
    /// FE spaces of this problem.
364 365
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
366

367
    /// A FileWriter object
368
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
369

370
    /// Pointer to the adaptation markers
371
    std::list<std::shared_ptr<Marker<Traits> > > marker_;
372

373
    /// An object of the linearSolver Interface
374
    std::shared_ptr<LinearSolverType> linearSolver_;
375

376
    /// A block-matrix that is filled during assembling
377
    std::shared_ptr<SystemMatrix> systemMatrix_;
378

379
    /// A block-vector with the solution components
380
    std::shared_ptr<SystemVector> solution_;
381

382
    /// A vector with the local element error estimates
383 384
    /// for each node in the basis tree, indexed by [to_string(treePath)][element index]
    std::map<std::string, std::vector<double> > estimates_;
385

386
    /// A block-vector (load-vector) corresponding to the right.hand side
387
    /// of the equation, filled during assembling
388
    std::shared_ptr<SystemVector> rhs_;
389 390


391 392
  private: // some internal data-structures

393 394 395
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
396
  };
397

398 399 400 401 402 403 404 405 406 407 408 409 410 411 412
#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};
  }

413

414
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
415
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
416
#endif
417

418
} // end namespace AMDiS
419 420 421

#include "ProblemStat.inc.hpp"