BiLinearForm.hpp 7.21 KB
Newer Older
1
2
3
4
5
#pragma once

#include <cmath>

#include <amdis/LinearAlgebra.hpp>
6
#include <amdis/Observer.hpp>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include <amdis/OperatorList.hpp>
#include <amdis/common/FlatMatrix.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
#include <amdis/typetree/TreePath.hpp>

namespace AMDiS
{
  /**
   * Basis implementation of DOFMatrix, i.e. a sparse matrix storing all the
   * assembled Operators indexed with DOF indices. The matrix data is associated
   * to a row and column global basis.
   *
   * \tparam RB  Basis of the matrix rows
   * \tparam CB  Basis of matrix columns
   * \tparam T   Coefficient type to store in the matrix
   **/
Praetorius, Simon's avatar
Praetorius, Simon committed
24
  template <class RB, class CB, class T = double, class Traits = BackendTraits<RB>>
25
  class BiLinearForm
26
27
      : public MatrixFacade<T, Traits::template MatrixImpl>
      , private ObserverSequence<event::adapt,2>
28
29
  {
    using Self = BiLinearForm;
30
    using Super = MatrixFacade<T, Traits::template MatrixImpl>;
31
32
33
34
35
36
37
38
39
40
41

  public:
    /// The type of the finite element space / basis of the row
    using RowBasis = RB;
    using RowLocalView = typename RowBasis::LocalView;

    /// The type of the finite element space / basis of the column
    using ColBasis = CB;
    using ColLocalView = typename ColBasis::LocalView;

    /// The type of the elements of the DOFVector
Praetorius, Simon's avatar
Praetorius, Simon committed
42
    using CoefficientType = T;
43
44
45
46

    /// The type of the matrix filled on an element with local contributions
    using ElementMatrix = FlatMatrix<CoefficientType>;

47
48
49
    /// Type of the sparsity pattern of the backend
    using SparsityPattern = typename Traits::SparsityPattern;

50
  public:
51
    /// Constructor. Stores the row and column basis in a local `shared_ptr` to const
Praetorius, Simon's avatar
Praetorius, Simon committed
52
53
    BiLinearForm(std::shared_ptr<RB> const& rowBasis, std::shared_ptr<CB> const& colBasis)
      : Super(*rowBasis, *colBasis)
54
      , ObserverSequence<event::adapt,2>(*rowBasis, *colBasis)
Praetorius, Simon's avatar
Praetorius, Simon committed
55
56
      , rowBasis_(rowBasis)
      , colBasis_(colBasis)
57
    {
58
      pattern_.init(*rowBasis_, *colBasis_);
59
      operators_.init(*rowBasis_, *colBasis_);
60

61
62
      auto const rowSize = rowBasis_->localView().maxSize();
      auto const colSize = colBasis_->localView().maxSize();
63
64
65
      elementMatrix_.resize(rowSize, colSize);
    }

66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
    /// Wraps the passed global bases into (non-destroying) shared_ptr
    template <class RB_, class CB_,
      REQUIRES(Concepts::Similar<RB_,RB>),
      REQUIRES(Concepts::Similar<CB_,CB>)>
    BiLinearForm(RB_&& rowBasis, CB_&& colBasis)
      : BiLinearForm(Dune::wrap_or_move(FWD(rowBasis)), Dune::wrap_or_move(FWD(colBasis)))
    {}

    /// Constructor for rowBasis == colBasis
    template <class RB_ = RB, class CB_ = CB,
      REQUIRES(std::is_same<RB_,CB_>::value)>
    explicit BiLinearForm(std::shared_ptr<RB> const& rowBasis)
      : BiLinearForm(rowBasis, rowBasis)
    {}

    /// Wraps the passed row-basis into a (non-destroying) shared_ptr
    template <class RB_,
      REQUIRES(Concepts::Similar<RB_,RB>)>
    explicit BiLinearForm(RB_&& rowBasis)
      : BiLinearForm(Dune::wrap_or_move(FWD(rowBasis)))
    {}

Praetorius, Simon's avatar
Praetorius, Simon committed
88
89
90
    std::shared_ptr<RowBasis const> const& rowBasis() const { return rowBasis_; }
    std::shared_ptr<ColBasis const> const& colBasis() const { return colBasis_; }

91
    /// \brief Associate a local operator with this BiLinearForm
92
93
94
95
96
97
98
99
    /**
     * Stores an operator in a list that gets assembled during a call to \ref assemble().
     * The operator may be assigned to a specific context, i.e. either an element
     * operator, an intersection operator, or a boundary operator.
     * The \p row and \p col tree paths specify the sub-basis for test and trial
     * functions the operator is applied to.
     *
     * \tparam ContextTag  One of \ref tag::element_operator, \ref tag::intersection_operator
100
101
     *                     or \ref tag::boundary_operator indicating where to assemble
     *                     this operator.
102
     * \tparam Expr        An pre-operator that can be bound to a gridView, or a valid
103
104
105
     *                     GridOperator.
     * \tparam row         A tree-path for the RowBasis
     * \tparam col         A tree-path for the ColBasis
106
107
108
     *
     * [[expects: row is valid tree-path in RowBasis]]
     * [[expects: col is valid tree-path in ColBasis]]
109
     * @{
110
     **/
111
112
    template <class ContextTag, class Expr, class Row, class Col>
    void addOperator(ContextTag contextTag, Expr const& expr, Row row, Col col);
113
114

    // Add an operator to be assembled on the elements of the grid
115
116
    template <class Expr, class Row = RootTreePath, class Col = RootTreePath>
    void addOperator(Expr const& expr, Row row = {}, Col col = {})
117
118
119
120
121
122
    {
      using E = typename RowLocalView::Element;
      addOperator(tag::element_operator<E>{}, expr, row, col);
    }

    // Add an operator to be assembled on the intersections of the grid
123
124
    template <class Expr, class Row = RootTreePath, class Col = RootTreePath>
    void addIntersectionOperator(Expr const& expr, Row row = {}, Col col = {})
125
126
127
128
129
    {
      using I = typename RowLocalView::GridView::Intersection;
      addOperator(tag::intersection_operator<I>{}, expr, row, col);
    }
    /// @}
130
131
132
133
134

    /// Assemble the matrix operators on the bound element.
    void assemble(RowLocalView const& rowLocalView,
                  ColLocalView const& colLocalView);

Praetorius, Simon's avatar
Praetorius, Simon committed
135
    /// Assemble all matrix operators, TODO: incorporate boundary conditions
136
137
    void assemble(SymmetryStructure symmetry = SymmetryStructure::unknown);

138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
    void init(SymmetryStructure symmetry = SymmetryStructure::unknown)
    {
      Super::init(pattern_, symmetry);
    }

    using Super::init;

    void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(e,i); }
    void updateImpl(event::adapt e, index_t<1> i) override { updateImpl2(e,i); }

