DUNE-FEM (unstable)

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