dune-localfunctions  2.2.1
basisevaluator.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 
4 #ifndef DUNE_BASISEVALUATOR_HH
5 #define DUNE_BASISEVALUATOR_HH
6 
7 #include <vector>
8 
9 #include <dune/common/fmatrix.hh>
10 #include <dune/common/fvector.hh>
11 #include <dune/common/typetraits.hh>
12 
13 #include <dune/geometry/genericgeometry/topologytypes.hh>
14 
18 
19 namespace Dune
20 {
21  /*******************************************
22  * Should be removed as soon as the Tensor
23  * classes have been revisited. See remarks
24  * in tensor.hh (also hold true here).
25  *******************************************/
26 
27 
28  template <class B>
30  {
31  typedef B Basis;
32  typedef typename Basis::Field Field;
33  typedef typename Basis::DomainVector DomainVector;
34  static const int dimension = Basis::dimension;
35  static const int dimRange = Basis::dimRange;
36 
37  typedef std::vector<Field> Container;
38 
39  template< class Deriv >
40  struct BaseIterator;
41 
42  template <unsigned int deriv>
43  struct Iterator
44  {
47  };
48 
49  unsigned int size() const
50  {
51  return size_;
52  }
53 
54  protected:
55  MonomialEvaluator(const Basis &basis,unsigned int order,unsigned int size)
56  : basis_(basis),
57  order_(order),
58  size_(size),
59  container_(0)
60  {
61  }
62  template <int deriv>
63  void resize()
64  {
66  container_.resize(totalSize);
67  }
69  const Basis &basis_;
70  unsigned int order_,size_;
72  };
73 
74 
75  template< class B >
76  template< class Deriv >
77  struct MonomialEvaluator< B >::BaseIterator
78  {
79  typedef Deriv Derivatives;
80  typedef typename Deriv::Field Field;
81  static const unsigned int blockSize = Deriv::size;
82  typedef Dune::FieldVector<Field,blockSize> Block;
83  static const DerivativeLayout layout = Deriv::layout;
84  static const unsigned int dimDomain = Deriv::dimDomain;
85  static const unsigned int dimRange = Deriv::dimRange;
86 
87  typedef std::vector<Field> Container;
88  typedef typename Container::iterator CIter;
89 
90  explicit BaseIterator ( Container &container )
91  : pos_( container.begin() ),
92  end_( container.end() )
93  {}
94 
95  const Deriv &operator*() const
96  {
97  assert(!done());
98  return reinterpret_cast<const Deriv&>(*pos_);
99  }
100 
101  const Deriv *operator->() const
102  {
103  return &(operator*());
104  }
105 
106  bool done () const
107  {
108  return pos_ == end_;
109  }
110 
111  BaseIterator &operator++ ()
112  {
113  pos_ += blockSize;
114  return *this;
115  }
116 
117  BaseIterator &operator+= ( unsigned int skip )
118  {
119  pos_ += skip*blockSize;
120  return *this;
121  }
122 
123  private:
124  CIter pos_;
125  const CIter end_;
126  };
127 
128  template< class B >
130  : public MonomialEvaluator< B >
131  {
132  typedef B Basis;
133  typedef typename Basis::Field Field;
134  typedef typename Basis::DomainVector DomainVector;
135  typedef std::vector<Field> Container;
136  static const int dimension = Basis::dimension;
137  static const int dimRange = Basis::dimRange;
139 
140  template <unsigned int deriv>
141  struct Iterator : public Base::template Iterator<deriv>
142  {
143  };
144 
145  StandardEvaluator(const Basis &basis)
146  : Base(basis,basis.order(),basis.size())
147  {
148  }
149  template <unsigned int deriv,class DVector>
150  typename Iterator<deriv>::All evaluate(const DVector &x)
151  {
152  Base::template resize<deriv>();
153  basis_.template evaluate<deriv>(x,&(container_[0]));
154  return typename Iterator<deriv>::All(container_);
155  }
157  {
158  Base::template resize<0>();
159  basis_.integrate(&(container_[0]));
160  return typename Iterator<0>::Integrate(container_);
161  }
162 
163  protected:
164  StandardEvaluator ( const Basis &basis, unsigned int size )
165  : Base( basis, basis.order(), size )
166  {}
167 
168  private:
170  using Base::basis_;
171  using Base::container_;
172  };
173 
174 #if 0 // OLD OLD
175  template< class B, class Fill >
176  struct VecEvaluator
177  : public StandardEvaluator< B >
178  {
179  typedef B Basis;
180  typedef typename Basis::Field Field;
181  static const int dimension = Basis::dimension;
182  static const int dimRange = Basis::dimRange*Fill::dimRange;
183  typedef typename Basis::DomainVector DomainVector;
184  typedef std::vector<Field> Container;
185  typedef StandardEvaluator<B> Base;
186 
187  template <unsigned int deriv>
188  struct Iterator
189  {
190  typedef typename Base::template BaseIterator<Derivatives<Field,dimension,dimRange,deriv,Fill::layout> > All;
191  };
192 
193  VecEvaluator ( const Basis &basis, const Fill &fill )
194  : Base( basis, basis.size() ),
195  fill_( fill ),
196  size_( basis.size()*dimRange )
197  {
198  }
199  template <unsigned int deriv>
200  typename Iterator<deriv>::All evaluate(const DomainVector &x)
201  {
202  resize< deriv >();
203  fill_.template apply<deriv>( x,Base::template evaluate<deriv>(x), vecContainer_ );
204  std::vector<Derivatives<Field,dimension,dimRange,deriv,Fill::layout> >& derivContainer =
205  reinterpret_cast<std::vector<Derivatives<Field,dimension,dimRange,deriv,Fill::layout> >&>(vecContainer_);
206  return typename Iterator<deriv>::All(derivContainer);
207  }
208  template <unsigned int deriv,class DVector>
209  typename Iterator<deriv>::All evaluate(const DVector &x)
210  {
211  resize< deriv >();
212  fill_.template apply<deriv>( x,Base::template evaluate<deriv>(x), vecContainer_ );
213  std::vector<Derivatives<Field,dimension,dimRange,deriv,Fill::layout> >& derivContainer =
214  reinterpret_cast<std::vector<Derivatives<Field,dimension,dimRange,deriv,Fill::layout> >&>(vecContainer_);
215  return typename Iterator<deriv>::All(derivContainer);
216  }
217  unsigned int size() const
218  {
219  return size_;
220  }
221 
222  protected:
223  VecEvaluator ( const Basis &basis, const Fill &fill, unsigned int size )
224  : Base( basis, basis.size() ),
225  fill_( fill ),
226  size_( size )
227  {
228  resize< 2 >();
229  }
230 
231  template <int deriv>
232  void resize()
233  {
234  const int totalSize = Derivatives<Field,dimension,dimRange,deriv,derivative>::size*size_;
235  vecContainer_.resize(totalSize);
236  }
237 
238  VecEvaluator(const VecEvaluator&);
239 
240  Container vecContainer_;
241  const Fill &fill_;
242  unsigned int size_;
243  };
244 
245  template <int dimR,DerivativeLayout layout>
246  struct DiagonalFill;
247 
248  template <int dimR>
249  struct DiagonalFill<dimR,derivative>
250  {
251  static const DerivativeLayout layout = derivative;
252  static const int dimRange = dimR;
253  template <int deriv, class Domain, class Iter,class Field>
254  void apply(const Domain &x,
255  Iter iter,std::vector<Field> &vecContainer) const
256  {
257  typedef std::vector<Field> Container;
258  typename Container::iterator vecIter = vecContainer.begin();
259  for ( ; !iter.done(); ++iter)
260  {
261  const typename Iter::Block &block = iter->block();
262  for (int r1=0;r1<dimR;++r1)
263  {
264  unsigned int b = 0;
265  apply<Field>(r1,x,block,b,vecIter);
266  }
267  }
268  }
269  template <class Field, class Domain, class Block,class VecIter>
270  void apply(int r1, const Domain &x,
271  const Block &block,unsigned int &b,
272  VecIter &vecIter) const
273  {
274  unsigned int bStart = b;
275  unsigned int bEnd = b+Block::size;
276  apply<Field>(r1,x,block,bStart,bEnd,vecIter);
277  b=bEnd;
278  }
279  template <class Field, class Domain, class Block,class VecIter>
280  void apply(int r1, const Domain &x,const Block &block,
281  unsigned int bStart, unsigned int bEnd,
282  VecIter &vecIter) const
283  {
284  for (int r2=0;r2<dimR;++r2)
285  {
286  for (unsigned int bb=bStart;bb<bEnd;++bb)
287  {
288  *vecIter = (r1==r2?block[bb]:Field(0));
289  ++vecIter;
290  }
291  }
292  }
293  };
294  template <int dimR>
295  struct DiagonalFill<dimR,value>
296  {
297  static const DerivativeLayout layout = value;
298  static const int dimRange = dimR;
299  template <int deriv, class Domain, class Iter,class Field>
300  void apply(const Domain &x,
301  Iter iter,std::vector<Field> &vecContainer) const
302  {
303  typedef std::vector<Field> Container;
304  typename Container::iterator vecIter = vecContainer.begin();
305  for ( ; !iter.done(); ++iter)
306  {
307  const typename Iter::Block &block = iter->block();
308  for (int r1=0;r1<dimR;++r1)
309  {
310  unsigned int b = 0;
311  apply<Field>(integral_constant<int,deriv>(),r1,x,block,b,vecIter);
312  }
313  }
314  }
315  template <class Field, class Domain, class Block,class VecIter,int deriv>
316  void apply(const integral_constat<int,deriv>&, int r1, const Domain &x,
317  const Block &block,unsigned int &b,
318  VecIter &vecIter) const
319  {
320  apply<Field>(integral_constant<int,deriv-1>(),r1,x,block,b,vecIter);
321  unsigned int bStart = b;
322  unsigned int bEnd = b+LFETensor<Field,Domain::dimension,deriv>::size;
323  apply<Field>(r1,x,block,bStart,bEnd,vecIter);
324  b=bEnd;
325  }
326  template <class Field, class Domain, class Block,class VecIter>
327  void apply(const integral_constant<int,0>&, int r1, const Domain &x,
328  const Block &block,unsigned int &b,
329  VecIter &vecIter) const
330  {
331  apply<Field>(r1,x,block,b,b+1,vecIter);
332  ++b;
333  }
334  template <class Field, class Domain, class Block,class VecIter>
335  void apply(int r1, const Domain &x,const Block &block,
336  unsigned int bStart, unsigned int bEnd,
337  VecIter &vecIter) const
338  {
339  for (int r2=0;r2<dimR;++r2)
340  {
341  for (unsigned int bb=bStart;bb<bEnd;++bb)
342  {
343  *vecIter = (r1==r2?block[bb]:Field(0));
344  ++vecIter;
345  }
346  }
347  }
348  };
349 
350  template <class B,int dimR,DerivativeLayout layout>
351  struct VectorialEvaluator
352  : public VecEvaluator<B,DiagonalFill<dimR,layout> >
353  {
354  typedef DiagonalFill<dimR,layout> Fill;
355  typedef VecEvaluator< B,Fill > Base;
356  VectorialEvaluator(const B &basis)
357  : Base(basis,fill_,basis.size()*dimR)
358  {
359  }
360  private:
361  Fill fill_;
362  };
363 #endif // OLD OLD
364 
365 }
366 
367 #endif