DOFMatrix.hpp 3.87 KB
Newer Older
1
2
#pragma once

Praetorius, Simon's avatar
Praetorius, Simon committed
3
#include <list>
4
5
6
7
8
9
#include <string>
#include <memory>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>

10
11
#include <amdis/Output.hpp>
#include <amdis/common/ClonablePtr.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
12
#include <amdis/linear_algebra/Common.hpp>
13
14
#include <amdis/linear_algebra/DOFMatrixBase.hpp>
#include <amdis/utility/MultiIndex.hpp>
15
16
17

namespace AMDiS
{
18
19
20
21
22
23
24
25
26
27
28
29
30
31
  template <class T, class = void>
  struct BlockMatrixType
  {
    using type = Dune::FieldMatrix<T,1,1>;
  };

  template <class T>
  struct BlockMatrixType<T, typename T::field_type>
  {
    using type = T;
  };

  template <class RowBasisType, class ColBasisType,
            class ValueType = double>
32
  class DOFMatrix
Praetorius, Simon's avatar
Praetorius, Simon committed
33
      : public DOFMatrixBase<RowBasisType, ColBasisType>
34
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
35
    using Super = DOFMatrixBase<RowBasisType, ColBasisType>;
36

Praetorius, Simon's avatar
Praetorius, Simon committed
37
  public:
38
39
    /// The type of the elements of the DOFMatrix
    using value_type = typename BlockMatrixType<ValueType>::type;
40

41
42
    /// The matrix type of the underlying base matrix
    using BaseMatrix = Dune::BCRSMatrix<value_type>;
43

44
  public:
45
    /// Constructor. Constructs new BaseVector.
Praetorius, Simon's avatar
Praetorius, Simon committed
46
47
    DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis)
      : Super(rowBasis, colBasis)
48
      , matrix_(ClonablePtr<BaseMatrix>::make())
49
    {}
50

51
    /// Constructor. Takes reference to a base matrix.
Praetorius, Simon's avatar
Praetorius, Simon committed
52
53
    DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis, BaseMatrix& matrix)
      : Super(rowBasis, colBasis)
54
      , matrix_(matrix)
55
    {}
56

57
    /// Return the data-vector \ref vector
58
    BaseMatrix const& matrix() const
59
    {
60
      return *matrix_;
61
    }
62

63
    /// Return the data-vector \ref vector
64
    BaseMatrix& matrix()
65
    {
66
      return *matrix_;
67
    }
68
69


70
    /// Access the entry \p i of the \ref vector with write-access.
71
72
    template <class Index>
    value_type& operator()(Index row, Index col)
73
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
74
      auto r = flatMultiIndex(row), c = flatMultiIndex(col);
75
      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
Praetorius, Simon's avatar
Praetorius, Simon committed
76
77
      test_exit_dbg( r < this->rows() && c < this->cols() ,
          "Indices out of range [0,", this->rows(), ")x[0,", this->cols(), ")" );
78
      return (*matrix_)[r][c];
79
    }
80

81
    /// create occupation pattern and apply it to the matrix
82
    void init(bool prepareForInsertion)
83
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
84
      occupationPattern_.resize(this->rows(), this->cols());
85
86

      // A loop over all elements of the grid
Praetorius, Simon's avatar
Praetorius, Simon committed
87
88
      auto rowLocalView = this->rowBasis_->localView();
      auto colLocalView = this->colBasis_->localView();
89

Praetorius, Simon's avatar
Praetorius, Simon committed
90
      for (const auto& element : elements(this->rowBasis_->gridView())) {
91
92
        rowLocalView.bind(element);
        colLocalView.bind(element);
93

Praetorius, Simon's avatar
Praetorius, Simon committed
94
        for (std::size_t i = 0; i < rowLocalView.size(); ++i) {
95
          // The global index of the i-th vertex of the element
Praetorius, Simon's avatar
Praetorius, Simon committed
96
97
          auto row = rowLocalView.index(i);
          for (std::size_t j = 0; j < colLocalView.size(); ++j) {
98
            // The global index of the j-th vertex of the element
Praetorius, Simon's avatar
Praetorius, Simon committed
99
            auto col = colLocalView.index(j);
100
            occupationPattern_.add(row, col);
101
102
103
          }
        }
      }
104
      occupationPattern_.exportIdx(*matrix_);
105

106
      initialized_ = true;
107
108
109
110
    }

    void finish()
    {
111
      initialized_ = false;
112
    }
113

114
    std::size_t nnz() const
115
    {
116
      return matrix_->nonzeroes();
117
    }
118

119
    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
120
    {
121
      // loop over the matrix rows
122
      for (std::size_t i = 0; i < matrix_->N(); ++i) {
123
        if (dirichletNodes[i]) {
124
125
            auto cIt = (*matrix_)[i].begin();
            auto cEndIt = (*matrix_)[i].end();
126
127
            // loop over nonzero matrix entries in current row
            for (; cIt != cEndIt; ++cIt)
128
                *cIt = (i == cIt.index() ? 1 : 0);
129
130
        }
      }
131

Praetorius, Simon's avatar
Praetorius, Simon committed
132
      return std::list<Triplet<value_type>>{};
133
    }
134

135
  private:
136
137
    ClonablePtr<BaseMatrix> matrix_;
    Dune::MatrixIndexSet occupationPattern_;
138

139
    bool initialized_ = false;
140
141
142
  };

} // end namespace AMDiS