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 14.9 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
#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/DirichletBC.hpp>
19
//#include <amdis/Estimator.hpp>
20 21 22
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
23
#include <amdis/Marker.cpp>
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 mesh
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))
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 96
    }

97

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

108 109

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

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


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

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

134

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

141
  public:
142

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


149 150 151 152 153 154 155 156 157 158 159 160
    /// Implementation of \ref ProblemStatBase::estimate.
    virtual void estimate(AdaptInfo& adaptInfo) override;


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


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


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

167 168

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

Praetorius, Simon's avatar
Praetorius, Simon committed
171

172
  public: // get-methods
173

174 175 176
    /// Returns a pointer to system-matrix, \ref systemMatrix_
    std::shared_ptr<SystemMatrix> getSystemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
177

178
    /// Returns a pointer to the solution vector, \ref solution_
179
    std::shared_ptr<SystemVector> getSolutionVector() const
180
    {
181
      return solution_;
182 183 184 185 186 187 188
    }

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

192 193 194
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
195
    {
196
      auto&& tp = makeTreePath(path);
197
      return makeDOFVectorView(*solution_, tp);
198
    }
199 200


201
    /// Return a point to the rhs system-vector, \ref rhs
202 203
    std::shared_ptr<SystemVector> getRhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
204 205


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

Praetorius, Simon's avatar
Praetorius, Simon committed
209
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
210
    {
211
      linearSolver_ = solver;
212 213
    }

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

217
    /// Set the grid. Stores pointer to passed reference and initializes feSpaces
218
    /// matrices and vectors, as well as the file-writer.
219
    void setGrid(Grid& grid)
220
    {
221 222
      grid_ = Dune::stackobject_to_shared_ptr(grid);

223
      createGlobalBasis();
224 225 226 227 228
      createMatricesAndVectors();
      createFileWriter();
    }


229
    /// Return the gridView of the leaf-level
230 231 232 233
    auto leafGridView() { return grid_->leafGridView(); }

    /// Return the gridView of levle `level`
    auto levelGridView(int level) { return grid_->levelGridView(level); }
234

235
    /// Return the \ref feSpaces
236
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
237 238


239
    /// Implementation of \ref ProblemStatBase::getName
240 241
    virtual std::string getName() const override
    {
242
      return name_;
243
    }
244

245
  protected: // initialization methods
246

247
    void createGrid()
248
    {
249 250
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
251
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
252

253
      grid_ = MeshCreator<Grid>::create(gridName_);
254

255
      msg("Create grid:");
256 257 258
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
259
      msg("");
260
    }
261

262
    void createGlobalBasis()
263
    {
264 265
      globalBasis_ = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid_->leafGridView()));
      initGlobalBasis(*globalBasis_);
266 267 268 269
    }

    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
270 271 272 273
      localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
      matrixOperators_.init(localView_->tree(), tag::store{});
      rhsOperators_.init(localView_->tree(), tag::store{});
      constraints_.init(localView_->tree(), tag::store{});
274
    }
275

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

      estimates.resize(leafGridView().indexSet().size(0));
      for (std::size_t i = 0; i < estimates.size(); i++) {
        estimates[i].resize(nComponents);
        for (std::size_t j = 0; j < estimates[i].size(); j++) {
          estimates[i][j] = 0.0; // TODO: Remove when estimate() is implemented
        }
      }
289
    }
290

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

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

299
      linearSolver_ = solverCreator->create(name_ + "->solver");
300
    }
301

302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332
    void createEstimator()
    {/*
      for (std::size_t i = 0, i < nComponents, i++) {
        std::string estName = "";
        Parameters::get(name + "->estimator->name[" + std::to_string(i) + "]", estName);
        auto estimatorCreator
          = named(CreatorMap<Estimator>::getCreator(estName, name + "->estimator->name[" + std::to_string(i) + "]"));
        if (estimatorCreator) {
	        estimatorCreator->setName(estName);
	        estimatorCreator->setRow(i);
	        estimatorCreator->setSolution(solution->getDOFVector(i));
	        estimator[i] = estimatorCreator->create();
        }

        if (estimator[i]) {
          for (std::size_t j = 0; j < nComponents; j++)
	          estimator[i]->addSystem((*systemMatrix)[i][j],
              solution->getDOFVector(j),
				      rhs->getDOFVector(i));
        }
      }*/
    }


    void createMarker()
    {
      int nMarkersCreated = 0;

      marker.resize(nComponents);

      for (std::size_t i = 0; i < nComponents; i++) {
333 334
        //treePath =
        marker[i] = Marker<Traits>:: // must use the same treePath as corresponding estimator
335
          createMarker(name + "->marker[" + std::to_string(i) + "]", i,
336
                       estimates[i]/*[treePath]*/, componentGrids[i]);
337 338 339 340 341

        if (marker[i]) {
	        nMarkersCreated++;

	        // If there is more than one marker, and all components are defined
342 343
	        // on the same grid, the maximum marking has to be enabled.
          // TODO: What about two markers each for two grids?
344 345 346 347 348 349 350
 	        if (nMarkersCreated > 1 && nGrids == 1)
 	          marker[i]->setMaximumMarking(true);
        }
      }
    }


351
    void createFileWriter();
352

Praetorius, Simon's avatar
Praetorius, Simon committed
353

354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388
  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.
    virtual Flag markElements(AdaptInfo& adaptInfo) override { return 0; }


389
  private:
390
    /// Name of this problem.
391
    std::string name_;
392

393
    /// Grid of this problem.
394
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
395

396 397 398
    /// Number of grids
    int nGrids = 1;

399
    /// Name of the grid
400
    std::string gridName_ = "none";
401

402
    /// FE spaces of this problem.
403 404
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
405

406
    /// A FileWriter object
407
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
408

409
    /// Pointer to the adaptation markers
410
    std::vector<Marker<Traits>* > marker;
411 412 413 414

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

415
    /// An object of the linearSolver Interface
416
    std::shared_ptr<LinearSolverType> linearSolver_;
417

418
    /// A block-matrix that is filled during assembling
419
    std::shared_ptr<SystemMatrix> systemMatrix_;
420

421
    /// A block-vector with the solution components
422
    std::shared_ptr<SystemVector> solution_;
423

424
    /// A vector with the local element error estimates
425
    /// reverse indexed by [component index][element index]
426 427
    std::vector<std::vector<double> > estimates;

428
    /// A block-vector (load-vector) corresponding to the right.hand side
429
    /// of the equation, filled during assembling
430
    std::shared_ptr<SystemVector> rhs_;
431 432


433 434
  private: // some internal data-structures

435 436 437
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
438 439

//    VectorData<GlobalBasis, TODO: NodeData> estimates_;
440
  };
441

442 443 444 445 446 447 448 449 450 451 452 453 454 455 456
#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};
  }

457

458
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
459
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
460
#endif
461

462
} // end namespace AMDiS
463 464 465

#include "ProblemStat.inc.hpp"