3#ifndef DUNE_COEFFMATRIX_HH
4#define DUNE_COEFFMATRIX_HH
9#include <dune/localfunctions/utility/field.hh>
10#include <dune/localfunctions/utility/tensor.hh>
22 template <
class Field,
class Field2>
25 typedef Field2 BasisEntry;
26 static void add(
const Field &vec1,
const BasisEntry &vec2,
33 template <
class Field,
class Field2,
int dimRange>
34 struct Mult< Field,FieldVector<Field2,dimRange> >
36 typedef FieldVector<Field2,dimRange> BasisEntry;
37 static void add(
const Field &vec1,
const BasisEntry &vec2,
44 template<
class F ,
unsigned int bSize >
45 class SparseCoeffMatrix
49 static const unsigned int blockSize = bSize;
50 typedef SparseCoeffMatrix<Field,blockSize> This;
67 unsigned int size ()
const
69 return numRows_/blockSize;
71 unsigned int baseSize ()
const
76 template<
class BasisIterator,
class FF>
77 void mult (
const BasisIterator &x,
81 typedef typename BasisIterator::Derivatives XDerivatives;
82 assert( numLsg*blockSize <= (
size_t)numRows_ );
84 Field *pos = rows_[ 0 ];
85 unsigned int *skipIt = skip_;
87 for(
size_t i = 0; i < numLsg; ++i)
89 for(
unsigned int r = 0; r < blockSize; ++r, ++row )
92 BasisIterator itx = x;
93 for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
98 DerivativeAssign<XDerivatives,FF>::apply(r,val,*(y+i*XDerivatives::size*blockSize));
102 template<
class BasisIterator,
class Vector>
103 void mult (
const BasisIterator &x,
106 typedef typename Vector::value_type YDerivatives;
107 typedef typename BasisIterator::Derivatives XDerivatives;
108 size_t numLsg = y.size();
109 assert( numLsg*blockSize <= (
size_t)numRows_ );
110 unsigned int row = 0;
111 Field *pos = rows_[ 0 ];
112 unsigned int *skipIt = skip_;
114 for(
size_t i = 0; i < numLsg; ++i)
116 for(
unsigned int r = 0; r < blockSize; ++r, ++row )
119 BasisIterator itx = x;
120 for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
125 DerivativeAssign<XDerivatives,YDerivatives>::apply(r,val,y[i]);
129 template <
unsigned int deriv,
class BasisIterator,
class Vector>
130 void mult (
const BasisIterator &x,
133 typedef typename Vector::value_type YDerivatives;
134 typedef typename BasisIterator::Derivatives XDerivatives;
135 typedef FieldVector<typename XDerivatives::Field,YDerivatives::dimension> XLFETensor;
136 size_t numLsg = y.size();
137 assert( numLsg*blockSize <= (
size_t)numRows_ );
138 unsigned int row = 0;
139 Field *pos = rows_[ 0 ];
140 unsigned int *skipIt = skip_;
141 for(
size_t i = 0; i < numLsg; ++i)
143 XLFETensor val(
typename XDerivatives::Field(0));
144 for(
unsigned int r = 0; r < blockSize; ++r, ++row )
146 BasisIterator itx = x;
147 for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
150 LFETensorAxpy<XDerivatives,XLFETensor,deriv>::apply(r,*pos,*itx,val);
157 template<
class RowMatrix >
158 void fill (
const RowMatrix &mat,
bool verbose=
false )
160 numRows_ = mat.rows();
161 numCols_ = mat.cols();
162 unsigned int size = numRows_*numCols_;
168 Field* coeff =
new Field[ size ];
172 unsigned int *skip =
new unsigned int[ size+1 ];
173 rows_ =
new Field*[ numRows_+1 ];
174 std::vector<Field> row( numCols_ );
178 unsigned int *sit = skip;
179 for(
unsigned int r = 0; r < numRows_; ++r )
183 for(
unsigned int c = 0; c < numCols_; ++c )
185 const Field &val = row[c];
186 if (val < Zero<Field>() || Zero<Field>() < val)
199 assert(
size_t(rows_[numRows_]-rows_[0]) <=
size_t(size) );
200 size = rows_[numRows_]-rows_[0];
201 coeff_ =
new Field[ size ];
202 skip_ =
new unsigned int[ size ];
203 for (
unsigned int i=0; i<size; ++i)
205 coeff_[i] = coeff[i];
208 for (
unsigned int i=0; i<=numRows_; ++i)
209 rows_[ i ] = coeff_ + (rows_[ i ] - coeff);
215 std::cout <<
"Entries: " << (rows_[numRows_]-rows_[0])
216 <<
" full: " << numCols_*numRows_
220 template <
class Vector>
221 void addRow(
unsigned int k,
const Field &a, Vector &b)
const
225 unsigned int *skipIt = skip_ + (rows_[ k ]-rows_[ 0 ]);
226 for( Field *pos = rows_[ k ];
231 assert( j < b.size() );
232 b[j] += field_cast<typename Vector::value_type>( (*pos)*a );
236 SparseCoeffMatrix (
const This &other )
237 : numRows_( other.numRows_ ),
238 numCols_( other.numCols_ )
240 const unsigned int size = other.rows_[numRows_]-other.rows_[0];
241 coeff_ =
new Field[ size ];
242 rows_ =
new Field*[ numRows_+1 ];
243 skip_ =
new unsigned int[ size ];
244 for (
unsigned int i=0; i<size; ++i)
246 coeff_[i] = other.coeff_[i];
247 skip_[i] = other.skip_[i];
249 for (
unsigned int i=0; i<=numRows_; ++i)
250 rows_[ i ] = coeff_ + (other.rows_[ i ] - other.coeff_);
253 This &operator= (
const This&);
257 unsigned int numRows_,numCols_;
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune namespace.
Definition: alignedallocator.hh:11
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157