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

#include <cmath>

#include <amdis/LinearAlgebra.hpp>
6
#include <amdis/Observer.hpp>
7
#include <amdis/OperatorList.hpp>
8
9
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ConceptsBase.hpp>
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#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
26
  template <class RB, class CB, class T = double, class Traits = BackendTraits<RB>>
27
  class BiLinearForm
28
29
      : public MatrixFacade<T, Traits::template MatrixImpl>
      , private ObserverSequence<event::adapt,2>
30
31
  {
    using Self = BiLinearForm;
32
    using Super = MatrixFacade<T, Traits::template MatrixImpl>;
33
34
35
36
37
38
39
40
41

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

    /// The type of the finite element space / basis of the column
    using ColBasis = CB;

    /// 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)
55
      , symmetry_(SymmetryStructure::unknown)
Praetorius, Simon's avatar
Praetorius, Simon committed
56
57
      , rowBasis_(rowBasis)
      , colBasis_(colBasis)
58
      , updatePattern_(true)
59
    {
60
61
      auto const rowSize = rowBasis_->localView().maxSize();
      auto const colSize = colBasis_->localView().maxSize();
62
63
64
      elementMatrix_.resize(rowSize, colSize);
    }

65
66
67
68
69
70
71
72
73
74
    /// 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,
75
      REQUIRES(std::is_same_v<RB_,CB_>)>
76
77
78
79
80
81
82
83
84
85
86
    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)))
    {}

87
88
    RowBasis const& rowBasis() const { return *rowBasis_; }
    ColBasis const& colBasis() const { return *colBasis_; }
Praetorius, Simon's avatar
Praetorius, Simon committed
89

90
    /// \brief Associate a local operator with this BiLinearForm
91
92
93
94
95
96
97
98
    /**
     * 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
99
100
     *                     or \ref tag::boundary_operator indicating where to assemble
     *                     this operator.
101
     * \tparam Operator    A (pre-)operator that can be bound to a gridView, or a valid
102
103
104
     *                     GridOperator.
     * \tparam row         A tree-path for the RowBasis
     * \tparam col         A tree-path for the ColBasis
105
106
107
     *
     * [[expects: row is valid tree-path in RowBasis]]
     * [[expects: col is valid tree-path in ColBasis]]
108
     * @{
109
     **/
110
111
    template <class ContextTag, class Operator, class Row, class Col>
    void addOperator(ContextTag contextTag, Operator const& op, Row row, Col col);
112
113

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

    // Add an operator to be assembled on the intersections of the grid
122
123
    template <class Operator, class Row = RootTreePath, class Col = RootTreePath>
    void addIntersectionOperator(Operator const& op, Row row = {}, Col col = {})
124
    {
125
      using I = typename RowBasis::LocalView::GridView::Intersection;
126
      addOperator(tag::intersection_operator<I>{}, op, row, col);
127
128
    }
    /// @}
129

130
131
132
133
134
    /// Provide some additional information about the structure of the matrix pattern
    /**
     * If the symmetry structure is changed, it forces the matrix to
     * rebuild its patterns.
     **/
135
136
    void setSymmetryStructure(SymmetryStructure symm)
    {
137
      updatePattern_ = (symmetry_ != symm);
138
139
140
      symmetry_ = symm;
    }

141
    /// Set the symmetry structure using a string
142
143
144
145
    void setSymmetryStructure(std::string str)
    {
      setSymmetryStructure(symmetryStructure(str));
    }
146

147
    /// Assemble the matrix operators on the bound element.
148
    template <class RowLocalView, class ColLocalView, class LocalOperators,
149
150
      REQUIRES(Concepts::LocalView<RowLocalView>),
      REQUIRES(Concepts::LocalView<ColLocalView>)>
151
    void assemble(RowLocalView const& rowLocalView,
152
153
                  ColLocalView const& colLocalView,
                  LocalOperators& localOperators);
154

Praetorius, Simon's avatar
Praetorius, Simon committed
155
    /// Assemble all matrix operators, TODO: incorporate boundary conditions
156
    void assemble();
157

158
159
160
161
162
163
164
165
166
167
    /// Prepare the underlying matrix for insertion and rebuild the sparsity pattern
    /// if required.
    /**
     * The pattern is rebuild if either the parameter `forceRebuild` is set to true
     * or if for some other reasons the internal flag `updatePattern_` is set.
     * This might be the case if the symmetry structure is changed or the basis
     * is updated.
     *
     * If the pattern is not updated, just the values are set to zero.
     **/
168
    void init(bool forcePatternRebuild = false)
169
    {
170
171
172
173
174
175
176
177
      if (updatePattern_ || forcePatternRebuild)
      {
        Super::init(SparsityPattern{*rowBasis_, *colBasis_, symmetry_});
        updatePattern_ = false;
      }
      else
      {
        Super::init();
178
179
180
      }
    }

181
182
    /// Set the flag that forces an update of the pattern since the underlying
    /// basis that defines the indexset has been changed
183
184
185
186
187
188
189
190
191
192
193
194
195
    void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(*rowBasis_); }
    void updateImpl(event::adapt e, index_t<1> i) override { updateImpl2(*colBasis_); }

    auto& operators() { return operators_; }

  private:
    template <class Basis>
    void updateImpl2(Basis const& basis)
    {
      if (!updatePattern_)
        Recursive::forEach(operators_, [&](auto& op) { op.update(basis.gridView()); });
      updatePattern_ = true;
    }
196

197
  protected:
198
199
200
    /// The symmetry property if the bilinear form
    SymmetryStructure symmetry_;

201
202
203
204
205
    /// 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
206
207
208

    std::shared_ptr<RowBasis const> rowBasis_;
    std::shared_ptr<ColBasis const> colBasis_;
209
210

  private:
211
    // The pattern is rebuilt on calling init if this flag is set
212
    bool updatePattern_;
213
214
  };

215

216
  // deduction guides
217
218
  template <class RB, class CB>
  BiLinearForm(RB&&, CB&&)
219
    -> BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>;
220
221
222
223

  template <class RB>
  BiLinearForm(RB&&)
    -> BiLinearForm<Underlying_t<RB>, Underlying_t<RB>>;
224

225

226
  template <class T = double, class RB, class CB>
227
  auto makeBiLinearForm(RB&& rowBasis, CB&& colBasis)
228
  {
229
    using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<CB>, T>;
230
    return BLF{FWD(rowBasis), FWD(colBasis)};
231
232
  }

233
234
235
236
237
238
239
  template <class T = double, class RB>
  auto makeBiLinearForm(RB&& rowBasis)
  {
    using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<RB>, T>;
    return BLF{FWD(rowBasis)};
  }

240
241
242
} // end namespace AMDiS

#include <amdis/BiLinearForm.inc.hpp>