ProblemStat.hpp 14.3 KB
Newer Older
1 2
#pragma once

3 4 5 6
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

7 8 9 10 11 12
#include <tuple>
#include <string>
#include <memory>
#include <list>
#include <map>

13 14
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
15
#include <dune/istl/matrix.hh>
16
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
17 18
#include <dune/grid/common/grid.hh>

19
#include <dune/amdis/AdaptInfo.hpp>
20
#include <dune/amdis/Assembler.hpp>
21 22 23 24 25 26 27 28 29
#include <dune/amdis/CreatorInterface.hpp>
#include <dune/amdis/DirichletBC.hpp>
#include <dune/amdis/FileWriter.hpp>
#include <dune/amdis/Flag.hpp>
#include <dune/amdis/Initfile.hpp>
#include <dune/amdis/LinearAlgebra.hpp>
#include <dune/amdis/Mesh.hpp>
#include <dune/amdis/Operator.hpp>
#include <dune/amdis/ProblemStatBase.hpp>
30
#include <dune/amdis/ProblemStatTraits.hpp>
31 32
#include <dune/amdis/StandardProblemIteration.hpp>

33
#include <dune/amdis/common/OptionalDelete.hpp>
34 35
#include <dune/amdis/common/Timer.hpp>
#include <dune/amdis/common/TupleUtility.hpp>
36
#include <dune/amdis/common/TypeDefs.hpp>
37
#include <dune/amdis/common/Utility.hpp>
38

39
#include <dune/amdis/utility/RangeType.hpp>
40
#include <dune/amdis/utility/TreePath.hpp>
41

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

51
  public: // typedefs and static constants
52

53 54
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
55

56
    using Element  = typename GridView::template Codim<0>::Entity;
57

58
    static const int nComponents = 1; // TODO: count leaf nodes in GlobalBasis
59

60
    /// Dimension of the mesh
61
    static constexpr int dim = Grid::dimension;
62

63
    /// Dimension of the world
64
    static constexpr int dow = Grid::dimensionworld;
65

66
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
67

68 69
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
    using SystemVector = DOFVector<GlobalBasis, double>;
70

71 72
    using ElementOperator = Operator<GridView, Element>;
    using IntersectionOperator = Operator<GridView, typename GridView::Intersection>;
73

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

76
  public:
77 78
    /**
     * \brief Constructor. Takes the name of the problem that is used to
79
     * access values correpsonding to this püroblem in the parameter file.
80
     *
81
     * Parameters read by ProblemStat, with name 'PROB'
82
     *   PROB->names: \ref componentNames
83
     **/
84 85 86
    ProblemStat(std::string name)
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      , name(name)
87
    {
88
      Parameters::get(name + "->names", componentNames);
89
      for (std::size_t i = componentNames.size(); i < nComponents; ++i)
90
        componentNames.push_back("solution[" + std::to_string(i) + "]");
91
    }
92

93 94
    /// Constructor taking additionally a reference to a mesh that is used
    /// instead of the default created mesh, \ref ProblemStat
95
    ProblemStat(std::string name, Grid& grid)
96
      : ProblemStat(name)
97
    {
98
      this->grid = std::shared_ptr<Grid>(&grid, optional_delete(false));
99
      componentGrids.resize(nComponents, this->grid.get());
100

101 102
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
103 104
    }

105

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

116 117

    /// Adds an operator to \ref A.
118
    /** @{ */
119 120
    template <class RowTreePath, class ColTreePath>
    void addMatrixOperator(ElementOperator const& op, RowTreePath, ColTreePath,
121 122
                           double* factor = nullptr, double* estFactor = nullptr);

123 124
    template <class RowTreePath, class ColTreePath>
    void addMatrixOperator(IntersectionOperator const& op, RowTreePath, ColTreePath,
125
                           double* factor = nullptr, double* estFactor = nullptr);
126

127
    // operator evaluated on the boundary
