Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist über den Reiter "Standard" erreichbar. Die Administratoren

Dear Gitlab user, it is now possible to log in to our service using the ZIH login/LDAP. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

BiLinearForm.hpp 8.13 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)))
    {}

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

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>