5#ifndef DUNE_COEFFMATRIX_HH
6#define DUNE_COEFFMATRIX_HH
11#include <dune/localfunctions/utility/field.hh>
12#include <dune/localfunctions/utility/tensor.hh>
24 template <
class Field,
class Field2>
27 typedef Field2 BasisEntry;
28 static void add(
const Field &vec1,
const BasisEntry &vec2,
35 template <
class Field,
class Field2,
int dimRange>
36 struct Mult< Field,FieldVector<Field2,dimRange> >
38 typedef FieldVector<Field2,dimRange> BasisEntry;
39 static void add(
const Field &vec1,
const BasisEntry &vec2,
46 template<
class F ,
unsigned int bSize >
47 class SparseCoeffMatrix
51 static const unsigned int blockSize = bSize;
52 typedef SparseCoeffMatrix<Field,blockSize> This;
69 unsigned int size ()
const
71 return numRows_/blockSize;
73 unsigned int baseSize ()
const
78 template<
class BasisIterator,
class FF>
79 void mult (
const BasisIterator &x,
83 typedef typename BasisIterator::Derivatives XDerivatives;
84 assert( numLsg*blockSize <= (
size_t)numRows_ );
86 Field *pos = rows_[ 0 ];
87 unsigned int *skipIt = skip_;
89 for(
size_t i = 0; i < numLsg; ++i)
91 for(
unsigned int r = 0; r < blockSize; ++r, ++row )
94 BasisIterator itx = x;
95 for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
100 DerivativeAssign<XDerivatives,FF>::apply(r,val,*(y+i*XDerivatives::size*blockSize));
104 template<
class BasisIterator,
class Vector>
105 void mult (
const BasisIterator &x,
108 typedef typename Vector::value_type YDerivatives;
109 typedef typename BasisIterator::Derivatives XDerivatives;
110 size_t numLsg = y.size();
111 assert( numLsg*blockSize <= (
size_t)numRows_ );
112 unsigned int row = 0;
113 Field *pos = rows_[ 0 ];
114 unsigned int *skipIt = skip_;
116 for(
size_t i = 0; i < numLsg; ++i)
118 for(
unsigned int r = 0; r < blockSize; ++r, ++row )
121 BasisIterator itx = x;
122 for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
127 DerivativeAssign<XDerivatives,YDerivatives>::apply(r,val,y[i]);
131 template <
unsigned int deriv,
class BasisIterator,
class Vector>
132 void mult (
const BasisIterator &x,
135 typedef typename Vector::value_type YDerivatives;
136 typedef typename BasisIterator::Derivatives XDerivatives;
137 typedef FieldVector<typename XDerivatives::Field,YDerivatives::dimension> XLFETensor;
138 size_t numLsg = y.size();
139 assert( numLsg*blockSize <= (
size_t)numRows_ );
140 unsigned int row = 0;
141 Field *pos = rows_[ 0 ];
142 unsigned int *skipIt = skip_;
143 for(
size_t i = 0; i < numLsg; ++i)
145 XLFETensor val(
typename XDerivatives::Field(0));
146 for(
unsigned int r = 0; r < blockSize; ++r, ++row )
148 BasisIterator itx = x;
149 for( ; pos != rows_[ row+1 ]; ++pos, ++skipIt )
152 LFETensorAxpy<XDerivatives,XLFETensor,deriv>::apply(r,*pos,*itx,val);
159 template<
class RowMatrix >
160 void fill (
const RowMatrix &mat,
bool verbose=
false )
162 numRows_ = mat.rows();
163 numCols_ = mat.cols();
164 unsigned int size = numRows_*numCols_;
170 Field* coeff =
new Field[
size ];
174 unsigned int *skip =
new unsigned int[
size+1 ];
175 rows_ =
new Field*[ numRows_+1 ];
176 std::vector<Field> row( numCols_ );
180 unsigned int *sit = skip;
181 for(
unsigned int r = 0; r < numRows_; ++r )
185 for(
unsigned int c = 0; c < numCols_; ++c )
187 const Field &val = row[c];
188 if (val < Zero<Field>() || Zero<Field>() < val)
201 assert(
size_t(rows_[numRows_]-rows_[0]) <=
size_t(
size) );
202 size = rows_[numRows_]-rows_[0];
203 coeff_ =
new Field[
size ];
204 skip_ =
new unsigned int[
size ];
205 for (
unsigned int i=0; i<
size; ++i)
207 coeff_[i] = coeff[i];
210 for (
unsigned int i=0; i<=numRows_; ++i)
211 rows_[ i ] = coeff_ + (rows_[ i ] - coeff);
217 std::cout <<
"Entries: " << (rows_[numRows_]-rows_[0])
218 <<
" full: " << numCols_*numRows_
222 template <
class Vector>
223 void addRow(
unsigned int k,
const Field &a, Vector &b)
const
227 unsigned int *skipIt = skip_ + (rows_[ k ]-rows_[ 0 ]);
228 for( Field *pos = rows_[ k ];
233 assert( j < b.size() );
234 b[j] += field_cast<typename Vector::value_type>( (*pos)*a );
238 SparseCoeffMatrix (
const This &other )
239 : numRows_( other.numRows_ ),
240 numCols_( other.numCols_ )
242 const unsigned int size = other.rows_[numRows_]-other.rows_[0];
243 coeff_ =
new Field[
size ];
244 rows_ =
new Field*[ numRows_+1 ];
245 skip_ =
new unsigned int[
size ];
246 for (
unsigned int i=0; i<
size; ++i)
248 coeff_[i] = other.coeff_[i];
249 skip_[i] = other.skip_[i];
251 for (
unsigned int i=0; i<=numRows_; ++i)
252 rows_[ i ] = coeff_ + (other.rows_[ i ] - other.coeff_);
255 This &operator= (
const This&);
259 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:13
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:159
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75