Dune Core Modules (2.8.0)

multiindex.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_MULTIINDEX_HH
4#define DUNE_MULTIINDEX_HH
5
6#include <vector>
7#include <ostream>
8
10
11#include <dune/localfunctions/utility/field.hh>
12
13namespace Dune
14{
15 /****************************************************************
16 * Provide a Field class which can be used in evaluation methods
17 * to produce MultiIndex presentation of polynomials.
18 ****************************************************************/
19 // Internal Forward Declarations
20 // -----------------------------
21
22 template< int dim, class Field >
23 class MultiIndex;
24
25 template< int dim, class Field >
26 std::ostream &operator<< ( std::ostream &, const MultiIndex< dim,Field > & );
27
28
29
30 // MultiIndex
31 // ----------
32
33 template< int dim,class Field >
34 class MultiIndex
35 {
36 typedef MultiIndex< dim, Field > This;
37
38 friend std::ostream &operator<<<> ( std::ostream &, const This & );
39
40 public:
41 static const int dimension = dim;
42
43 MultiIndex ()
44 : vecZ_( 0 ),
45 vecOMZ_( 0 ),
46 factor_( 1. ),
47 next_( 0 )
48 {}
49 template <class F>
50 explicit MultiIndex (const F &f)
51 : vecZ_( 0 ),
52 vecOMZ_( 0 ),
53 factor_( field_cast<Field>(f) ),
54 next_( 0 )
55 {}
56
57 MultiIndex ( int, const This &other )
58 : vecZ_( other.vecOMZ_ ),
59 vecOMZ_( other.vecZ_ ),
60 factor_( other.factor_ )
61 {
62 assert(!other.next_);
63 if (other.next_)
64 {
65 next_ = new This( *(other.next_) );
66 }
67 else
68 next_ = 0;
69 }
70
71 MultiIndex ( const This &other )
72 : vecZ_( other.vecZ_ ),
73 vecOMZ_( other.vecOMZ_ ),
74 factor_( other.factor_ )
75 {
76 if (other.next_)
77 {
78 next_ = new This( *(other.next_) );
79 }
80 else
81 next_ = 0;
82 }
83
84 ~MultiIndex()
85 {
86 remove();
87 }
88
89 int z(int i) const
90 {
91 return vecZ_[i];
92 }
93 int omz(int i) const
94 {
95 return vecOMZ_[i];
96 }
97 const Field &factor() const
98 {
99 return factor_;
100 }
101
102 This &operator= ( const This &other )
103 {
104 remove();
105 vecZ_ = other.vecZ_;
106 vecOMZ_ = other.vecOMZ_;
107 factor_ = other.factor_;
108 if (other.next_)
109 next_ = new This(*(other.next_));
110 return *this;
111 }
112 This &operator= ( const Zero<This> &f )
113 {
114 remove();
115 vecZ_ = 0;
116 vecOMZ_ = 0;
117 factor_ = 0.;
118 return *this;
119 }
120 This &operator= ( const Unity<This> &f )
121 {
122 remove();
123 vecZ_ = 0;
124 vecOMZ_ = 0;
125 factor_ = 1.;
126 return *this;
127 }
128 template <class F>
129 This &operator= ( const F &f )
130 {
131 remove();
132 vecZ_ = 0;
133 vecOMZ_ = 0;
134 factor_ = field_cast<Field>(f);
135 return *this;
136 }
137
138 bool operator== (const This &other) const
139 {
140 assert(!next_ && !other.next_);
141 return (vecZ_==other.vecZ_ && vecOMZ_==other.vecOMZ_ && factor_==other.factor_);
142 }
143
144 template <class F>
145 This &operator*= ( const F &f )
146 {
147 factor_ *= field_cast<Field>(f);
148 if (next_)
149 (*next_) *= f;
150 return *this;
151 }
152 template <class F>
153 This &operator/= ( const F &f )
154 {
155 factor_ /= field_cast<Field>(f);
156 if (next_)
157 (*next_) /= f;
158 return *this;
159 }
160
161 This &operator*= ( const This &other )
162 {
163 assert(!other.next_);
164 vecZ_ += other.vecZ_;
165 vecOMZ_ += other.vecOMZ_;
166 factor_ *= other.factor_;
167 if (next_)
168 (*next_) *= other;
169 return *this;
170 }
171 This &operator/= ( const This &other )
172 {
173 assert(!other.next_);
174 vecZ_ -= other.vecZ_;
175 vecOMZ_ -= other.vecOMZ_;
176 factor_ /= other.factor_;
177 if (next_)
178 (*next_) /= other;
179 return *this;
180 }
181
182 This &operator+= ( const This &other )
183 {
184 assert(!other.next_);
185 if (std::abs(other.factor_)<1e-10)
186 return *this;
187 if (std::abs(factor_)<1e-10)
188 {
189 *this = other;
190 return *this;
191 }
192 if (!sameMultiIndex(other))
193 {
194 if (next_)
195 (*next_)+=other;
196 else
197 {
198 next_ = new This(other);
199 }
200 }
201 else
202 factor_ += other.factor_;
203 return *this;
204 }
205 This &operator-= ( const This &other )
206 {
207 assert(!other.next_);
208 if (!sameMultiIndex(other))
209 {
210 if (next_)
211 next_-=other;
212 else
213 {
214 next_ = new This(other);
215 }
216 }
217 else
218 factor_ -= other.factor_;
219 return *this;
220 }
221
222 template <class F>
223 This operator* ( const F &f ) const
224 {
225 This z = *this;
226 return (z *= f);
227 }
228 template <class F>
229 This operator/ ( const F &f ) const
230 {
231 This z = *this;
232 return (z /= f);
233 }
234
235 This operator* ( const This &other ) const
236 {
237 This z = *this;
238 return (z *= other);
239 }
240 This operator/ ( const This &other ) const
241 {
242 This z = *this;
243 return (z /= other);
244 }
245
246 This operator+ ( const This &other ) const
247 {
248 This z = *this;
249 return (z += other);
250 }
251 This operator- ( const This &other ) const
252 {
253 This z = *this;
254 return (z -= other);
255 }
256
257 void set ( int d, int power = 1 )
258 {
259 vecZ_[ d ] = power;
260 }
261
262 int absZ () const
263 {
264 int ret = 0;
265 for( int i = 0; i < dimension; ++i )
266 ret += std::abs( vecZ_[ i ] );
267 return ret;
268 }
269
270 int absOMZ() const
271 {
272 int ret = 0;
273 for( int i = 0; i < dimension; ++i )
274 ret += std::abs( vecOMZ_[ i ] );
275 return ret;
276 }
277
278 bool sameMultiIndex(const This &ind)
279 {
280 for( int i = 0; i < dimension; ++i )
281 {
282 if ( vecZ_[i] != ind.vecZ_[i] ||
283 vecOMZ_[i] != vecOMZ_[i] )
284 return false;
285 }
286 return true;
287 }
288
289 private:
290 void remove()
291 {
292 if (next_)
293 {
294 next_->remove();
295 delete next_;
296 next_ = 0;
297 }
298 }
299
301
302 Vector vecZ_;
303 Vector vecOMZ_;
304 Field factor_;
305
306 This *next_;
307 };
308
309 template <int dim, class Field, class F>
310 MultiIndex<dim,Field> operator* ( const F &f,
311 const MultiIndex<dim,Field> &m)
312 {
313 MultiIndex<dim,Field> z = m;
314 return (z *= f);
315 }
316 template <int dim, class Field, class F>
317 MultiIndex<dim,Field> operator/ ( const F &f,
318 const MultiIndex<dim,Field> &m)
319 {
320 MultiIndex<dim,Field> z = m;
321 return (z /= f);
322 }
323
324 template <int d, class F>
325 std::ostream &operator<<(std::ostream& out,const std::vector<MultiIndex<d,F> >& y) {
326 for (unsigned int r=0; r<y.size(); ++r) {
327 out << "f_{" << r << "}(" << char('a');
328 for (int i=1; i<d; ++i)
329 out << "," << char('a'+i);
330 out << ")=";
331 out << y[r] << std::endl;
332 }
333 return out;
334 }
335 template <int d,class F,int dimR>
336 std::ostream &operator<<(std::ostream& out,
337 const std::vector<Dune::FieldVector<MultiIndex<d,F>,dimR> >& y) {
338 out << "\\begin{eqnarray*}" << std::endl;
339 for (unsigned int k=0; k<y.size(); ++k) {
340 out << "f_{" << k << "}(" << char('a');
341 for (int i=1; i<d; ++i)
342 out << "," << char('a'+i);
343 out << ") &=& ( ";
344 out << y[k][0] ;
345 for (unsigned int r=1; r<dimR; ++r) {
346 out << " , " << y[k][r] ;
347 }
348 out << " ) \\\\" << std::endl;
349 }
350 out << "\\end{eqnarray*}" << std::endl;
351 return out;
352 }
353 template <int d,class F,int dimR1,int dimR2>
354 std::ostream &operator<<(std::ostream& out,
355 const std::vector<Dune::FieldMatrix<MultiIndex<d,F>,dimR1,dimR2> >& y) {
356 out << "\\begin{eqnarray*}" << std::endl;
357 for (unsigned int k=0; k<y.size(); ++k) {
358 for (int q=0; q<dimR2; q++) {
359 out << "d_{" << char('a'+q) << "}f_{" << k << "}(" << char('a');
360 for (int i=1; i<d; ++i)
361 out << "," << char('a'+i);
362 out << ") &=& ( ";
363 out << y[k][0][q] ;
364 for (unsigned int r=1; r<dimR1; ++r) {
365 out << " , " << y[k][r][q] ;
366 }
367 out << " ) \\\\" << std::endl;
368 }
369 }
370 out << "\\end{eqnarray*}" << std::endl;
371 return out;
372 }
373 template <int d, class F>
374 std::ostream &operator<<(std::ostream& out,const MultiIndex<d,F>& val)
375 {
376 bool first = true;
377 const MultiIndex<d,F> *m = &val;
378 do {
379 if (m->absZ()==0 && std::abs(m->factor())<1e-10)
380 {
381 if (!m->next_ || !first)
382 {
383 out << "0";
384 break;
385 }
386 else {
387 m = m->next_;
388 continue;
389 }
390 }
391 if (m->factor()>0 && !first)
392 out << " + ";
393 else if (m->factor()<0)
394 if (!first)
395 out << " - ";
396 else
397 out << "- ";
398 else
399 out << " ";
400 first = false;
401 F f = std::abs(m->factor());
402 if (m->absZ()==0)
403 out << f;
404 else {
405 if ( std::abs(f)<1e-10)
406 out << 0;
407 else
408 {
409 F f_1(f);
410 f_1 -= 1.; // better Unity<F>();
411 if ( std::abs(f_1)>1e-10)
412 out << f;
413 int absVal = 0;
414 for (int i=0; i<d; ++i) {
415 if (m->vecZ_[i]==0)
416 continue;
417 else if (m->vecZ_[i]==1)
418 out << char('a'+i);
419 else if (m->vecZ_[i]>0)
420 out << char('a'+i) << "^" << m->vecZ_[i] << "";
421 else if (m->vecZ_[i]<0)
422 out << char('a'+i) << "^" << m->vecZ_[i] << "";
423 absVal += m->vecZ_[i];
424 if (absVal<m->absZ()) out << "";
425 }
426 }
427 }
428 /*
429 if (mi.absOMZ()>0) {
430 for (int i=0;i<=mi.absZ();++i) {
431 if (mi.vecOMZ_[i]==0)
432 continue;
433 else if (mi.vecOMZ_[i]==1)
434 out << (1-char('a'+i));
435 else if (mi.vecOMZ_[i]>0)
436 out << (1-char('a'+i)) << "^" << mi.vecOMZ_[i];
437 else if (mi.vecOMZ_[i]<0)
438 out << (1-char('a'+i)) << "^" << mi.vecOMZ_[i];
439 if (i==mi.absZ()+1) out << "*";
440 }
441 }
442 */
443 m = m->next_;
444 } while (m);
445 return out;
446 }
447
448 template< int dim, class F>
449 struct Unity< MultiIndex< dim, F > >
450 {
451 typedef MultiIndex< dim, F > Field;
452
453 operator Field () const
454 {
455 return Field();
456 }
457
458 Field operator- ( const Field &other ) const
459 {
460 return Field( 1, other );
461 }
462
463 Field operator/ ( const Field &other ) const
464 {
465 return Field() / other;
466 }
467 };
468
469
470
471 template< int dim, class F >
472 struct Zero< MultiIndex< dim,F > >
473 {
474 typedef MultiIndex< dim,F > Field;
475
476 // zero does not acutally exist
477 operator Field ()
478 {
479 return Field(0);
480 }
481 };
482
483 template< int dim, class Field >
484 bool operator< ( const Zero< MultiIndex< dim,Field > > &, const MultiIndex< dim,Field > & )
485 {
486 return true;
487 }
488
489 template< int dim, class Field >
490 bool operator< ( const MultiIndex< dim, Field > &f, const Zero< MultiIndex< dim,Field > > & )
491 {
492 return true;
493 }
494
495}
496
497#endif // #ifndef DUNE_MULTIINDEX_HH
A dense n x m matrix.
Definition: fmatrix.hh:69
Implements a vector constructed from a given type representing a field and a compile-time given size.
EnableIfInterOperable< T1, T2, bool >::type operator<(const RandomAccessIteratorFacade< T1, V1, R1, D > &lhs, const RandomAccessIteratorFacade< T2, V2, R2, D > &rhs)
Comparison operator.
Definition: iteratorfacades.hh:635
EnableIfInterOperable< T1, T2, bool >::type operator==(const ForwardIteratorFacade< T1, V1, R1, D > &lhs, const ForwardIteratorFacade< T2, V2, R2, D > &rhs)
Checks for equality.
Definition: iteratorfacades.hh:235
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
constexpr Mantissa power(Mantissa m, Exponent p)
Power method for integer exponents.
Definition: math.hh:73
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)