dune-localfunctions  2.3.0
multiindex.hh
Go to the documentation of this file.
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 
9 #include <dune/common/fvector.hh>
10 
12 
13 namespace 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  {
37 
38  friend std::ostream &operator<<<> ( std::ostream &, const This & );
39 
40  public:
41  static const int dimension = dim;
42 
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 
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  }
113  {
114  remove();
115  vecZ_ = 0;
116  vecOMZ_ = 0;
117  factor_ = 0.;
118  return *this;
119  }
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 
300  typedef Dune::FieldVector< int, dimension > Vector;
301 
302  Vector vecZ_;
303  Vector vecOMZ_;
304  Field factor_;
305 
306  This *next_;
307  };
308 
309  template <int dim, class Field, class 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>
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  {
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  {
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
This operator/(const F &f) const
Definition: multiindex.hh:229
int z(int i) const
Definition: multiindex.hh:89
MultiIndex(int, const This &other)
Definition: multiindex.hh:57
This operator*(const F &f) const
Definition: multiindex.hh:223
MultiIndex(const F &f)
Definition: multiindex.hh:50
This operator+(const This &other) const
Definition: multiindex.hh:246
Field operator/(const Unity< Field > &u, const Field &f)
Definition: field.hh:54
MultiIndex< dim, F > Field
Definition: multiindex.hh:474
int absZ() const
Definition: multiindex.hh:262
Field operator*(const Unity< Field > &u, const Field &f)
Definition: field.hh:48
std::ostream & operator<<(std::ostream &out, const LFEMatrix< Field, aligned > &mat)
Definition: lfematrix.hh:152
This & operator+=(const This &other)
Definition: multiindex.hh:182
This & operator*=(const F &f)
Definition: multiindex.hh:145
~MultiIndex()
Definition: multiindex.hh:84
static const int dimension
Definition: multiindex.hh:41
Definition: multiindex.hh:23
const Field & factor() const
Definition: multiindex.hh:97
int absOMZ() const
Definition: multiindex.hh:270
bool sameMultiIndex(const This &ind)
Definition: multiindex.hh:278
MultiIndex()
Definition: multiindex.hh:43
bool operator==(const This &other) const
Definition: multiindex.hh:138
This & operator/=(const F &f)
Definition: multiindex.hh:153
int omz(int i) const
Definition: multiindex.hh:93
A class representing the zero of a given Field.
Definition: field.hh:76
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
MultiIndex(const This &other)
Definition: multiindex.hh:71
void set(int d, int power=1)
Definition: multiindex.hh:257
This & operator-=(const This &other)
Definition: multiindex.hh:205
Field operator-(const Unity< Field > &u, const Field &f)
Definition: field.hh:42
MultiIndex< dim, F > Field
Definition: multiindex.hh:451
A class representing the unit of a given Field.
Definition: field.hh:27
This & operator=(const This &other)
Definition: multiindex.hh:102
This operator-(const This &other) const
Definition: multiindex.hh:251