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