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

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

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

13
#include <amdis/Boundary.hpp>
14
15
16
17
18
#include <amdis/Output.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ValueCategory.hpp>
#include <amdis/utility/RangeType.hpp>
#include <amdis/utility/TreeData.hpp>
19
#include <amdis/utility/Visitor.hpp>
20

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


51
52
53
    /// 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.
54
55
56
    template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
    void init(Matrix& /*matrix*/, VectorX& /*solution*/, VectorB& rhs,
              RowNode const& rowNode, ColNode const&)
57
    {
58
      if (!initialized_) {
59
60
61
62
        forEachLeafNode_(rowNode, [&](auto const& /*node*/, auto const& tp) {
          auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), cat(rowNode.treePath(),tp));
          this->initImpl(subBasis);
        });
63
64
        initialized_ = true;
      }
65
66
    }

67

68
69
    /// \brief Apply dirichlet BC to matrix and vector, i.e., add a unit-row
    /// to the matrix and optionally delete the corresponding matrix-column.
70
71
72
    template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
    void finish(Matrix& matrix, VectorX& solution, VectorB& rhs,
                RowNode const& rowNode, ColNode const& colNode)
73
    {
74
75
76
      test_exit_dbg(initialized_, "Boundary condition not initialized!");
      auto columns = matrix.applyDirichletBC(dirichletNodes_);

77
      Dune::Hybrid::ifElse(std::is_same<RangeType_t<RowNode>, Range>{},
78
        [&](auto id) {
79
80
          auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), rowNode.treePath());
          Dune::Functions::interpolate(subBasis, rhs, values_, dirichletNodes_);
81
82
        });

83
      Dune::Hybrid::ifElse(std::is_same<RangeType_t<ColNode>, Range>{},
84
        [&](auto id) {
85
86
          auto subBasis = Dune::Functions::subspaceBasis(solution.basis(), colNode.treePath());
          Dune::Functions::interpolate(subBasis, solution, values_, dirichletNodes_);
87
        });
88
89

      // subtract columns of dirichlet nodes from rhs
90
91
92
      for (auto const& triplet : columns)
        rhs[triplet.row] -= triplet.value * solution[triplet.col];

93
      initialized_ = false;
94
    }
95
96

  protected:
97
98
99
100
101
102
103
104
105
106
107
108
    // fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the boundary or not.
    template <class SubBasis>
    void initImpl(SubBasis const& subBasis)
    {
      using Dune::Functions::forEachBoundaryDOF;
      using LV = typename SubBasis::LocalView;
      using I = typename SubBasis::GridView::Intersection;
      dirichletNodes_.resize(subBasis.dimension());
      forEachBoundaryDOF(subBasis, [&](int localIndex, LV const& localView, I const& intersection) {
        dirichletNodes_[localView.index(localIndex)] = predicate_(intersection.geometry().center());
      });
    }
109

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

114
115
    bool initialized_ = false;
    std::vector<char> dirichletNodes_;
116
117
  };

118

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

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

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

131
} // end namespace AMDiS