DirichletBC.hpp 3.98 KB
Newer Older
1
2
3
4
5
6
#pragma once

#include <functional>
#include <type_traits>
#include <vector>

7
#include <dune/amdis/Output.hpp>
8
9
#include <dune/amdis/common/Concepts.hpp>
#include <dune/amdis/utility/TreeData.hpp>
10

11
12
namespace AMDiS
{
13
14
  /// Implements a boundary condition of Dirichlet-type.
  /**
15
   * By calling the methods \ref init() and \ref finish before and after
16
   * assembling the system-matrix, respectively, dirichlet boundary conditions
17
   * can be applied to the matrix and system vector. Therefore, a predicate
18
   * functions indicates the DOFs where values should be enforced and a second
19
   * functor provided in the constructor is responsible for determining the
20
   * values to be set at the DOFs.
21
22
23
24
   *
   * In the \ref finish method the matrix is called with \ref applyDirichletBC
   * to erase the corresponding rows and columns for the DOF indices. This
   * application of boundary conditions can be symmetric if the matrix does
25
26
27
   * support this symmetric modification. As a result, this method returns a list
   * of columns values, that should be subtracted from the rhs.
   **/
28
  template <class Domain, class Range = double>
29
30
31
32
  class DirichletBC
  {
  public:
    template <class Predicate, class Values,
33
34
      REQUIRES(Concepts::Functor<Predicate, bool(Domain)> &&
				       Concepts::Functor<Values,  Range(Domain)>) >
35
36
    DirichletBC(Predicate&& predicate, Values&& values)
      : predicate(std::forward<Predicate>(predicate))
37
      , values(std::forward<Values>(values))
38
    {}
39
40


41
42
43
    /// Prepare the matrix, rhs and solution block for assembling dirichlet
    /// boundary conditions, e.g. collect the corresponding indices of boundary
    /// DOFS for later elimination in \ref finish.
44
45
46
47
48
49
50
51
52
    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
    void init(Matrix& /*matrix*/, VectorX& /*rhs*/, VectorB& /*solution*/,
              RowBasis const& rowBasis,
              ColBasis const& /*colBasis*/)
    {
      using RowNode = typename RowBasis::LocalView::Tree;
      initImpl(rowBasis, typename RowNode::NodeTag{});
    }

53
54


55
56
57
58
59
60
61
    /// \brief Apply dirichlet BC to matrix and vector, i.e., add a unit-row
    /// to the matrix and optionally delete the corresponding matrix-column.
    /** If DBC is set for matrix block {R,C}, then \p apply_row means that
      * the block-row R is assembled, and \p apply_col means that the block-col C
      * is assembled. The \p matrix, \p rhs, and \p solution correspond to the
      * currently visited blocks of system-matrix and system-vector.
      **/
62
63
    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
    void finish(Matrix& matrix,
64
                VectorX& rhs,
65
66
67
68
69
70
71
72
73
74
75
76
77
                VectorB& solution,
                RowBasis const& rowBasis,
                ColBasis const& colBasis);

  protected:
    template <class RowBasis>
    void initImpl(RowBasis const& rowBasis, Dune::TypeTree::LeafNodeTag);

    template <class RowBasis>
    void initImpl(RowBasis const& rowBasis, Dune::TypeTree::PowerNodeTag);

    template <class RB, class RowNodeTag>
    void initImpl(RB const&, RowNodeTag) {}
78

79
  private:
80
81
    std::function<bool(Domain)> predicate;
    std::function<Range(Domain)> values;
82

83
84
85
86
    bool initialized = false;
    std::vector<char> dirichletNodes;
  };

87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

  template <class GlobalBasis>
  struct DirichletData
  {
    using WorldVector = typename GlobalBasis::GridView::template Codim<0>::Geometry::GlobalCoordinate;

    template <class Node>
    struct type
    {
      std::list<std::shared_ptr<DirichletBC<WorldVector, double>>> scalar;
      std::list<std::shared_ptr<DirichletBC<WorldVector, WorldVector>>> vector;

      void push_back(std::shared_ptr<DirichletBC<WorldVector, double>> const& bc) { scalar.push_back(bc); }
      void push_back(std::shared_ptr<DirichletBC<WorldVector, WorldVector>> const& bc) { vector.push_back(bc); }
    };
  };

  template <class GlobalBasis>
  using Constraints = MatrixData<GlobalBasis, DirichletData<GlobalBasis>::template type>;

107
108
109
} // end namespace AMDiS

#include "DirichletBC.inc.hpp"