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