128 129
    template <class RowTreePath, class ColTreePath>
    void addMatrixOperator(BoundaryType b, IntersectionOperator const& op, RowTreePath, ColTreePath,
130
                           double* factor = nullptr, double* estFactor = nullptr);
131

132 133 134 135 136 137 138
#if 0
    template <class RowTreePath, class ColTreePath>
    void addMatrixOperator(std::map< std::pair<RowTreePath,ColTreePath>, ElementOperator> ops);

    template <class RowTreePath, class ColTreePath>
    void addMatrixOperator(std::map< std::pair<RowTreePath,ColTreePath>, std::pair<ElementOperator, bool> > ops);
#endif
139
    /** @} */
140 141 142


    /// Adds an operator to \ref rhs.
143
    /** @{ */
144 145
    template <class TreePath>
    void addVectorOperator(ElementOperator const& op, TreePath,
146 147
                           double* factor = nullptr, double* estFactor = nullptr);

148 149
    template <class TreePath>
    void addVectorOperator(IntersectionOperator const& op, TreePath,
150
                           double* factor = nullptr, double* estFactor = nullptr);
151

152
    // operator evaluated on the boundary
153 154
    template <class TreePath>
    void addVectorOperator(BoundaryType b, IntersectionOperator const& op, TreePath,
155
                           double* factor = nullptr, double* estFactor = nullptr);
156

157 158 159 160 161 162 163
#if 0
    template <class TreePath>
    void addVectorOperator(std::map<TreePath, ElementOperator> ops);

    template <class TreePath>
    void addVectorOperator(std::map<TreePath, std::pair<ElementOperator, bool> > ops);
#endif
164
    /** @} */
165

166

167
    /// Adds a Dirichlet boundary condition
168
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
169
    void addDirichletBC(Predicate const& predicate,
170
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
171
                        Values const& values);
172

173 174 175 176 177 178 179 180
    template <class Predicate, class RowTreePath, class ColTreePath>
    void addDirichletBC(Predicate const& predicate,
                        RowTreePath row, ColTreePath col,
                        double constant)
    {
      addDirichletBC(predicate, row, col, [constant](auto const&) { return constant; });
    }

181

182
    /// Implementation of \ref ProblemStatBase::solve
183
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
184 185
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
186 187


188
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
189
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
190 191 192
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
193

194 195

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

Praetorius, Simon's avatar
Praetorius, Simon committed
198

199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
  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; }


240
  public: // get-methods
241

242
    /// Returns a pointer to system-matrix, \ref systemMatrix
243 244
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
245 246


247
    /// Return a pointer to the solution, \ref solution
248 249
    std::shared_ptr<SystemVector>       getSolution()       { return solution; }
    std::shared_ptr<SystemVector const> getSolution() const { return solution; }
250

251
    /// Return a view to a solution component
252 253
    template <class TreePath>
    auto getSolution(TreePath path)
254
    {
255 256 257
      using namespace Dune::Functions;

      auto tp = makeTreePath(path);
258
      using TP = decltype(tp);
259

260 261 262 263
      using Tree = typename GlobalBasis::LocalView::Tree;
      using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TP>;
      using Range = RangeType<Node>;
      using NTRM = Dune::Functions::DefaultNodeToRangeMap<Node>;
264

265 266 267 268
      using DiscreteFunction = Dune::Functions::DiscreteGlobalBasisFunction<GlobalBasis,TP,SystemVector,NTRM,Range>;

      auto nodeToRange = std::make_shared<NTRM>(makeDefaultNodeToRangeMap(*globalBasis, tp));
      return std::make_shared<DiscreteFunction>(globalBasis, tp, solution, nodeToRange);
269
    }
270 271


272
    /// Return a point to the rhs system-vector, \ref rhs
273 274
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
275 276


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

Praetorius, Simon's avatar
Praetorius, Simon committed
280
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
281
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
282
      linearSolver = solver;
