DOFMatrixBase.hpp 5.77 KB
Newer Older
1
2
#pragma once

3
#include <cmath>
4
5
6

#include <dune/common/timer.hh>

7
#include <amdis/OperatorList.hpp>
8
#include <amdis/common/FlatMatrix.hpp>
9
#include <amdis/common/Math.hpp>
10
11
#include <amdis/typetree/MultiIndex.hpp>
#include <amdis/typetree/TreePath.hpp>
12

13
14
namespace AMDiS
{
15
16
17
18
19
20
21
22
23
  /**
   * 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 Backend  A linear-algebra backend for the matrix storage
   **/
24
  template <class RB, class CB, class Backend>
Praetorius, Simon's avatar
Praetorius, Simon committed
25
  class DOFMatrixBase
26
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
27
28
  public:
    /// The type of the finite element space / basis of the row
29
    using RowBasis = RB;
30
    using RowLocalView = typename RowBasis::LocalView;
Praetorius, Simon's avatar
Praetorius, Simon committed
31
32

    /// The type of the finite element space / basis of the column
33
    using ColBasis = CB;
34
35
36
37
    using ColLocalView = typename ColBasis::LocalView;

    using Element = typename RowLocalView::Element;
    using Geometry = typename Element::Geometry;
Praetorius, Simon's avatar
Praetorius, Simon committed
38
39
40

    /// The index/size - type
    using size_type = typename RowBasis::size_type;
41
    using value_type = typename Backend::value_type;
Praetorius, Simon's avatar
Praetorius, Simon committed
42

43
    /// The type of the data matrix used in the backend
44
    using BaseMatrix = typename Backend::BaseMatrix;
45
46

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

Praetorius, Simon's avatar
Praetorius, Simon committed
49
  public:
50
51
52
53
    /// Constructor. Stores the shared_ptr to the bases.
    DOFMatrixBase(std::shared_ptr<RowBasis> rowBasis, std::shared_ptr<ColBasis> colBasis)
      : rowBasis_(rowBasis)
      , colBasis_(colBasis)
54
    {
55
      operators_.init(*rowBasis_, *colBasis_);
56
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
57

58
59
60
61
62
63
    /// Constructor. Wraps the reference into a non-destroying shared_ptr or moves the basis into a new shared_ptr.
    template <class RB_, class CB_>
    DOFMatrixBase(RB_&& rowBasis, CB_&& colBasis)
      : DOFMatrixBase(Dune::wrap_or_move(FWD(rowBasis)), Dune::wrap_or_move(FWD(colBasis)))
    {}

Praetorius, Simon's avatar
Praetorius, Simon committed
64
    /// Return the row-basis \ref rowBasis of the matrix
65
    std::shared_ptr<RowBasis const> rowBasis() const
Praetorius, Simon's avatar
Praetorius, Simon committed
66
    {
67
      return rowBasis_;
Praetorius, Simon's avatar
Praetorius, Simon committed
68
69
70
    }

    /// Return the col-basis \ref colBasis of the matrix
71
    std::shared_ptr<ColBasis const> colBasis() const
Praetorius, Simon's avatar
Praetorius, Simon committed
72
    {
73
      return colBasis_;
Praetorius, Simon's avatar
Praetorius, Simon committed
74
75
    }

76
77
78
79
80
81
82
83
84
85
86
87
    /// Return the data-matrix
    BaseMatrix const& matrix() const
    {
      return backend_.matrix();
    }

    /// Return the data-matrix
    BaseMatrix& matrix()
    {
      return backend_.matrix();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
88
89
90
91
92
93
94
95
96
97
98
99
    /// Return the size of the \ref rowBasis_
    size_type rows() const
    {
      return rowBasis_->dimension();
    }

    /// Return the size of the \ref colBasis_
    size_type cols() const
    {
      return colBasis_->dimension();
    }

100
101
    /// Initialize the matrix for insertion, e.g. allocate the non-zero pattern
    /// If \p setToZero is true, the matrix is set to 0
102
    void init(bool asmMatrix);
103
104

    /// Finish the matrix insertion, e.g. cleanup or final insertion
105
    void finish(bool asmMatrix);
106
107

    /// Insert a block of values into the matrix (add to existing values)
108
    /// The global matrix indices are determined by the corresponding localviews.
109
110
    void insert(RowLocalView const& rowLocalView,
                ColLocalView const& colLocalView,
111
                ElementMatrix const& elementMatrix);
112
113
114
115
116
117
118
119

    /// Insert a single value into the matrix (add to existing value)
    template <class RowIndex, class ColIndex>
    void insert(RowIndex row, ColIndex col, typename Backend::value_type const& value)
    {
      backend_.insert(flatMultiIndex(row), flatMultiIndex(col), value);
    }

120
121
122
123
124
125
126
127
128
129
    /// \brief Associate a local operator with this DOFMatrix
    /**
     * 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
     *                     or \ref tag::boundary_operator indicating where to assemble this operator.
130
     * \tparam Expr        An pre-operator that can be bound to a gridView, or a valid
131
132
133
134
135
136
137
     *                      GridOperator.
     * \tparam row  A tree-path for the RowBasis
     * \tparam col  A tree-path for the ColBasis
     *
     * [[expects: row is valid tree-path in RowBasis]]
     * [[expects: col is valid tree-path in ColBasis]]
     **/
138
    // TODO: add method without contextTag.
139
    template <class ContextTag, class Expr,
140
              class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
141
    void addOperator(ContextTag contextTag, Expr const& expr,
142
                     RowTreePath row = {}, ColTreePath col = {});
143

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

148
149
    /// Assemble all matrix operators, TODO: incooperate boundary conditions
    void assemble();
150

151
    /// Number of nonzeros in the matrix
152
153
154
155
156
    size_type nnz() const
    {
      return backend_.nnz();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
157
  protected:
158
    /// The finite element space / basis associated with the rows
159
    std::shared_ptr<RowBasis> rowBasis_;
160
161

    /// The finite element space / basis associated with the columns
162
    std::shared_ptr<ColBasis> colBasis_;
163
164

    /// Data backend
165
    Backend backend_;
166

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

    /// List of operators associated to row/col node
171
    MatrixOperators<RowBasis,ColBasis,ElementMatrix> operators_;
172
173
174
  };

} // end namespace AMDiS
175
176

#include "DOFMatrixBase.inc.hpp"