  private:
    // Track for each basis wether updated and if both are, update the sparsity pattern
    template <std::size_t I>
    void updateImpl2(event::adapt, index_t<I>)
    {
      assert(!updateCounter_.test(I));
      updateCounter_.set(I);

      if (updateCounter_.all()) {
        pattern_.init(*rowBasis_, *colBasis_);
        updateCounter_.reset();
      }
    }

162
  protected:
163
164
165
    /// The structure of the non-zeros in the matrix
    SparsityPattern pattern_;

166
167
168
169
170
    /// Dense matrix to store coefficients during \ref assemble()
    ElementMatrix elementMatrix_;

    /// List of operators associated to row/col node
    MatrixOperators<RowBasis,ColBasis,ElementMatrix> operators_;
Praetorius, Simon's avatar
Praetorius, Simon committed
171
172
173

    std::shared_ptr<RowBasis const> rowBasis_;
    std::shared_ptr<ColBasis const> colBasis_;
174
175
176

  private:
    std::bitset<2> updateCounter_ = 0;
177
178
  };

179
180

#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
181
182
  template <class RB, class CB>
  BiLinearForm(RB&&, CB&&)
183
    -> BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>;
184
185
186
187

  template <class RB>
  BiLinearForm(RB&&)
    -> BiLinearForm<Underlying_t<RB>, Underlying_t<RB>>;
188
189
#endif

190
  template <class T = double, class RB, class CB>
191
  auto makeBiLinearForm(RB&& rowBasis, CB&& colBasis)
192
  {
193
    using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<CB>, T>;
194
    return BLF{FWD(rowBasis), FWD(colBasis)};
195
196
  }

197
198
199
200
201
202
203
  template <class T = double, class RB>
  auto makeBiLinearForm(RB&& rowBasis)
  {
    using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<RB>, T>;
    return BLF{FWD(rowBasis)};
  }

204
205
206
} // end namespace AMDiS

#include <amdis/BiLinearForm.inc.hpp>