dune-fem  2.8-git
eigenmatrix.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_EIGENMATRIX_HH
2 #define DUNE_FEM_EIGENMATRIX_HH
3 
4 #ifdef HAVE_EIGEN
5 
6 // system includes
7 #include <iostream>
8 #include <string>
9 #include <utility>
10 
11 // local includes
13 #include <dune/fem/misc/functor.hh>
17 #include <dune/fem/io/parameter.hh>
23 #include <Eigen/Sparse>
24 
25 namespace Dune
26 {
27 
28  namespace Fem
29  {
30 
32  template <class T>
33  class EigenMatrix
34  {
35  static constexpr int defaultCol = -1;
36  static constexpr int firstCol = defaultCol + 1;
37 
38  public:
40  typedef T field_type;
42  typedef int size_type;
43  typedef Eigen::SparseMatrix<field_type,Eigen::RowMajor> MatrixStorageType;
44  typedef EigenMatrix<field_type> ThisType;
47  typedef ThisType MatrixBaseType;
48 
49  EigenMatrix(const ThisType& ) = delete;
50 
52  explicit EigenMatrix() :
53  matrix_()
54  {}
55 
58  EigenMatrix(size_type rows, size_type cols, size_type nz) :
59  matrix_(rows,cols)
60  {
61  reserve(rows,cols,nz);
62  }
63 
65  void reserve(size_type rows, size_type cols, size_type nz)
66  {
67  matrix_.resize(rows,cols);
68  matrix_.reserve(Eigen::VectorXi::Constant(rows,nz));;
69  }
70 
72  size_type rows() const
73  {
74  return matrix_.rows();
75  }
76 
78  size_type cols() const
79  {
80  return matrix_.cols();
81  }
82 
84  void set(size_type row, size_type col, field_type val)
85  {
86  matrix_.coeffRef(row,col) = val;
87  }
88 
90  void add(size_type row, size_type col, field_type val)
91  {
92  matrix_.coeffRef(row,col) += val;
93  }
94 
96  template<class ArgDFType, class DestDFType>
97  void apply(const ArgDFType& f, DestDFType& ret) const
98  {
99  ret.dofVector().array().coefficients() =
100  matrix_ * f.dofVector().array().coefficients();
101  }
102 
104  field_type operator()(size_type row, size_type col) const
105  {
106  return matrix_.coeff(row,col);
107  }
108 
110  void clear()
111  {
112  matrix_.setZero();
113  }
114 
116  void clearRow (size_type row)
117  {
118  matrix_.prune([row](const size_type& r, const size_type& c, const field_type& v)
119  { return (row != r); }); // || r == c); });
120  }
121 
122  template <class Vector>
123  void setUnitRows( const Vector &rows )
124  {
125  for (auto r : rows)
126  {
127  clearRow(r);
128  set(r,r,1.);
129  }
130  }
131 
134  size_type numNonZeros() const
135  {
136  return matrix_.nonZeros();
137  }
138 
141  size_type numNonZeros(size_type i) const
142  {
143  std::cout << "EigenMatrix::numNonZeros not yet implemented" << std::endl;
144  abort();
145  }
146 
149  std::pair<const field_type, size_type> realValue(size_type index) const
150  {
151  std::cout << "EigenMatrix::realValue not yet implemented" << std::endl;
152  abort();
153  }
154 
155  MatrixStorageType& data()
156  {
157  return matrix_;
158  }
159 
160  const MatrixStorageType& data() const
161  {
162  return matrix_;
163  }
164 
165  protected:
166  MatrixStorageType matrix_;
167  };
168 
169  template< class DomainSpace, class RangeSpace >
170  struct EigenMatrixObject
171  : public SparseRowMatrixObject< DomainSpace, RangeSpace, EigenMatrix< typename DomainSpace::RangeFieldType> >
172  {
173  typedef EigenMatrix< typename DomainSpace::RangeFieldType > MatrixType;
174  typedef SparseRowMatrixObject< DomainSpace, RangeSpace, MatrixType > BaseType;
175 
176  inline EigenMatrixObject( const DomainSpace &domainSpace,
177  const RangeSpace &rangeSpace,
178  const SolverParameter& param = SolverParameter() )
179  : BaseType( domainSpace, rangeSpace, param )
180  {}
181  };
182 
183  } // namespace Fem
184 
185 } // namespace Dune
186 
187 #endif
188 
189 #endif // #ifndef DUNE_FEM_SPMATRIX_HH
Definition: bindguard.hh:11