BiLinearForm.hpp 8.09 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 Expr        An 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 Expr, class Row, class Col>
    void addOperator(ContextTag contextTag, Expr const& expr, Row row, Col col);
112
113

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

    // Add an operator to be assembled on the intersections of the grid
122
123
    template <class Expr, class Row = RootTreePath, class Col = RootTreePath>
    void addIntersectionOperator(Expr const& expr, Row row = {}, Col col = {})
124
    {
125
      using I = typename RowBasis::LocalView::GridView::Intersection;
126
127
128
      addOperator(tag::intersection_operator<I>{}, expr, row, col);
    }
    /// @}
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
149
150
    template <class RowLocalView, class ColLocalView,
      REQUIRES(Concepts::LocalView<RowLocalView>),
      REQUIRES(Concepts::LocalView<ColLocalView>)>
151
152
153
    void assemble(RowLocalView const& rowLocalView,
                  ColLocalView const& colLocalView);

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

157
158
159
160
161
162
163
164
165
166
    /// 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.
     **/
167
    void init(bool forcePatternRebuild = false)
168
    {
169
170
171
172
173
174
175
176
      if (updatePattern_ || forcePatternRebuild)
      {
        Super::init(SparsityPattern{*rowBasis_, *colBasis_, symmetry_});
        updatePattern_ = false;
      }
      else
      {
        Super::init();
177
178
179
      }
    }

180
181
    /// Set the flag that forces an update of the pattern since the underlying
    /// basis that defines the indexset has been changed
182
183
    void updateImpl(event::adapt e, index_t<0> i) override { updatePattern_ = true; }
    void updateImpl(event::adapt e, index_t<1> i) override { updatePattern_ = true; }
184

185
  protected:
186
187
188
    /// The symmetry property if the bilinear form
    SymmetryStructure symmetry_;

189
190
191
192
193
    /// 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
194
195
196

    std::shared_ptr<RowBasis const> rowBasis_;
    std::shared_ptr<ColBasis const> colBasis_;
197
198

  private:
199
    // The pattern is rebuilt on calling init if this flag is set
200
    bool updatePattern_;
201
202
  };

203

204
  // deduction guides
205
206
  template <class RB, class CB>
  BiLinearForm(RB&&, CB&&)
207
    -> BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>;
208
209
210
211

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

213

214
  template <class T = double, class RB, class CB>
215
  auto makeBiLinearForm(RB&& rowBasis, CB&& colBasis)
216
  {
217
    using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<CB>, T>;
218
    return BLF{FWD(rowBasis), FWD(colBasis)};
219
220
  }

221
222
223
224
225
226
227
  template <class T = double, class RB>
  auto makeBiLinearForm(RB&& rowBasis)
  {
    using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<RB>, T>;
    return BLF{FWD(rowBasis)};
  }

228
229
230
} // end namespace AMDiS

#include <amdis/BiLinearForm.inc.hpp>