DOFMatrixBase.hpp 5.35 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/Math.hpp>
9 10
#include <amdis/typetree/MultiIndex.hpp>
#include <amdis/typetree/TreePath.hpp>
11

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

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

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

    /// The index/size - type
    using size_type = typename RowBasis::size_type;

41
    using BaseMatrix = typename Backend::BaseMatrix;
42
    using ElementMatrix = Dune::DynamicMatrix<double>;
43

Praetorius, Simon's avatar
Praetorius, Simon committed
44 45 46 47 48
  public:
    /// Constructor. Constructs new BaseVector.
    DOFMatrixBase(RowBasis const& rowBasis, ColBasis const& colBasis)
      : rowBasis_(&rowBasis)
      , colBasis_(&colBasis)
49
    {
50
      operators_.init(rowBasis, colBasis);
51
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
52

53 54 55
    DOFMatrixBase(DOFMatrixBase const&) = delete;
    DOFMatrixBase(DOFMatrixBase&&) = delete;

Praetorius, Simon's avatar
Praetorius, Simon committed
56 57 58 59 60 61 62 63 64 65 66 67
    /// Return the row-basis \ref rowBasis of the matrix
    RowBasis const& rowBasis() const
    {
      return *rowBasis_;
    }

    /// Return the col-basis \ref colBasis of the matrix
    ColBasis const& colBasis() const
    {
      return *colBasis_;
    }

68 69 70 71 72 73 74 75 76 77 78 79
    /// 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
80 81 82 83 84 85 86 87 88 89 90 91
    /// 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();
    }

92 93
    /// Initialize the matrix for insertion, e.g. allocate the non-zero pattern
    /// If \p setToZero is true, the matrix is set to 0
94
    void init(bool asmMatrix);
95 96

    /// Finish the matrix insertion, e.g. cleanup or final insertion
97
    void finish(bool asmMatrix);
98 99

    /// Insert a block of values into the matrix (add to existing values)
100 101
    void insert(RowLocalView const& rowLocalView,
                ColLocalView const& colLocalView,
102
                ElementMatrix const& elementMatrix);
103 104 105 106 107 108 109 110

    /// 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);
    }

111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
    /// \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.
     * \tparam PreOperator  An pre-operator that can be bound to a gridView, or a valid
     *                      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]]
     **/
129 130 131
    template <class ContextTag, class Operator,
              class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addOperator(ContextTag contextTag, Operator const& preOp,
132
                     RowTreePath row = {}, ColTreePath col = {});
133

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

138 139
    /// Assemble all matrix operators, TODO: incooperate boundary conditions
    void assemble();
140

141 142 143 144 145 146 147 148 149 150 151 152
    /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
    /// a one on the diagonal of the matrix.
    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
    {
      return backend_.applyDirichletBC(dirichletNodes);
    }

    size_type nnz() const
    {
      return backend_.nnz();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
153
  protected:
154 155 156
    /// The finite element space / basis associated with the rows/columns
    RowBasis const* rowBasis_;
    ColBasis const* colBasis_;
157 158

    /// Data backend
159
    Backend backend_;
160

161 162 163 164 165
    /// Dense matrix to store coefficients during \ref assemble()
    ElementMatrix elementMatrix_;

    /// List of operators associated to row/col node
    MatrixOperators<RowBasis,ColBasis> operators_;
166 167 168
  };

} // end namespace AMDiS
169 170

#include "DOFMatrixBase.inc.hpp"