DirichletBC.hpp 4.87 KB
Newer Older
1
2
3
#pragma once

#include <functional>
4
#include <list>
5
6
7
#include <type_traits>
#include <vector>

8
9
10
#include <dune/common/hybridutilities.hh>
#include <dune/functions/functionspacebases/interpolate.hh>

11
#include <amdis/Boundary.hpp>
12
13
14
15
16
#include <amdis/Output.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ValueCategory.hpp>
#include <amdis/utility/RangeType.hpp>
#include <amdis/utility/TreeData.hpp>
17

18
19
namespace AMDiS
{
20
21
  /// Implements a boundary condition of Dirichlet-type.
  /**
22
   * By calling the methods \ref init() and \ref finish before and after
23
   * assembling the system-matrix, respectively, dirichlet boundary conditions
24
   * can be applied to the matrix and system vector. Therefore, a predicate
25
   * functions indicates the DOFs where values should be enforced and a second
26
   * functor provided in the constructor is responsible for determining the
27
   * values to be set at the DOFs.
28
29
30
31
   *
   * 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
32
33
34
   * support this symmetric modification. As a result, this method returns a list
   * of columns values, that should be subtracted from the rhs.
   **/
35
  template <class Domain, class Range = double>
36
37
38
39
  class DirichletBC
  {
  public:
    template <class Predicate, class Values,
40
41
      REQUIRES(Concepts::Functor<Predicate, bool(Domain)> &&
				       Concepts::Functor<Values,  Range(Domain)>) >
42
    DirichletBC(Predicate&& predicate, Values&& values)
43
44
      : predicate_(std::forward<Predicate>(predicate))
      , values_(std::forward<Values>(values))
45
    {}
46
47


48
49
50
    /// 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.
51
52
53
54
55
    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*/)
    {
56
57
58
59
60
      if (!initialized_) {
        using RowNode = typename RowBasis::LocalView::Tree;
        initImpl(rowBasis, typename RowNode::NodeTag{});
        initialized_ = true;
      }
61
62
    }

63
64


65
66
67
68
69
70
71
    /// \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.
      **/
72
73
    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
    void finish(Matrix& matrix,
74
                VectorX& rhs,
75
76
                VectorB& solution,
                RowBasis const& rowBasis,
77
78
                ColBasis const& colBasis)
    {
79
80
81
      test_exit_dbg(initialized_, "Boundary condition not initialized!");
      auto columns = matrix.applyDirichletBC(dirichletNodes_);

82
83
84
      using Dune::Functions::interpolate;
      Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename RowBasis::LocalView::Tree>, Range>{},
        [&](auto id) {
85
          interpolate(id(rowBasis), rhs, values_, dirichletNodes_);
86
87
88
89
        });

      Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename ColBasis::LocalView::Tree>, Range>{},
        [&](auto id) {
90
          interpolate(id(colBasis), solution, values_, dirichletNodes_);
91
        });
92
93

      // subtract columns of dirichlet nodes from rhs
94
95
      // for (auto const& triplet : columns)
      //   rhs[triplet.row] -= triplet.value * solution[triplet.col];
96
      initialized_ = false;
97
    }
98
99
100
101
102
103
104
105

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

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

106
107
108
    template <class RowBasis>
    void initImpl(RowBasis const& rowBasis, Dune::TypeTree::CompositeNodeTag);

109
110
    template <class RB, class RowNodeTag>
    void initImpl(RB const&, RowNodeTag) {}
111

112
  private:
113
114
    std::function<bool(Domain)> predicate_;
    std::function<Range(Domain)> values_;
115

116
117
    bool initialized_ = false;
    std::vector<char> dirichletNodes_;
118
119
  };

120

121
  template <class Basis>
122
123
  struct DirichletData
  {
124
    using WorldVector = typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
125

126
    template <class RowNode, class ColNode>
127
    using type = std::list< std::shared_ptr<DirichletBC<WorldVector, RangeType_t<RowNode>>> >;
128
129
  };

130
131
  template <class RowBasis, class ColBasis>
  using Constraints = MatrixData<RowBasis, ColBasis, DirichletData<RowBasis>::template type>;
132

133
134
135
} // end namespace AMDiS

#include "DirichletBC.inc.hpp"