dune-localfunctions  2.2.1
coeffmatrix.hh
Go to the documentation of this file.
1 #ifndef DUNE_COEFFMATRIX_HH
2 #define DUNE_COEFFMATRIX_HH
3 #include <cassert>
4 #include <iostream>
5 #include <vector>
6 #include <dune/common/fvector.hh>
9 
10 namespace Dune
11 {
12  /*************************************************
13  * Default class for storing a coefficient matrix
14  * for the PolynomialBasis. Basically a simple
15  * CRS structure is used. The additional complexity
16  * is due to the storage and efficient evaluation
17  * of higher order derivatives. See the remarks
18  * in tensor.hh which also hold true for this file.
19  *************************************************/
20  template <class Field, class Field2>
21  struct Mult
22  {
23  typedef Field2 BasisEntry;
24  static void add(const Field &vec1, const BasisEntry &vec2,
25  BasisEntry &res)
26  {
27  res += vec1*vec2;
28  }
29  };
30 
31  template <class Field,class Field2, int dimRange>
32  struct Mult< Field,FieldVector<Field2,dimRange> >
33  {
34  typedef FieldVector<Field2,dimRange> BasisEntry;
35  static void add(const Field &vec1, const BasisEntry &vec2,
36  BasisEntry &res)
37  {
38  res.axpy(vec1,vec2);
39  }
40  };
41 
42  template< class F , unsigned int bSize >
44  {
45  public:
46  typedef F Field;
47  static const unsigned int blockSize = bSize;
49 
51  : coeff_(0),
52  rows_(0),
53  skip_(0),
54  numRows_(0),
55  numCols_(0)
56  {}
57 
59  {
60  delete [] coeff_;
61  delete [] rows_;
62  delete [] skip_;
63  }
64 
65  const unsigned int size () const
66  {
67  return numRows_/blockSize;
68  }
69  const unsigned int baseSize () const
70  {
71  return numCols_;
72  }
73 
74  template< class BasisIterator, class FF>
75  void mult ( const BasisIterator &x,
76  unsigned int numLsg,
77  FF *y ) const
78  {
79  typedef typename BasisIterator::Derivatives XDerivatives;
80  assert( numLsg*blockSize <= (size_t)numRows_ );
81  unsigned int row = 0;
82  Field *pos = rows_[ 0 ];
83  unsigned int *skipIt = skip_;
84  XDerivatives val;
85  for( size_t i = 0; i < numLsg; ++i)
86  {
87  for( unsigned int r = 0; r < blockSize; ++r, ++row )
88  {
89  val = 0;
90  BasisIterator itx = x;
91  for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
92  {
93  itx += *skipIt;
94  val.axpy(*pos,*itx);
95  }
96  DerivativeAssign<XDerivatives,FF>::apply(r,val,*(y+i*XDerivatives::size*blockSize));
97  }
98  }
99  }
100  template< class BasisIterator, class Vector>
101  void mult ( const BasisIterator &x,
102  Vector &y ) const
103  {
104  typedef typename Vector::value_type YDerivatives;
105  typedef typename BasisIterator::Derivatives XDerivatives;
106  size_t numLsg = y.size();
107  assert( numLsg*blockSize <= (size_t)numRows_ );
108  unsigned int row = 0;
109  Field *pos = rows_[ 0 ];
110  unsigned int *skipIt = skip_;
111  XDerivatives val;
112  for( size_t i = 0; i < numLsg; ++i)
113  {
114  for( unsigned int r = 0; r < blockSize; ++r, ++row )
115  {
116  val = 0;
117  BasisIterator itx = x;
118  for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
119  {
120  itx += *skipIt;
121  val.axpy(*pos,*itx);
122  }
124  }
125  }
126  }
127  template <unsigned int deriv, class BasisIterator, class Vector>
128  void mult ( const BasisIterator &x,
129  Vector &y ) const
130  {
131  typedef typename Vector::value_type YDerivatives;
132  typedef typename BasisIterator::Derivatives XDerivatives;
133  typedef FieldVector<typename XDerivatives::Field,YDerivatives::dimension> XLFETensor;
134  size_t numLsg = y.size();
135  assert( numLsg*blockSize <= (size_t)numRows_ );
136  unsigned int row = 0;
137  Field *pos = rows_[ 0 ];
138  unsigned int *skipIt = skip_;
139  for( size_t i = 0; i < numLsg; ++i)
140  {
141  XLFETensor val(typename XDerivatives::Field(0));
142  for( unsigned int r = 0; r < blockSize; ++r, ++row )
143  {
144  BasisIterator itx = x;
145  for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
146  {
147  itx += *skipIt;
149  }
150  }
151  field_cast(val,y[i]);
152  }
153  }
154 
155  template< class RowMatrix >
156  void fill ( const RowMatrix &mat, bool verbose=false )
157  {
158  numRows_ = mat.rows();
159  numCols_ = mat.cols();
160  unsigned int size = numRows_*numCols_;
161 
162  delete [] coeff_;
163  delete [] rows_;
164  delete [] skip_;
165 
166  Field* coeff = new Field[ size ];
167  // we always initialize the next skip entry to zero,
168  // including the one following the end, so allocate
169  // size+1 entries so we will stay within the bounds.
170  unsigned int *skip = new unsigned int[ size+1 ];
171  rows_ = new Field*[ numRows_+1 ];
172  std::vector<Field> row( numCols_ );
173 
174  rows_[ 0 ] = coeff;
175  Field *cit = coeff;
176  unsigned int *sit = skip;
177  for( unsigned int r = 0; r < numRows_; ++r )
178  {
179  *sit = 0;
180  mat.row( r, row );
181  for( unsigned int c = 0; c < numCols_; ++c )
182  {
183  const Field &val = row[c];
184  if (val < Zero<Field>() || Zero<Field>() < val)
185  {
186  *cit = val;
187  ++sit;
188  ++cit;
189  *sit = 1;
190  } else
191  {
192  ++(*sit);
193  }
194  }
195  rows_[ r+1 ] = cit;
196  }
197  assert( size_t(rows_[numRows_]-rows_[0]) <= size_t(size) );
198  size = rows_[numRows_]-rows_[0];
199  coeff_ = new Field[ size ];
200  skip_ = new unsigned int[ size ];
201  for (unsigned int i=0;i<size;++i)
202  {
203  coeff_[i] = coeff[i];
204  skip_[i] = skip[i];
205  }
206  for (unsigned int i=0;i<=numRows_;++i)
207  rows_[ i ] = coeff_ + (rows_[ i ] - coeff);
208 
209  delete [] coeff;
210  delete [] skip;
211 
212  if (verbose)
213  std::cout << "Entries: " << (rows_[numRows_]-rows_[0])
214  << " full: " << numCols_*numRows_
215  << std::endl;
216  }
217  // b += a*C[k]
218  template <class Vector>
219  void addRow( unsigned int k, const Field &a, Vector &b) const
220  {
221  assert(k<numRows_);
222  unsigned int j=0;
223  unsigned int *skipIt = skip_ + (rows_[ k ]-rows_[ 0 ]);
224  for( Field *pos = rows_[ k ];
225  pos != rows_[ k+1 ];
226  ++pos, ++skipIt )
227  {
228  j += *skipIt;
229  assert( j < b.size() );
230  b[j] += field_cast<typename Vector::value_type>( (*pos)*a ); // field_cast
231  }
232  }
233  private:
234  SparseCoeffMatrix ( const This &other )
235  : numRows_( other.numRows_ ),
236  numCols_( other.numCols_ )
237  {
238  const unsigned int size = other.rows_[numRows_]-other.rows_[0];
239  coeff_ = new Field[ size ];
240  rows_ = new Field*[ numRows_+1 ];
241  skip_ = new unsigned int[ size ];
242  for (unsigned int i=0;i<size;++i)
243  {
244  coeff_[i] = other.coeff_[i];
245  skip_[i] = other.skip_[i];
246  }
247  for (unsigned int i=0;i<=numRows_;++i)
248  rows_[ i ] = coeff_ + (other.rows_[ i ] - other.coeff_);
249  }
250 
251  This &operator= (const This&);
252  Field *coeff_;
253  Field **rows_;
254  unsigned int *skip_;
255  unsigned int numRows_,numCols_;
256  };
257 
258 }
259 
260 #endif // DUNE_COEFFMATRIX_HH
261