dune-localfunctions  2.2.1
basismatrix.hh
Go to the documentation of this file.
1 #ifndef DUNE_BASISMATRIX_HH
2 #define DUNE_BASISMATRIX_HH
3 
4 #include <fstream>
5 #include <dune/common/exceptions.hh>
6 
10 
11 namespace Dune
12 {
13  /****************************************
14  * A dense matrix representation of a ''polynomial''
15  * basis. Its represent a basis as a linear
16  * combination of a second basis, i.e., a
17  * monomial basis. It is simular to the PolynomialBasis
18  * but it not derived from the LocalBasis class.
19  * It is used to define a ''pre basis''.
20  ****************************************/
21  template< class PreBasis, class Interpolation,
22  class Field >
23  struct BasisMatrix;
24 
25  template< class PreBasis, class Interpolation,
26  class Field >
27  struct BasisMatrixBase : public LFEMatrix<Field>
28  {
30 
31  BasisMatrixBase( const PreBasis& preBasis,
32  const Interpolation& localInterpolation )
33  : cols_(preBasis.size())
34  {
35  localInterpolation.interpolate( preBasis, *this );
36 
37  if ( !Matrix::invert() )
38  {
39  DUNE_THROW(MathError, "While computing basis a singular matrix was constructed!");
40  }
41  }
42  unsigned int cols () const
43  {
44  return cols_;
45  }
46  unsigned int rows () const
47  {
48  return Matrix::rows();
49  }
50  private:
51  unsigned int cols_;
52  };
53 
54  template< class Topology, class F,
55  class Interpolation,
56  class Field >
57  struct BasisMatrix< const MonomialBasis< Topology, F >, Interpolation, Field >
58  : public BasisMatrixBase< const MonomialBasis< Topology, F >, Interpolation, Field >
59  {
62  typedef typename Base::Matrix Matrix;
63 
64  BasisMatrix( const PreBasis& preBasis,
65  const Interpolation& localInterpolation )
66  : Base(preBasis, localInterpolation)
67  {
68  }
69  template <class Vector>
70  void row( const unsigned int row, Vector &vec ) const
71  {
72  const unsigned int N = Matrix::rows();
73  assert( Matrix::cols() == N && vec.size() == N );
74  // note: that the transposed matrix is computed,
75  // and is square
76  for (unsigned int i=0;i<N;++i)
77  field_cast(Matrix::operator()(i,row),vec[i]);
78  }
79  };
80  template< int dim, class F,
81  class Interpolation,
82  class Field >
83  struct BasisMatrix< const Dune::VirtualMonomialBasis< dim, F >, Interpolation, Field >
84  : public BasisMatrixBase< const VirtualMonomialBasis< dim, F >, Interpolation, Field >
85  {
88  typedef typename Base::Matrix Matrix;
89 
90  BasisMatrix( const PreBasis& preBasis,
91  const Interpolation& localInterpolation )
92  : Base(preBasis, localInterpolation)
93  {
94  }
95  template <class Vector>
96  void row( const unsigned int row, Vector &vec ) const
97  {
98  const unsigned int N = Matrix::rows();
99  assert( Matrix::cols() == N && vec.size() == N );
100  // note: that the transposed matrix is computed,
101  // and is square
102  for (unsigned int i=0;i<N;++i)
103  field_cast(Matrix::operator()(i,row),vec[i]);
104  }
105  };
106  template< class Eval, class CM, class D, class R,
107  class Interpolation,
108  class Field >
109  struct BasisMatrix< const PolynomialBasis<Eval,CM,D,R>, Interpolation, Field >
110  : public BasisMatrixBase< const PolynomialBasis<Eval,CM,D,R>, Interpolation, Field >
111  {
114  typedef typename Base::Matrix Matrix;
115 
116  BasisMatrix( const PreBasis& preBasis,
117  const Interpolation& localInterpolation )
118  : Base(preBasis, localInterpolation),
119  preBasis_(preBasis)
120  {
121  }
122  unsigned int cols() const
123  {
124  return preBasis_.matrix().baseSize() ;
125  }
126  template <class Vector>
127  void row( const unsigned int row, Vector &vec ) const
128  {
129  assert( Matrix::rows() == Matrix::cols() );
130  assert( vec.size() == preBasis_.matrix().baseSize() );
131  assert( Matrix::cols() == preBasis_.size() );
132  for (unsigned int j=0;j<Matrix::cols();++j)
133  vec[j] = 0;
134  for (unsigned int i=0;i<Matrix::rows();++i)
135  preBasis_.matrix().
136  addRow(i,Base::Matrix::operator()(i,row),vec);
137  }
138  private:
139  const PreBasis& preBasis_;
140  };
141  template< class Eval, class CM,
142  class Interpolation,
143  class Field >
144  struct BasisMatrix< const PolynomialBasisWithMatrix<Eval,CM>, Interpolation, Field >
145  : public BasisMatrixBase< const PolynomialBasisWithMatrix<Eval,CM>, Interpolation, Field >
146  {
149  typedef typename Base::Matrix Matrix;
150 
151  BasisMatrix( const PreBasis& preBasis,
152  const Interpolation& localInterpolation )
153  : Base(preBasis, localInterpolation),
154  preBasis_(preBasis)
155  {
156  }
157  unsigned int cols() const
158  {
159  return preBasis_.matrix().baseSize() ;
160  }
161  unsigned int rows () const
162  {
163  assert( Matrix::rows() == preBasis_.matrix().size() );
164  return preBasis_.matrix().size()*CM::blockSize ;
165  }
166  template <class Vector>
167  void row( const unsigned int row, Vector &vec ) const
168  {
169  unsigned int r = row / CM::blockSize;
170  assert( r < Matrix::rows() );
171  assert( Matrix::rows() == Matrix::cols() );
172  assert( vec.size() == preBasis_.matrix().baseSize() );
173  assert( Matrix::cols() == preBasis_.size() );
174  for (unsigned int j=0;j<vec.size();++j)
175  vec[j] = 0;
176  for (unsigned int i=0;i<Matrix::rows();++i)
177  preBasis_.matrix().
178  addRow(i*CM::blockSize+row%CM::blockSize,Base::Matrix::operator()(i,r),vec);
179  }
180  private:
181  const PreBasis& preBasis_;
182  };
183 }
184 
185 #endif // DUNE_BASISMATRIX_HH
186