Dune Core Modules (2.8.0)

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