ProblemStat.hpp 14.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
#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 38
#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>
39

40
namespace AMDiS
41
{
42
  template <class Traits>
43
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
44
      : public ProblemStatBase
45
      , public StandardProblemIteration
46
  {
47
    using Self = ProblemStat;
48

49
  public: // typedefs and static constants
50

51
    using GlobalBasis = typename Traits::GlobalBasis;
52 53 54
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
    using Element  = typename GridView::template Codim<0>::Entity;
55
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
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 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 corresponding to this problem in the parameter file.
72
     **/
73
    explicit ProblemStat(std::string name)
74
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
75
      , name_(std::move(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(std::move(name))
82
    {
83 84 85 86 87 88 89 90 91 92 93
      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_);
94 95
    }

96

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

107 108

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

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


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

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

133

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

140
  public:
141

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


148 149 150 151 152 153 154 155 156 157 158 159
    /// 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;


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

166 167

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

Praetorius, Simon's avatar
Praetorius, Simon committed
170

171
  public: // get-methods
172

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

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

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

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


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


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

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

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

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

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


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

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

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


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

244
  protected: // initialization methods
245

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

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

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

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

    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
269 270 271 272
      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{});
273
    }
274

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

      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
        }
      }
288
    }
289

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

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

298
      linearSolver_ = solverCreator->create(name_ + "->solver");
299
    }
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 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347
    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++) {
        marker[i] = Marker<Grid>::
          createMarker(name + "->marker[" + std::to_string(i) + "]", i,
                       estimates, componentGrids[i]);

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

	        // If there is more than one marker, and all components are defined
	        // on the same mesh, the maximum marking has to be enabled.
 	        if (nMarkersCreated > 1 && nGrids == 1)
 	          marker[i]->setMaximumMarking(true);
        }
      }
    }


348
    void createFileWriter();
349

Praetorius, Simon's avatar
Praetorius, Simon committed
350

351 352 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
  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; }


386
  private:
387
    /// Name of this problem.
388
    std::string name_;
389

390
    /// Grid of this problem.
391
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
392

393 394 395
    /// Number of grids
    int nGrids = 1;

396
    /// Name of the mesh
397
    std::string gridName_ = "none";
398

399
    /// FE spaces of this problem.
400 401
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
402

403
    /// A FileWriter object
404
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
405

406 407 408 409 410 411
    /// Pointer to the adaptation markers
    std::vector<Marker<Grid>* > marker;

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

412
    /// An object of the linearSolver Interface
413
    std::shared_ptr<LinearSolverType> linearSolver_;
414

415
    /// A block-matrix that is filled during assembling
416
    std::shared_ptr<SystemMatrix> systemMatrix_;
417

418
    /// A block-vector with the solution components
419
    std::shared_ptr<SystemVector> solution_;
420

421 422 423
    /// A vector with the local element error estimates
    std::vector<std::vector<double> > estimates;

424
    /// A block-vector (load-vector) corresponding to the right.hand side
425
    /// of the equation, filled during assembling
426
    std::shared_ptr<SystemVector> rhs_;
427 428


429 430
  private: // some internal data-structures

431 432 433
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
434
  };
435

436 437 438 439 440 441 442 443 444 445 446 447 448 449 450
#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};
  }

451

452
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
453
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
454
#endif
455

456
} // end namespace AMDiS
457 458 459

#include "ProblemStat.inc.hpp"