283 284
    }

285 286
    /// Return a pointer to the grid, \ref grid
    std::shared_ptr<Grid> getGrid() { return grid; }
287

288 289
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
290
    void setGrid(Grid& grid_)
291
    {
292 293
      grid = std::shared_ptr<Grid>(&grid_, optional_delete(false));
      std::fill(componentGrids.begin(), componentGrids.end(), this->grid.get());
294

295
      createGlobalBasis();
296 297 298 299 300
      createMatricesAndVectors();
      createFileWriter();
    }


301
    /// Return the gridView of the leaf-level
302
    auto leafGridView() { return grid->leafGridView(); }
303 304

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

307
    /// Return the \ref feSpaces
308
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis; }
309 310


311
    /// Implementation of \ref ProblemStatBase::getName
312 313 314 315
    virtual std::string getName() const override
    {
      return name;
    }
316

317
    /// Return a vector of names of the problem components
318 319 320 321
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
322

323 324 325 326 327
    int getNumComponents() const
    {
      return nComponents;
    }

328
  protected: // initialization methods
329

330
    void createGrid()
331
    {
332 333 334
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
      test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
335

336
      grid = MeshCreator<Grid>::create(gridName);
337

338 339 340 341
      msg("Create grid:");
      msg("#elements = "   , grid->size(0));
      msg("#faces/edges = ", grid->size(1));
      msg("#vertices = "   , grid->size(dim));
342
      msg("");
343
    }
344

345
    void createGlobalBasis()
346
    {
347
      globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(leafGridView()));
348 349 350 351
      globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis->localView().tree());
      matrixOperators.init(*globalTree);
      rhsOperators.init(*globalTree);
      constraints.init(*globalTree);
352
    }
353

354 355
    void createMatricesAndVectors()
    {
356 357
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
      solution = std::make_shared<SystemVector>(*globalBasis, componentNames[0]);
358

359
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
360
      rhs = std::make_shared<SystemVector>(*globalBasis, rhsNames[0]);
361
    }
362

363
    void createSolver()
364
    {
365 366
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
367

368 369
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
370

371
      linearSolver = solverCreator->create(name + "->solver");
372
    }
373

374 375
    void createFileWriter()
    {
376
      filewriter = std::make_shared<FileWriter<GlobalBasis>>(name + "->output", leafGridView(), componentNames);
377
    }
378

Praetorius, Simon's avatar
Praetorius, Simon committed
379

380
  private:
381 382 383 384 385 386
    /// 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;
387

388 389
    /// Grid of this problem.
    std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
390

391
    /// Name of the mesh
392
    std::string gridName = "none";
393

394
    /// Pointer to the meshes for the different problem components
395
    std::vector<Grid*> componentGrids;
396

397
    /// FE spaces of this problem.
398
    std::shared_ptr<GlobalBasis> globalBasis; // eventuell const
399
    std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
400

401
    /// A FileWriter object
402
    std::shared_ptr<FileWriter<GlobalBasis>> filewriter;
403

404
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
405
    std::shared_ptr<LinearSolverType> linearSolver;
406

407
    /// A block-matrix that is filled during assembling
408
    std::shared_ptr<SystemMatrix> systemMatrix;
409

410
    /// A block-vector with the solution components
411
    std::shared_ptr<SystemVector> solution;
412 413

    /// A block-vector (load-vector) corresponding to the right.hand side
414
    /// of the equation, filled during assembling
415
    std::shared_ptr<SystemVector> rhs;
416 417


418 419
  private: // some internal data-structures

420 421
    MatrixOperators<GlobalBasis> matrixOperators;
    VectorOperators<GlobalBasis> rhsOperators;
422

423
    Constraints<GlobalBasis> constraints;
424
  };
425 426


427
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
428
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
429
#endif
430

431
} // end namespace AMDiS
432 433 434

#include "ProblemStat.inc.hpp"