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 15.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
#include <amdis/AdaptInfo.hpp>
#include <amdis/CreatorInterface.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
17
#include <amdis/CreatorMap.hpp>
18
#include <amdis/DirichletBC.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
19
//#include <amdis/Estimator.hpp>
20 21 22
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
23
#include <amdis/LinearSolvers.hpp>
24
#include <amdis/LocalAssemblerList.hpp>
25
#include <amdis/Marker.hpp>
26 27 28 29 30 31 32 33 34
#include <amdis/Mesh.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>

#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/Utility.hpp>

#include <amdis/GridFunctions.hpp>
35
#include <amdis/gridfunctions/DiscreteFunction.hpp>
36 37 38 39
#include <amdis/gridfunctions/DOFVectorView.hpp>

#include <amdis/io/FileWriterInterface.hpp>

40
#include <amdis/utility/TreeData.hpp>
41
#include <amdis/utility/TreePath.hpp>
42

43
namespace AMDiS
44
{
45 46 47 48
  // forward declaration
  template <class Traits>
  class ProblemInstat;

49
  template <class Traits>
50
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
51
      : public ProblemStatBase
52
      , public StandardProblemIteration
53
  {
54
    using Self = ProblemStat;
55

56 57
    friend class ProblemInstat<Traits>;

58
  public: // typedefs and static constants
59

60
    using GlobalBasis = typename Traits::GlobalBasis;
Praetorius, Simon's avatar
Praetorius, Simon committed
61 62 63
    using GridView    = typename GlobalBasis::GridView;
    using Grid        = typename GridView::Grid;
    using Element     = typename GridView::template Codim<0>::Entity;
64
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
65

66
    /// Dimension of the grid
67
    static constexpr int dim = Grid::dimension;
68

69
    /// Dimension of the world
70
    static constexpr int dow = Grid::dimensionworld;
71

72
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
73
    using SystemVector = DOFVector<GlobalBasis, double>;
74

75
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
76

77
  public:
78 79
    /**
     * \brief Constructor. Takes the name of the problem that is used to
80
     * access values corresponding to this problem in the parameter file.
81
     **/
82
    explicit ProblemStat(std::string const& name)
83
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
84
      , name_(name)
85
    {}
86

87 88
    /// Constructor taking additionally a reference to a grid that is used
    /// instead of the default created grid, \ref ProblemStat
89 90
    ProblemStat(std::string const& name, Grid& grid)
      : ProblemStat(name)
91
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
92
      adoptGrid(Dune::stackobject_to_shared_ptr(grid));
93 94 95 96
    }

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
97 98
    ProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
      : ProblemStat(name, grid)
99
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
100
      adoptGlobalBasis(Dune::stackobject_to_shared_ptr(globalBasis));
101
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
102

103

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

112

Praetorius, Simon's avatar
Praetorius, Simon committed
113
    /// Add an operator to \ref A.
114
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
    /// Operator evaluated on the whole element
    /**
     * Adds an operator to the list of element operators to be assembled in
     * quadrature points inside the element.
     *
     * \param op   A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator
     * \param row  TreePath identifying the sub-basis in the global basis tree
     *             corresponding to the row basis. \see treepath()
     * \param col  TreePath identifying the sub-basis in the global basis tree
     *             corresponding to the column basis. \see treepath()
     *
     * Example:
     * ```
     * auto op = makeOperator(tag::test_trial{}, 1.0/tau);
     * prob.addMatrixOperator(op, _0, _0);
     * ```
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
132
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
133 134 135 136
    void addMatrixOperator(Operator const& op, RowTreePath row = {}, ColTreePath col = {})
    {
      systemMatrix_->addOperator(tag::element_operator<Element>{}, op, row, col);
    }
137

Praetorius, Simon's avatar
Praetorius, Simon committed
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
    /// Operator evaluated on the boundary of the domain with boundary index `b`
    /**
     * Adds an operator to the list of boundary operators to be assembled in
     * quadrature points on the boundary intersections.
     *
     * \param b    Boundary indentifier where to assemble this operator. Can be
     *             constructed from an integer. \see BoundaryType
     * \param op   A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator
     * \param row  TreePath identifying the sub-basis in the global basis tree
     *             corresponding to the row basis. \see treepath()
     * \param col  TreePath identifying the sub-basis in the global basis tree
     *             corresponding to the column basis. \see treepath()
     *
     * Example:
     * ```
     * auto op = makeOperator(tag::test_trial{}, alpha);
     * prob.addMatrixOperator(BoundaryType{1}, op, _0, _0);
     * ```
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
157
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
158 159
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath row = {}, ColTreePath col = {})
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
160 161
      using I = typename GridView::Intersection;
      systemMatrix_->addOperator(tag::boundary_operator<I>{b}, op, row, col);
162
    }
163
    /** @} */
164 165


Praetorius, Simon's avatar
Praetorius, Simon committed
166
    /// Add an operator to \ref rhs.
167
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
    /// Operator evaluated on the whole element
    /**
     * Adds an operator to the list of element operators to be assembled in
     * quadrature points inside the element.
     *
     * \param op    A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator
     * \param path  TreePath identifying the sub-basis in the global basis tree
     *              corresponding to the row basis. \see treepath()
     *
     * Example:
     * ```
     * auto op = makeOperator(tag::test{}, probInstat.getOldSolution(0) / tau);
     * prob.addVectorOperator(op, _0);
     * ```
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
183
    template <class Operator, class TreePath = RootTreePath>
184 185 186 187
    void addVectorOperator(Operator const& op, TreePath path = {})
    {
      rhs_->addOperator(tag::element_operator<Element>{}, op, path);
    }
188

Praetorius, Simon's avatar
Praetorius, Simon committed
189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
    /// Operator evaluated on the boundary of the domain with boundary index `b`
    /**
     * Adds an operator to the list of boundary operators to be assembled in
     * quadrature points on the boundary intersections.
     *
     * \param b     Boundary indentifier where to assemble this operator. Can be
     *              constructed from an integer. \see BoundaryType
     * \param op    A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator
     * \param path  TreePath identifying the sub-basis in the global basis tree
     *              corresponding to the row basis. \see treepath()
     *
     * Example:
     * ```
     * auto op = makeOperator(tag::test{}, [g](auto const& x) { return g(x); });
     * prob.addVectorOperator(BoundaryType{1}, op, _0);
     * ```
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
206
    template <class Operator, class TreePath = RootTreePath>
207 208
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath path = {})
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
209 210
      using I = typename GridView::Intersection;
      rhs_->addOperator(tag::boundary_operator<I>{b}, op, path);
211
    }
212
    /** @} */
213

214

Praetorius, Simon's avatar
Praetorius, Simon committed
215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236
    /// Add boundary conditions to the system
    /** @{ */
    /// Dirichlet boundary condition
    /**
     * Enforce Dirichlet boundary values for the solution vector on boundary
     * regions identified by the predicate.
     *
     * \param predicate  Functor `bool(WorldVector)` returning true for all
     *                   DOFs on the boundary that should be assigned a value.
     * \param row        TreePath identifying the sub-basis in the global basis tree
     *                   corresponding to the row basis. \see treepath()
     * \param col        TreePath identifying the sub-basis in the global basis tree
     *                   corresponding to the column basis. \see treepath()
     * \param values     Functor `Range(WorldVector)` or any \ref GridFunction
     *                   that is evaluated in the DOFs identified by the predicate.
     *
     * Example:
     * ```
     * prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8; }, 0, 0,
     *                     [](auto const& x) { return 0.0; });
     * ```
     **/
237
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
238
    void addDirichletBC(Predicate const& predicate,
239
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
240
                        Values const& values);
Praetorius, Simon's avatar
Praetorius, Simon committed
241 242
    /** @} */

243

244
  public:
245

246
    /// Implementation of \ref ProblemStatBase::solve
247
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
248 249
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
250

251
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
252 253 254 255
    virtual void buildAfterAdapt(AdaptInfo& adaptInfo,
                                 Flag flag,
                                 bool asmMatrix = true,
                                 bool asmVector = true) override;
256

Praetorius, Simon's avatar
Praetorius, Simon committed
257 258 259 260 261 262 263
    /// \brief Assemble the linear system by calling \ref buildAfterAdapt with
    /// `asmMatrix` and `asmVector` set to true.
    void assemble(AdaptInfo& adaptInfo)
    {
      buildAfterAdapt(adaptInfo, Flag{0}, true, true);
    }

264
    /// Writes output files.
265
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
266

Praetorius, Simon's avatar
Praetorius, Simon committed
267

268
  public: // get-methods
269

Praetorius, Simon's avatar
Praetorius, Simon committed
270 271 272 273 274 275 276 277 278
    /// Implementation of \ref ProblemStatBase::name
    virtual std::string const& name() const override { return name_; }


    /// Return a reference to the grid, \ref grid
    Grid&       grid()       { return *grid_; }
    Grid const& grid() const { return *grid_; }

    /// Return the gridView of the leaf-level
279
    GridView const& gridView() const { return globalBasis_->gridView(); }
Praetorius, Simon's avatar
Praetorius, Simon committed
280 281 282 283 284 285 286 287 288

    /// Return the \ref globalBasis_
    GlobalBasis&       globalBasis()       { return *globalBasis_; }
    GlobalBasis const& globalBasis() const { return *globalBasis_; }

    /// Return a reference to the linear solver, \ref linearSolver
    LinearSolverType&       solver()       { return *linearSolver_; }
    LinearSolverType const& solver() const { return *linearSolver_; }

289
    /// Returns a reference to system-matrix, \ref systemMatrix_
Praetorius, Simon's avatar
Praetorius, Simon committed
290 291
    SystemMatrix&       systemMatrix()       { return *systemMatrix_; }
    SystemMatrix const& systemMatrix() const { return *systemMatrix_; }
292 293

    /// Returns a reference to the solution vector, \ref solution_
Praetorius, Simon's avatar
Praetorius, Simon committed
294 295
    SystemVector&       solutionVector()       { return *solution_; }
    SystemVector const& solutionVector() const { return *solution_; }
296 297

    /// Return a reference to the rhs system-vector, \ref rhs
Praetorius, Simon's avatar
Praetorius, Simon committed
298 299
    SystemVector&       rhsVector()       { return *rhs_; }
    SystemVector const& rhsVector() const { return *rhs_; }
300

301 302 303

    /// Return a mutable view to a solution component
    template <class TreePath = RootTreePath>
Praetorius, Simon's avatar
Praetorius, Simon committed
304
    auto solution(TreePath path = {})
305 306
    {
      auto&& tp = makeTreePath(path);
307
      return makeDOFVectorView(*solution_, tp);
308
    }
309

310 311
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
Praetorius, Simon's avatar
Praetorius, Simon committed
312
    auto solution(TreePath path = {}) const
313
    {
314
      auto&& tp = makeTreePath(path);
315
      return makeDiscreteFunction(*solution_, tp);
316
    }
317 318


Praetorius, Simon's avatar
Praetorius, Simon committed
319
  public: // set-methods
320

321
    /// Set a new linear solver for the problem
Praetorius, Simon's avatar
Praetorius, Simon committed
322
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
323
    {
324
      linearSolver_ = solver;
325 326
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
327 328 329 330 331
    void setSolver(LinearSolverType& solver)
    {
      setSolver(Dune::stackobject_to_shared_ptr(solver));
    }

332

333
    /// Set the grid. Stores pointer and initializes feSpaces
Praetorius, Simon's avatar
Praetorius, Simon committed
334
    /// matrices and vectors, as well as markers and file-writers.
335
    void setGrid(std::shared_ptr<Grid> const& grid)
336
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
337
      adoptGrid(grid);
338
      createGlobalBasis();
339
      createMatricesAndVectors();
Praetorius, Simon's avatar
Praetorius, Simon committed
340
      createMarker();
341 342 343
      createFileWriter();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
344 345 346 347 348
    void setGrid(Grid& grid)
    {
      setGrid(Dune::stackobject_to_shared_ptr(grid));
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
349 350 351 352 353 354 355 356 357 358 359 360
    void addMarker(std::shared_ptr<Marker<Grid>> const& marker)
    {
      marker_.push_back(marker);
      if (marker_.size() > 1)
        marker_.back()->setMaximumMarking(true);
    }

    void addMarker(Marker<Grid>& marker)
    {
      addMarker(Dune::stackobject_to_shared_ptr(marker));
    }

361

362
  protected: // initialization methods
363

Praetorius, Simon's avatar
Praetorius, Simon committed
364 365 366 367 368 369
    void createGlobalBasis();
    void createGrid();
    void createMatricesAndVectors();
    void createSolver();
    void createMarker();
    void createFileWriter();
370

Praetorius, Simon's avatar
Praetorius, Simon committed
371
    void adoptGlobalBasis(std::shared_ptr<GlobalBasis> const& globalBasis)
372
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
373
      globalBasis_ = globalBasis;
374 375 376
      initGlobalBasis(*globalBasis_);
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
377
    void adoptGrid(std::shared_ptr<Grid> const& grid)
378
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
379
      grid_ = grid;
380
      Parameters::get(name_ + "->mesh", gridName_);
381
    }
382

Praetorius, Simon's avatar
Praetorius, Simon committed
383
  private:
384

Praetorius, Simon's avatar
Praetorius, Simon committed
385 386
    void createGlobalBasisImpl(std::true_type);
    void createGlobalBasisImpl(std::false_type);
387

Praetorius, Simon's avatar
Praetorius, Simon committed
388
    void initGlobalBasis(GlobalBasis const& globalBasis);
389

Praetorius, Simon's avatar
Praetorius, Simon committed
390

391 392
  public: // implementation of iteration interface methods

Praetorius, Simon's avatar
Praetorius, Simon committed
393
    /// Implementation of \ref StandardProblemIteration::oneIteration.
394 395 396 397 398 399
    virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
    {
      return StandardProblemIteration::oneIteration(adaptInfo, toDo);
    }

    /// Implementation of \ref ProblemStatBase::estimate.
Praetorius, Simon's avatar
Praetorius, Simon committed
400
    virtual void estimate(AdaptInfo& adaptInfo) override { /* do nothing. */ }
401 402

    /// Implementation of \ref ProblemStatBase::refineMesh.
403
    virtual Flag adaptGrid(AdaptInfo& adaptInfo) override;
404 405

    /// Implementation of \ref ProblemStatBase::markElements.
406
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
407 408


409
  private:
410
    /// Name of this problem.
411
    std::string name_;
412

413
    /// Grid of this problem.
Praetorius, Simon's avatar
Praetorius, Simon committed
414
    std::shared_ptr<Grid> grid_;
415

416
    /// Name of the grid
417
    std::string gridName_ = "mesh";
418

419
    /// FE spaces of this problem.
420
    std::shared_ptr<GlobalBasis> globalBasis_;
421

422
    /// A FileWriter object
423
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
424

425
    /// Pointer to the adaptation markers
Praetorius, Simon's avatar
Praetorius, Simon committed
426
    std::list<std::shared_ptr<Marker<Grid>>> marker_;
Praetorius, Simon's avatar
Praetorius, Simon committed
427 428 429

    /// Pointer to the estimators for this problem
//    std::vector<Estimator*> estimator;
430

431
    /// An object of the linearSolver Interface
432
    std::shared_ptr<LinearSolverType> linearSolver_;
433

Praetorius, Simon's avatar
Praetorius, Simon committed
434
    /// Matrix that is filled during assembling
435
    std::shared_ptr<SystemMatrix> systemMatrix_;
436

Praetorius, Simon's avatar
Praetorius, Simon committed
437
    /// Vector with the solution components
438
    std::shared_ptr<SystemVector> solution_;
439

Praetorius, Simon's avatar
Praetorius, Simon committed
440
    /// Vector (load-vector) corresponding to the right-hand side
441
    /// of the equation, filled during assembling
442
    std::shared_ptr<SystemVector> rhs_;
443

Praetorius, Simon's avatar
Praetorius, Simon committed
444 445 446 447
    /// A vector with the local element error estimates
    /// for each node in the basis tree, indexed by [to_string(treePath)][element index]
    std::map<std::string, std::vector<double>> estimates_;

448

449 450
  private: // some internal data-structures

451
    Constraints<GlobalBasis, GlobalBasis> constraints_;
452
  };
453

Praetorius, Simon's avatar
Praetorius, Simon committed
454

455 456 457
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
  // Deduction rule
  template <class Grid, class GlobalBasis>
458
  ProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
459 460 461 462 463 464
    -> ProblemStat<DefaultProblemTraits<GlobalBasis>>;
#endif

  // Generator for ProblemStat with given Grid and GlobalBasis
  template <class Grid, class GlobalBasis>
  ProblemStat<DefaultProblemTraits<GlobalBasis>>
465
  makeProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
466
  {
467
    return {name, grid, globalBasis};
468 469
  }

470
} // end namespace AMDiS
471 472 473

#include "ProblemStat.inc.hpp"