DUNE PDELab (2.8)

gridfunctionspaceutilities.hh
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_PDELAB_GRIDFUNCTIONSPACE_GRIDFUNCTIONSPACEUTILITIES_HH
4 #define DUNE_PDELAB_GRIDFUNCTIONSPACE_GRIDFUNCTIONSPACEUTILITIES_HH
5 
6 #include <cstdlib>
7 #include<vector>
8 
9 #include<dune/common/exceptions.hh>
10 #include<dune/common/fvector.hh>
11 #include<dune/common/shared_ptr.hh>
12 
13 #include <dune/localfunctions/common/interfaceswitch.hh>
14 
15 #include"../common/function.hh"
16 #include <dune/pdelab/common/jacobiantocurl.hh>
17 #include"gridfunctionspace.hh"
18 #include <dune/pdelab/gridfunctionspace/localfunctionspace.hh>
19 #include <dune/pdelab/gridfunctionspace/lfsindexcache.hh>
20 
21 namespace Dune {
22  namespace PDELab {
23 
27 
28  //===============================================================
29  // output: convert grid function space to discrete grid function
30  //===============================================================
31 
32 
53  template<typename T, typename X>
54  class DiscreteGridFunction
55  : public GridFunctionBase<
56  GridFunctionTraits<
57  typename T::Traits::GridViewType,
58  typename BasisInterfaceSwitch<
59  typename FiniteElementInterfaceSwitch<
60  typename T::Traits::FiniteElementType
61  >::Basis
62  >::RangeField,
63  BasisInterfaceSwitch<
64  typename FiniteElementInterfaceSwitch<
65  typename T::Traits::FiniteElementType
66  >::Basis
67  >::dimRange,
68  typename BasisInterfaceSwitch<
69  typename FiniteElementInterfaceSwitch<
70  typename T::Traits::FiniteElementType
71  >::Basis
72  >::Range
73  >,
74  DiscreteGridFunction<T,X>
75  >
76  {
77  typedef T GFS;
78 
79  typedef typename Dune::BasisInterfaceSwitch<
80  typename FiniteElementInterfaceSwitch<
81  typename T::Traits::FiniteElementType
82  >::Basis
83  > BasisSwitch;
84  typedef GridFunctionBase<
85  GridFunctionTraits<
86  typename T::Traits::GridViewType,
87  typename BasisSwitch::RangeField,
88  BasisSwitch::dimRange,
89  typename BasisSwitch::Range
90  >,
91  DiscreteGridFunction<T,X>
92  > BaseT;
93 
94  public:
95  typedef typename BaseT::Traits Traits;
96 
102  DiscreteGridFunction (const GFS& gfs, const X& x_)
103  : pgfs(stackobject_to_shared_ptr(gfs))
104  , lfs(gfs)
105  , lfs_cache(lfs)
106  , x_view(x_)
107  , xl(gfs.maxLocalSize())
108  , yb(gfs.maxLocalSize())
109  {
110  }
111 
117  DiscreteGridFunction (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_)
118  : pgfs(gfs)
119  , lfs(*gfs)
120  , lfs_cache(lfs)
121  , x_view(*x_)
122  , xl(gfs->maxLocalSize())
123  , yb(gfs->maxLocalSize())
124  , px(x_) // FIXME: The LocalView should handle a shared_ptr correctly!
125  {
126  }
127 
128  // Evaluate
129  inline void evaluate (const typename Traits::ElementType& e,
130  const typename Traits::DomainType& x,
131  typename Traits::RangeType& y) const
132  {
133  typedef FiniteElementInterfaceSwitch<
134  typename Dune::PDELab::LocalFunctionSpace<GFS>::Traits::FiniteElementType
135  > FESwitch;
136  lfs.bind(e);
137  lfs_cache.update();
138  x_view.bind(lfs_cache);
139  x_view.read(xl);
140  x_view.unbind();
141  FESwitch::basis(lfs.finiteElement()).evaluateFunction(x,yb);
142  y = 0;
143  for (unsigned int i=0; i<yb.size(); i++)
144  {
145  y.axpy(xl[i],yb[i]);
146  }
147  }
148 
150  inline const typename Traits::GridViewType& getGridView () const
151  {
152  return pgfs->gridView();
153  }
154 
155  private:
156 
157  typedef LocalFunctionSpace<GFS> LFS;
158  typedef LFSIndexCache<LFS> LFSCache;
159  typedef typename X::template ConstLocalView<LFSCache> XView;
160 
161  std::shared_ptr<GFS const> pgfs;
162  mutable LFS lfs;
163  mutable LFSCache lfs_cache;
164  mutable XView x_view;
165  mutable std::vector<typename Traits::RangeFieldType> xl;
166  mutable std::vector<typename Traits::RangeType> yb;
167  std::shared_ptr<const X> px; // FIXME: dummy pointer to make sure we take ownership of X
168  };
169 
179  template<typename T, typename X>
180  class DiscreteGridFunctionCurl :
181  public GridFunctionBase<
182  GridFunctionTraits<
183  typename T::Traits::GridViewType,
184  typename JacobianToCurl<typename T::Traits::FiniteElementType::
185  Traits::LocalBasisType::Traits::JacobianType>::CurlField,
186  JacobianToCurl<typename T::Traits::FiniteElementType::Traits::LocalBasisType::
187  Traits::JacobianType>::dimCurl,
188  typename JacobianToCurl<typename T::Traits::FiniteElementType::
189  Traits::LocalBasisType::Traits::JacobianType>::Curl
190  >,
191  DiscreteGridFunctionCurl<T,X>
192  >
193  {
194  typedef T GFS;
195  typedef typename T::Traits::FiniteElementType::Traits::LocalBasisType::Traits::
196  JacobianType Jacobian;
197  typedef JacobianToCurl<Jacobian> J2C;
198 
199  public:
200  typedef GridFunctionTraits<
201  typename T::Traits::GridViewType,
202  typename J2C::CurlField, J2C::dimCurl, typename J2C::Curl
203  > Traits;
204 
210  DiscreteGridFunctionCurl(const GFS& gfs, const X& x_)
211  : pgfs(stackobject_to_shared_ptr(gfs))
212  , lfs(gfs)
213  , lfs_cache(lfs)
214  , x_view(x_)
215  , xl(gfs.maxLocalSize())
216  , jacobian(gfs.maxLocalSize())
217  , yb(gfs.maxLocalSize())
218  , px(stackobject_to_shared_ptr(x_))
219  {}
220 
226  DiscreteGridFunctionCurl(std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_)
227  : pgfs(gfs)
228  , lfs(*gfs)
229  , lfs_cache(lfs)
230  , x_view(*x_)
231  , xl(gfs->maxLocalSize())
232  , jacobian(gfs->maxLocalSize())
233  , yb(gfs->maxLocalSize())
234  , px(x_)
235  {}
236 
237  // Evaluate
238  void evaluate (const typename Traits::ElementType& e,
239  const typename Traits::DomainType& x,
240  typename Traits::RangeType& y) const
241  {
242  static const J2C& j2C = J2C();
243 
244  lfs.bind();
245  lfs_cache.update();
246  x_view.bind(lfs_cache);
247  x_view.read(xl);
248  x_view.unbind();
249 
250  lfs.finiteElement().basis().evaluateJacobian(x,jacobian);
251 
252  y = 0;
253  for (std::size_t i=0; i < lfs.size(); i++) {
254  j2C(jacobian[i], yb);
255  y.axpy(xl[i], yb);
256  }
257  }
258 
260  const typename Traits::GridViewType& getGridView() const
261  { return pgfs->gridView(); }
262 
263  private:
264  typedef GridFunctionBase<Traits, DiscreteGridFunctionCurl<T,X> >
265  BaseT;
266  typedef LocalFunctionSpace<GFS> LFS;
267  typedef LFSIndexCache<LFS> LFSCache;
268  typedef typename X::template ConstLocalView<LFSCache> XView;
269 
270  std::shared_ptr<GFS const> pgfs;
271  mutable LFS lfs;
272  mutable LFSCache lfs_cache;
273  mutable XView x_view;
274  mutable std::vector<typename Traits::RangeFieldType> xl;
275  mutable std::vector<Jacobian> jacobian;
276  mutable std::vector<typename Traits::RangeType> yb;
277  std::shared_ptr<const X> px; // FIXME: dummy pointer to make sure we take ownership of X
278  };
279 
281 
293  template<typename GV, typename RangeFieldType, int dimRangeOfBasis>
294  struct DiscreteGridFunctionCurlTraits {
295  static_assert(AlwaysFalse<GV>::value,
296  "DiscreteGridFunctionCurl (and friends) work in 2D "
297  "and 3D only");
298  };
300 
306  template<typename GV, typename RangeFieldType>
307  struct DiscreteGridFunctionCurlTraits<GV, RangeFieldType, 1>
308  : public GridFunctionTraits<GV,
309  RangeFieldType, 2,
310  FieldVector<RangeFieldType, 2> >
311  {
312  static_assert(GV::dimensionworld == 2,
313  "World dimension of grid must be 2 for the curl of a "
314  "scalar (1D) quantity");
315  };
317 
322  template<typename GV, typename RangeFieldType>
323  struct DiscreteGridFunctionCurlTraits<GV, RangeFieldType, 2>
324  : public GridFunctionTraits<GV,
325  RangeFieldType, 1,
326  FieldVector<RangeFieldType, 1> >
327  {
328  static_assert(GV::dimensionworld == 2,
329  "World dimension of grid must be 2 for the curl of a"
330  "2D quantity");
331  };
333 
338  template<typename GV, typename RangeFieldType>
339  struct DiscreteGridFunctionCurlTraits<GV, RangeFieldType, 3>
340  : public GridFunctionTraits<GV,
341  RangeFieldType, 3,
342  FieldVector<RangeFieldType, 3> >
343  {
344  static_assert(GV::dimensionworld == 3,
345  "World dimension of grid must be 3 for the curl of a"
346  "3D quantity");
347  };
348 
351 
365  template<typename T, typename X>
366  class DiscreteGridFunctionGlobalCurl
367  : public GridFunctionBase<
368  DiscreteGridFunctionCurlTraits<
369  typename T::Traits::GridViewType,
370  typename T::Traits::FiniteElementType::Traits::
371  LocalBasisType::Traits::RangeFieldType,
372  T::Traits::FiniteElementType::Traits::LocalBasisType::Traits::
373  dimRange>,
374  DiscreteGridFunctionGlobalCurl<T,X> >
375  {
376  public:
377  typedef DiscreteGridFunctionCurlTraits<
378  typename T::Traits::GridViewType,
379  typename T::Traits::FiniteElementType::Traits::
380  LocalBasisType::Traits::RangeFieldType,
381  T::Traits::FiniteElementType::Traits::LocalBasisType::Traits::
382  dimRange> Traits;
383 
384  private:
385  typedef T GFS;
386  typedef GridFunctionBase<
387  Traits,
388  DiscreteGridFunctionGlobalCurl<T,X> > BaseT;
389  typedef typename T::Traits::FiniteElementType::Traits::
390  LocalBasisType::Traits LBTraits;
391 
392  public:
398  DiscreteGridFunctionGlobalCurl (const GFS& gfs, const X& x_)
399  : pgfs(stackobject_to_shared_ptr(gfs))
400  , lfs(gfs)
401  , lfs_cache(lfs)
402  , x_view(x_)
403  , xl(gfs.maxLocalSize())
404  , J(gfs.maxLocalSize())
405  , px(stackobject_to_shared_ptr(x_))
406  {}
407 
413  DiscreteGridFunctionGlobalCurl(std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_)
414  : pgfs(gfs)
415  , lfs(*gfs)
416  , lfs_cache(lfs)
417  , x_view(*x_)
418  , xl(gfs->maxLocalSize())
419  , J(gfs->maxLocalSize())
420  , px(x_)
421  {}
422 
423  // Evaluate
424  inline void evaluate (const typename Traits::ElementType& e,
425  const typename Traits::DomainType& x,
426  typename Traits::RangeType& y) const
427  {
428  lfs.bind(e);
429  lfs_cache.update();
430  x_view.bind(lfs_cache);
431  x_view.read(xl);
432  x_view.unbind();
433 
434  lfs.finiteElement().localBasis().
435  evaluateJacobianGlobal(x,J,e.geometry());
436  y = 0;
437  for (unsigned int i=0; i<J.size(); i++)
438  // avoid a "case label value exceeds maximum value for type"
439  // warning: since dimRange is an anonymous enum, its type may
440  // contain only the values 0 and 1, resulting in a warning.
441  switch(unsigned(Traits::dimRange)) {
442  case 1:
443  y[0] += xl[i] * J[i][0][1];
444  y[1] += xl[i] * -J[i][0][0];
445  break;
446  case 2:
447  y[0] += xl[i]*(J[i][1][0] - J[i][0][1]);
448  break;
449  case 3:
450  y[0] += xl[i]*(J[i][2][1] - J[i][1][2]);
451  y[1] += xl[i]*(J[i][0][2] - J[i][2][0]);
452  y[2] += xl[i]*(J[i][1][0] - J[i][0][1]);
453  break;
454  default:
455  //how did that pass all the static asserts?
456  std::abort();
457  }
458  }
459 
461  inline const typename Traits::GridViewType& getGridView () const
462  {
463  return pgfs->gridView();
464  }
465 
466  private:
467  typedef LocalFunctionSpace<GFS> LFS;
468  typedef LFSIndexCache<LFS> LFSCache;
469  typedef typename X::template ConstLocalView<LFSCache> XView;
470 
471  std::shared_ptr<GFS const> pgfs;
472  mutable LFS lfs;
473  mutable LFSCache lfs_cache;
474  mutable XView x_view;
475  mutable std::vector<typename Traits::RangeFieldType> xl;
476  mutable std::vector<typename T::Traits::FiniteElementType::Traits::LocalBasisType::Traits::JacobianType> J;
477  std::shared_ptr<const X> px; // FIXME: dummy pointer to make sure we take ownership of X
478  };
479 
482 
490  template<typename T, typename X>
491  class DiscreteGridFunctionGradient
492  : public GridFunctionBase<
493  GridFunctionTraits<
494  typename T::Traits::GridViewType,
495  typename T::Traits::FiniteElementType::Traits::LocalBasisType
496  ::Traits::RangeFieldType,
497  T::Traits::FiniteElementType::Traits::LocalBasisType::Traits
498  ::dimDomain,
499  FieldVector<
500  typename T::Traits::FiniteElementType::Traits
501  ::LocalBasisType::Traits::RangeFieldType,
502  T::Traits::FiniteElementType::Traits::LocalBasisType::Traits
503  ::dimDomain> >,
504  DiscreteGridFunctionGradient<T,X> >
505  {
506  typedef T GFS;
507  typedef typename GFS::Traits::FiniteElementType::Traits::
508  LocalBasisType::Traits LBTraits;
509 
510  public:
511  typedef GridFunctionTraits<
512  typename GFS::Traits::GridViewType,
513  typename LBTraits::RangeFieldType,
514  LBTraits::dimDomain,
515  FieldVector<
516  typename LBTraits::RangeFieldType,
517  LBTraits::dimDomain> > Traits;
518 
519  private:
520  typedef GridFunctionBase<
521  Traits,
522  DiscreteGridFunctionGradient<T,X> > BaseT;
523 
524  public:
530  DiscreteGridFunctionGradient (const GFS& gfs, const X& x_)
531  : pgfs(stackobject_to_shared_ptr(gfs))
532  , lfs(gfs)
533  , lfs_cache(lfs)
534  , x_view(x_)
535  , xl(lfs.size())
536  { }
537 
543  DiscreteGridFunctionGradient (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_)
544  : pgfs(gfs)
545  , lfs(*gfs)
546  , lfs_cache(lfs)
547  , x_view(*x_)
548  , xl(lfs.size())
549  { }
550 
551  // Evaluate
552  inline void evaluate (const typename Traits::ElementType& e,
553  const typename Traits::DomainType& x,
554  typename Traits::RangeType& y) const
555  {
556  // get and bind local functions space
557  lfs.bind(e);
558  lfs_cache.update();
559  x_view.bind(lfs_cache);
560 
561  // get local coefficients
562  xl.resize(lfs.size());
563  x_view.read(xl);
564  x_view.unbind();
565 
566  // get Jacobian of geometry
567  const typename Traits::ElementType::Geometry::JacobianInverseTransposed
568  JgeoIT = e.geometry().jacobianInverseTransposed(x);
569 
570  // get local Jacobians/gradients of the shape functions
571  std::vector<typename LBTraits::JacobianType> J(lfs.size());
572  lfs.finiteElement().localBasis().evaluateJacobian(x,J);
573 
574  typename Traits::RangeType gradphi;
575  y = 0;
576  for(unsigned int i = 0; i < lfs.size(); ++i) {
577  // compute global gradient of shape function i
578  gradphi = 0;
579  JgeoIT.umv(J[i][0], gradphi);
580 
581  // sum up global gradients, weighting them with the appropriate coeff
582  y.axpy(xl[i], gradphi);
583  }
584 
585  }
586 
588  inline const typename Traits::GridViewType& getGridView () const
589  {
590  return pgfs->gridView();
591  }
592 
593  private:
594  typedef LocalFunctionSpace<GFS> LFS;
595  typedef LFSIndexCache<LFS> LFSCache;
596  typedef typename X::template ConstLocalView<LFSCache> XView;
597 
598  std::shared_ptr<GFS const> pgfs;
599  mutable LFS lfs;
600  mutable LFSCache lfs_cache;
601  mutable XView x_view;
602  mutable std::vector<typename Traits::RangeFieldType> xl;
603  };
604 
609  template<typename T, typename X>
610  class DiscreteGridFunctionPiola
611  : public GridFunctionBase<
612  GridFunctionTraits<
613  typename T::Traits::GridViewType,
614  typename T::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
615  T::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimRange,
616  typename T::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType
617  >,
618  DiscreteGridFunctionPiola<T,X>
619  >
620  {
621  typedef T GFS;
622 
623  typedef GridFunctionBase<
624  GridFunctionTraits<
625  typename T::Traits::GridViewType,
626  typename T::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
627  T::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimRange,
628  typename T::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType
629  >,
630  DiscreteGridFunctionPiola<T,X>
631  > BaseT;
632 
633  public:
634  typedef typename BaseT::Traits Traits;
635 
640  DiscreteGridFunctionPiola (const GFS& gfs, const X& x_)
641  : pgfs(stackobject_to_shared_ptr(gfs))
642  , lfs(gfs)
643  , lfs_cache(lfs)
644  , x_view(x_)
645  , xl(pgfs->maxLocalSize())
646  , yb(pgfs->maxLocalSize())
647  {
648  }
649 
655  DiscreteGridFunctionPiola (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_)
656  : pgfs(gfs)
657  , lfs(*gfs)
658  , lfs_cache(lfs)
659  , x_view(*x_)
660  , xl(pgfs->maxLocalSize())
661  , yb(pgfs->maxLocalSize())
662  {
663  }
664 
665  inline void evaluate (const typename Traits::ElementType& e,
666  const typename Traits::DomainType& x,
667  typename Traits::RangeType& y) const
668  {
669  // evaluate shape function on the reference element as before
670  lfs.bind(e);
671  lfs_cache.update();
672  x_view.bind(lfs_cache);
673  x_view.read(xl);
674  x_view.unbind();
675 
676  lfs.finiteElement().localBasis().evaluateFunction(x,yb);
677  typename Traits::RangeType yhat;
678  yhat = 0;
679  for (unsigned int i=0; i<yb.size(); i++)
680  yhat.axpy(xl[i],yb[i]);
681 
682  // apply Piola transformation
683  typename Traits::ElementType::Geometry::JacobianInverseTransposed
684  J = e.geometry().jacobianInverseTransposed(x);
685  J.invert();
686  y = 0;
687  J.umtv(yhat,y);
688  y /= J.determinant();
689  }
690 
692  inline const typename Traits::GridViewType& getGridView () const
693  {
694  return pgfs->gridView();
695  }
696 
697  private:
698 
699  typedef LocalFunctionSpace<GFS> LFS;
700  typedef LFSIndexCache<LFS> LFSCache;
701  typedef typename X::template ConstLocalView<LFSCache> XView;
702 
703  std::shared_ptr<GFS const> pgfs;
704  mutable LFS lfs;
705  mutable LFSCache lfs_cache;
706  mutable XView x_view;
707  mutable std::vector<typename Traits::RangeFieldType> xl;
708  mutable std::vector<typename Traits::RangeType> yb;
709 
710  };
711 
723 #ifdef __clang__
724  // clang is too stupid to correctly apply the constexpr qualifier of staticDegree in this context
725  template<typename T, typename X, std::size_t dimR = TypeTree::StaticDegree<T>::value>
726 #else
727  template<typename T, typename X, std::size_t dimR = TypeTree::StaticDegree<T>::value>
728 #endif
729  class VectorDiscreteGridFunction
730  : public GridFunctionBase<
731  GridFunctionTraits<
732  typename T::Traits::GridViewType,
733  typename T::template Child<0>::Type::Traits::FiniteElementType
734  ::Traits::LocalBasisType::Traits::RangeFieldType,
735  dimR,
736  Dune::FieldVector<
737  typename T::template Child<0>::Type::Traits::FiniteElementType
738  ::Traits::LocalBasisType::Traits::RangeFieldType,
739  dimR
740  >
741  >,
742  VectorDiscreteGridFunction<T,X>
743  >
744  {
745  typedef T GFS;
746 
747  typedef GridFunctionBase<
748  GridFunctionTraits<
749  typename T::Traits::GridViewType,
750  typename T::template Child<0>::Type::Traits::FiniteElementType
751  ::Traits::LocalBasisType::Traits::RangeFieldType,
752  dimR,
753  Dune::FieldVector<
754  typename T::template Child<0>::Type::Traits::FiniteElementType
755  ::Traits::LocalBasisType::Traits::RangeFieldType,
756  dimR
757  >
758  >,
759  VectorDiscreteGridFunction<T,X,dimR>
760  > BaseT;
761 
762  public:
763  typedef typename BaseT::Traits Traits;
764  typedef typename T::template Child<0>::Type ChildType;
765  typedef typename ChildType::Traits::FiniteElementType
766  ::Traits::LocalBasisType::Traits::RangeFieldType RF;
767  typedef typename ChildType::Traits::FiniteElementType
768  ::Traits::LocalBasisType::Traits::RangeType RT;
769 
771 
776  VectorDiscreteGridFunction(const GFS& gfs, const X& x_,
777  std::size_t start = 0)
778  : pgfs(stackobject_to_shared_ptr(gfs))
779  , lfs(gfs)
780  , lfs_cache(lfs)
781  , x_view(x_)
782  , xl(gfs.maxLocalSize())
783  , yb(gfs.maxLocalSize())
784  {
785  for(std::size_t i = 0; i < dimR; ++i)
786  remap[i] = i + start;
787  }
788 
795  VectorDiscreteGridFunction (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_,
796  std::size_t start = 0)
797  : pgfs(gfs)
798  , lfs(*gfs)
799  , lfs_cache(lfs)
800  , x_view(*x_)
801  , xl(pgfs->maxLocalSize())
802  , yb(pgfs->maxLocalSize())
803  {
804  for(std::size_t i = 0; i < dimR; ++i)
805  remap[i] = i + start;
806  }
807 
809 
820  template<class Remap>
821  VectorDiscreteGridFunction(const GFS& gfs, const X& x_,
822  const Remap &remap_)
823  : pgfs(stackobject_to_shared_ptr(gfs))
824  , lfs(gfs)
825  , lfs_cache(lfs)
826  , x_view(x_)
827  , xl(gfs.maxLocalSize())
828  , yb(gfs.maxLocalSize())
829  , px(stackobject_to_shared_ptr(x_))
830  {
831  for(std::size_t i = 0; i < dimR; ++i)
832  remap[i] = remap_[i];
833  }
834 
847  template<class Remap>
848  VectorDiscreteGridFunction (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_,
849  const Remap &remap_)
850  : pgfs(gfs)
851  , lfs(*gfs)
852  , lfs_cache(lfs)
853  , x_view(*x_)
854  , xl(pgfs->maxLocalSize())
855  , yb(pgfs->maxLocalSize())
856  , px(x_)
857  {
858  for(std::size_t i = 0; i < dimR; ++i)
859  remap[i] = remap_[i];
860  }
861 
862  inline void evaluate (const typename Traits::ElementType& e,
863  const typename Traits::DomainType& x,
864  typename Traits::RangeType& y) const
865  {
866  lfs.bind(e);
867  lfs_cache.update();
868  x_view.bind(lfs_cache);
869  x_view.read(xl);
870  x_view.unbind();
871  for (unsigned int k=0; k < dimR; k++)
872  {
873  lfs.child(remap[k]).finiteElement().localBasis().
874  evaluateFunction(x,yb);
875  y[k] = 0.0;
876  for (unsigned int i=0; i<yb.size(); i++)
877  y[k] += xl[lfs.child(remap[k]).localIndex(i)]*yb[i];
878  }
879  }
880 
882  inline const typename Traits::GridViewType& getGridView () const
883  {
884  return pgfs->gridView();
885  }
886 
887  private:
888  typedef LocalFunctionSpace<GFS> LFS;
889  typedef LFSIndexCache<LFS> LFSCache;
890  typedef typename X::template ConstLocalView<LFSCache> XView;
891 
892  std::shared_ptr<GFS const> pgfs;
893  std::size_t remap[dimR];
894  mutable LFS lfs;
895  mutable LFSCache lfs_cache;
896  mutable XView x_view;
897  mutable std::vector<RF> xl;
898  mutable std::vector<RT> yb;
899  std::shared_ptr<const X> px; // FIXME: dummy pointer to make sure we take ownership of X
900  };
901 
907  template<typename T, typename X>
908  class VectorDiscreteGridFunctionGradient
909  : public GridFunctionBase<
910  GridFunctionTraits<
911  typename T::Traits::GridViewType,
912  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
913  //T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimDomain,
914  TypeTree::StaticDegree<T>::value,
915  Dune::FieldMatrix<
916  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
917  TypeTree::StaticDegree<T>::value,
918  T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimDomain
919  >
920  >,
921  VectorDiscreteGridFunctionGradient<T,X>
922  >
923  {
924  typedef T GFS;
925 
926  typedef GridFunctionBase<
927  GridFunctionTraits<
928  typename T::Traits::GridViewType,
929  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
930  //T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimDomain,
931  TypeTree::StaticDegree<T>::value,
932  Dune::FieldMatrix<
933  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
934  TypeTree::StaticDegree<T>::value,
935  T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimDomain>
936  >,
937  VectorDiscreteGridFunctionGradient<T,X>
938  > BaseT;
939 
940  public:
941  typedef typename BaseT::Traits Traits;
942  typedef typename T::template Child<0>::Type ChildType;
943  typedef typename ChildType::Traits::FiniteElementType::Traits::LocalBasisType::Traits LBTraits;
944 
945  typedef typename LBTraits::RangeFieldType RF;
946  typedef typename LBTraits::JacobianType JT;
947 
953  VectorDiscreteGridFunctionGradient (const GFS& gfs, const X& x_)
954  : pgfs(stackobject_to_shared_ptr(gfs))
955  , lfs(gfs)
956  , lfs_cache(lfs)
957  , x_view(x_)
958  , xl(gfs.maxLocalSize())
959  , J(gfs.maxLocalSize())
960  {
961  }
962 
968  VectorDiscreteGridFunctionGradient (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_)
969  : pgfs(gfs)
970  , lfs(*gfs)
971  , lfs_cache(lfs)
972  , x_view(*x_)
973  , xl(pgfs->maxLocalSize())
974  , J(pgfs->maxLocalSize())
975  {
976  }
977 
978  inline void evaluate(const typename Traits::ElementType& e,
979  const typename Traits::DomainType& x,
980  typename Traits::RangeType& y) const
981  {
982  // get and bind local functions space
983  lfs.bind(e);
984  lfs_cache.update();
985  x_view.bind(lfs_cache);
986  x_view.read(xl);
987  x_view.unbind();
988 
989  // get Jacobian of geometry
990  const typename Traits::ElementType::Geometry::JacobianInverseTransposed
991  JgeoIT = e.geometry().jacobianInverseTransposed(x);
992 
993  y = 0.0;
994 
995  // Loop over PowerLFS and calculate gradient for each child separately
996  for(unsigned int k = 0; k != TypeTree::degree(lfs); ++k)
997  {
998  // get local Jacobians/gradients of the shape functions
999  std::vector<typename LBTraits::JacobianType> J(lfs.child(k).size());
1000  lfs.child(k).finiteElement().localBasis().evaluateJacobian(x,J);
1001 
1002  Dune::FieldVector<RF,LBTraits::dimDomain> gradphi;
1003  for (typename LFS::Traits::SizeType i=0; i<lfs.child(k).size(); i++)
1004  {
1005  gradphi = 0;
1006  JgeoIT.umv(J[i][0], gradphi);
1007 
1008  y[k].axpy(xl[lfs.child(k).localIndex(i)], gradphi);
1009  }
1010  }
1011  }
1012 
1013 
1015  inline const typename Traits::GridViewType& getGridView () const
1016  {
1017  return pgfs->gridView();
1018  }
1019 
1020  private:
1021  typedef LocalFunctionSpace<GFS> LFS;
1022  typedef LFSIndexCache<LFS> LFSCache;
1023  typedef typename X::template ConstLocalView<LFSCache> XView;
1024 
1025  std::shared_ptr<GFS const> pgfs;
1026  mutable LFS lfs;
1027  mutable LFSCache lfs_cache;
1028  mutable XView x_view;
1029  mutable std::vector<RF> xl;
1030  mutable std::vector<JT> J;
1031  std::shared_ptr<const X> px; // FIXME: dummy pointer to make sure we take ownership of X
1032  };
1033 
1047  template<typename Mat, typename RF, std::size_t size>
1048  struct SingleDerivativeComputationHelper {
1053  template<typename T>
1054  static inline RF compute_derivative(const Mat& mat, const T& t, const unsigned int k)
1055  {
1056  // setting this to zero is just a test if the template specialization work
1057  Dune::FieldVector<RF,size> grad_phi(0.0);
1058  mat.umv(t,grad_phi);
1059  return grad_phi[k];
1060  // return 0.0;
1061  }
1062  };
1063 
1072  template<typename RF, std::size_t size>
1073  struct SingleDerivativeComputationHelper<Dune::FieldMatrix<RF,size,size>,RF,size> {
1079  template<typename T>
1080  static inline RF compute_derivative(const Dune::FieldMatrix<RF,size,size>& mat, const T& t, const unsigned int k)
1081  {
1082  return mat[k]*t;
1083  }
1084  };
1085 
1096  template<typename RF, std::size_t size>
1097  struct SingleDerivativeComputationHelper<Dune::DiagonalMatrix<RF,size>,RF,size> {
1103  template<typename T>
1104  static inline RF compute_derivative(const Dune::DiagonalMatrix<RF,size>& mat, const T& t, const unsigned int k)
1105  {
1106  return mat[k][k]*t[k];
1107  }
1108  };
1109 
1120  template<typename T, typename X>
1121  class VectorDiscreteGridFunctionDiv
1122  : public Dune::PDELab::GridFunctionBase<
1123  Dune::PDELab::GridFunctionTraits<
1124  typename T::Traits::GridViewType,
1125  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
1126  T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimRange,
1127  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType>,
1128  VectorDiscreteGridFunctionDiv<T,X> >
1129  {
1130  typedef T GFS;
1131 
1132  typedef Dune::PDELab::GridFunctionBase<
1133  Dune::PDELab::GridFunctionTraits<
1134  typename T::Traits::GridViewType,
1135  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
1136  T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimRange,
1137  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType>,
1138  VectorDiscreteGridFunctionDiv<T,X> > BaseT;
1139  public :
1140  typedef typename BaseT::Traits Traits;
1141  typedef typename T::template Child<0>::Type ChildType;
1142  typedef typename ChildType::Traits::FiniteElementType::Traits::LocalBasisType::Traits LBTraits;
1143 
1144  typedef typename LBTraits::RangeFieldType RF;
1145  typedef typename LBTraits::JacobianType JT;
1146 
1152  VectorDiscreteGridFunctionDiv(const GFS& gfs, const X& x_)
1153  : pgfs(stackobject_to_shared_ptr(gfs))
1154  , lfs(gfs)
1155  , lfs_cache(lfs)
1156  , x_view(x_)
1157  , xl(gfs.maxLocalSize())
1158  , J(gfs.maxLocalSize())
1159  {
1160  static_assert(LBTraits::dimDomain == TypeTree::StaticDegree<T>::value,
1161  "dimDomain and number of children has to be the same");
1162  }
1163 
1169  VectorDiscreteGridFunctionDiv (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_)
1170  : pgfs(gfs)
1171  , lfs(*gfs)
1172  , lfs_cache(lfs)
1173  , x_view(*x_)
1174  , xl(pgfs->maxLocalSize())
1175  , J(pgfs->maxLocalSize())
1176  {
1177  static_assert(LBTraits::dimDomain == TypeTree::StaticDegree<T>::value,
1178  "dimDomain and number of children has to be the same");
1179  }
1180 
1181  inline void evaluate(const typename Traits::ElementType& e,
1182  const typename Traits::DomainType& x,
1183  typename Traits::RangeType& y) const
1184  {
1185  // get and bind local functions space
1186  lfs.bind(e);
1187  lfs_cache.update();
1188  x_view.bind(lfs_cache);
1189  x_view.read(xl);
1190  x_view.unbind();
1191 
1192  // get Jacobian of geometry
1193  const typename Traits::ElementType::Geometry::JacobianInverseTransposed
1194  JgeoIT = e.geometry().jacobianInverseTransposed(x);
1195 
1196  const typename Traits::ElementType::Geometry::JacobianInverseTransposed::size_type N =
1197  Traits::ElementType::Geometry::JacobianInverseTransposed::rows;
1198 
1199  y = 0.0;
1200 
1201  // loop over VectorGFS and calculate k-th derivative of k-th child
1202  for(unsigned int k=0; k != TypeTree::degree(lfs); ++k) {
1203 
1204  // get local Jacobians/gradients of the shape functions
1205  std::vector<typename LBTraits::JacobianType> J(lfs.child(k).size());
1206  lfs.child(k).finiteElement().localBasis().evaluateJacobian(x,J);
1207 
1208  RF d_k_phi;
1209  for(typename LFS::Traits::SizeType i=0; i<lfs.child(k).size(); i++) {
1210  // compute k-th derivative of k-th child
1211  d_k_phi =
1212  SingleDerivativeComputationHelper<
1213  typename Traits::ElementType::Geometry::JacobianInverseTransposed,
1214  typename Traits::ElementType::Geometry::JacobianInverseTransposed::field_type,
1215  N>::template compute_derivative<typename LBTraits::JacobianType::row_type>
1216  (JgeoIT,J[i][0],k);
1217 
1218  y += xl[lfs.child(k).localIndex(i)] * d_k_phi;
1219  }
1220  }
1221  }
1222 
1224  inline const typename Traits::GridViewType& getGridView() const
1225  {
1226  return pgfs->gridView();
1227  }
1228 
1229  private :
1230  typedef Dune::PDELab::LocalFunctionSpace<GFS> LFS;
1231  typedef Dune::PDELab::LFSIndexCache<LFS> LFSCache;
1232  typedef typename X::template ConstLocalView<LFSCache> XView;
1233 
1234  std::shared_ptr<GFS const> pgfs;
1235  mutable LFS lfs;
1236  mutable LFSCache lfs_cache;
1237  mutable XView x_view;
1238  mutable std::vector<RF> xl;
1239  mutable std::vector<JT> J;
1240  std::shared_ptr<const X> px;
1241  }; // end class VectorDiscreteGridFunctionDiv
1242 
1255  template<typename T, typename X, std::size_t dimR = TypeTree::StaticDegree<T>::value>
1256  class VectorDiscreteGridFunctionCurl
1257  {
1258  typedef T GFS;
1259  public :
1260 
1261  VectorDiscreteGridFunctionCurl(const GFS& gfs, const X& x)
1262  {
1263  static_assert(AlwaysFalse<typename GFS::Traits::GridViewType>::value,
1264  "Curl computation can only be done in two or three dimensions");
1265  }
1266 
1267  VectorDiscreteGridFunctionCurl(std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_)
1268  {
1269  static_assert(AlwaysFalse<typename GFS::Traits::GridViewType>::value,
1270  "Curl computation can only be done in two or three dimensions");
1271  }
1272 
1273  };
1274 
1285  template<typename T, typename X>
1286  class VectorDiscreteGridFunctionCurl<T,X,3>
1287  : public Dune::PDELab::GridFunctionBase<
1288  Dune::PDELab::GridFunctionTraits<
1289  typename T::Traits::GridViewType,
1290  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
1291  //T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimRange,
1292  TypeTree::StaticDegree<T>::value,
1293  Dune::FieldVector<
1294  typename T::template Child<0>::Type::Traits::FiniteElementType
1295  ::Traits::LocalBasisType::Traits::RangeFieldType,
1296  //T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimRange
1297  TypeTree::StaticDegree<T>::value
1298  >
1299  >,
1300  VectorDiscreteGridFunctionCurl<T,X>
1301  >
1302  {
1303  typedef T GFS;
1304 
1305  typedef Dune::PDELab::GridFunctionBase<
1306  Dune::PDELab::GridFunctionTraits<
1307  typename T::Traits::GridViewType,
1308  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
1309  //T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimRange,
1310  TypeTree::StaticDegree<T>::value,
1311  Dune::FieldVector<
1312  typename T::template Child<0>::Type::Traits::FiniteElementType
1313  ::Traits::LocalBasisType::Traits::RangeFieldType,
1314  //T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimRange
1315  TypeTree::StaticDegree<T>::value
1316  >
1317  >,
1318  VectorDiscreteGridFunctionCurl<T,X> > BaseT;
1319 
1320  public :
1321  typedef typename BaseT::Traits Traits;
1322  typedef typename T::template Child<0>::Type ChildType;
1323  typedef typename ChildType::Traits::FiniteElementType::Traits::LocalBasisType::Traits LBTraits;
1324 
1325  typedef typename LBTraits::RangeFieldType RF;
1326  typedef typename LBTraits::JacobianType JT;
1327 
1333  VectorDiscreteGridFunctionCurl(const GFS& gfs, const X& x_)
1334  : pgfs(stackobject_to_shared_ptr(gfs))
1335  , lfs(gfs)
1336  , lfs_cache(lfs)
1337  , x_view(x_)
1338  , xl(gfs.maxLocalSize())
1339  , J(gfs.maxLocalSize())
1340  {
1341  static_assert(LBTraits::dimDomain == TypeTree::StaticDegree<T>::value,
1342  "dimDomain and number of children has to be the same");
1343  }
1344 
1350  VectorDiscreteGridFunctionCurl (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_)
1351  : pgfs(gfs)
1352  , lfs(*gfs)
1353  , lfs_cache(lfs)
1354  , x_view(*x_)
1355  , xl(pgfs->maxLocalSize())
1356  , J(pgfs->maxLocalSize())
1357  {
1358  static_assert(LBTraits::dimDomain == TypeTree::StaticDegree<T>::value,
1359  "dimDomain and number of children has to be the same");
1360  }
1361 
1362  inline void evaluate(const typename Traits::ElementType& e,
1363  const typename Traits::DomainType& x,
1364  typename Traits::RangeType& y) const
1365  {
1366  // get and bind local functions space
1367  lfs.bind(e);
1368  lfs_cache.update();
1369  x_view.bind(lfs_cache);
1370  x_view.read(xl);
1371  x_view.unbind();
1372 
1373  // get Jacobian of geometry
1374  const typename Traits::ElementType::Geometry::JacobianInverseTransposed
1375  JgeoIT = e.geometry().jacobianInverseTransposed(x);
1376 
1377  const typename Traits::ElementType::Geometry::JacobianInverseTransposed::size_type N =
1378  Traits::ElementType::Geometry::JacobianInverseTransposed::rows;
1379 
1380  y = 0.0;
1381 
1382  // some handy variables for the curl in 3D
1383  int i1, i2;
1384 
1385  // loop over childs of VectorGFS
1386  for(unsigned int k=0; k != TypeTree::degree(lfs); ++k) {
1387 
1388  // get local Jacobians/gradients of the shape functions
1389  std::vector<typename LBTraits::JacobianType> J(lfs.child(k).size());
1390  lfs.child(k).finiteElement().localBasis().evaluateJacobian(x,J);
1391 
1392  // pick out the right derivative and component for the curl computation
1393  i1 = (k+1)%3;
1394  i2 = (k+2)%3;
1395 
1396  RF d_k_phi;
1397  for(typename LFS::Traits::SizeType i=0; i<lfs.child(k).size(); i++) {
1398  // compute i2-th derivative of k-th child
1399  d_k_phi =
1400  SingleDerivativeComputationHelper<
1401  typename Traits::ElementType::Geometry::JacobianInverseTransposed,
1402  typename Traits::ElementType::Geometry::JacobianInverseTransposed::field_type,
1403  N>::template compute_derivative<typename LBTraits::JacobianType::row_type>
1404  (JgeoIT,J[i][0],i2);
1405 
1406  y[i1] += xl[lfs.child(k).localIndex(i)] * d_k_phi;
1407 
1408  // compute i1-th derivative of k-th child
1409  d_k_phi =
1410  SingleDerivativeComputationHelper<
1411  typename Traits::ElementType::Geometry::JacobianInverseTransposed,
1412  typename Traits::ElementType::Geometry::JacobianInverseTransposed::field_type,
1413  N>::template compute_derivative<typename LBTraits::JacobianType::row_type>
1414  (JgeoIT,J[i][0],i1);
1415 
1416  y[i2] -= xl[lfs.child(k).localIndex(i)] * d_k_phi;
1417  }
1418  }
1419  }
1420 
1422  inline const typename Traits::GridViewType& getGridView() const
1423  {
1424  return pgfs->gridView();
1425  }
1426 
1427  private :
1428  typedef Dune::PDELab::LocalFunctionSpace<GFS> LFS;
1429  typedef Dune::PDELab::LFSIndexCache<LFS> LFSCache;
1430  typedef typename X::template ConstLocalView<LFSCache> XView;
1431 
1432  std::shared_ptr<GFS const> pgfs;
1433  mutable LFS lfs;
1434  mutable LFSCache lfs_cache;
1435  mutable XView x_view;
1436  mutable std::vector<RF> xl;
1437  mutable std::vector<JT> J;
1438  std::shared_ptr<const X> px;
1439  }; // end class VectorDiscreteGridFunctionCurl (3D)
1440 
1451  template<typename T, typename X>
1452  class VectorDiscreteGridFunctionCurl<T,X,2>
1453  : public Dune::PDELab::GridFunctionBase<
1454  Dune::PDELab::GridFunctionTraits<
1455  typename T::Traits::GridViewType,
1456  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
1457  T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimRange,
1458  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType>,
1459  VectorDiscreteGridFunctionDiv<T,X> >
1460  {
1461  typedef T GFS;
1462 
1463  typedef Dune::PDELab::GridFunctionBase<
1464  Dune::PDELab::GridFunctionTraits<
1465  typename T::Traits::GridViewType,
1466  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
1467  T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::dimRange,
1468  typename T::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType>,
1469  VectorDiscreteGridFunctionDiv<T,X> > BaseT;
1470  public :
1471  typedef typename BaseT::Traits Traits;
1472  typedef typename T::template Child<0>::Type ChildType;
1473  typedef typename ChildType::Traits::FiniteElementType::Traits::LocalBasisType::Traits LBTraits;
1474 
1475  typedef typename LBTraits::RangeFieldType RF;
1476  typedef typename LBTraits::JacobianType JT;
1477 
1483  VectorDiscreteGridFunctionCurl(const GFS& gfs, const X& x_)
1484  : pgfs(stackobject_to_shared_ptr(gfs))
1485  , lfs(gfs)
1486  , lfs_cache(lfs)
1487  , x_view(x_)
1488  , xl(gfs.maxLocalSize())
1489  , J(gfs.maxLocalSize())
1490  {
1491  static_assert(LBTraits::dimDomain == TypeTree::StaticDegree<T>::value,
1492  "dimDomain and number of children has to be the same");
1493  }
1494 
1500  VectorDiscreteGridFunctionCurl (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_)
1501  : pgfs(gfs)
1502  , lfs(*gfs)
1503  , lfs_cache(lfs)
1504  , x_view(*x_)
1505  , xl(pgfs->maxLocalSize())
1506  , J(pgfs->maxLocalSize())
1507  {
1508  static_assert(LBTraits::dimDomain == TypeTree::StaticDegree<T>::value,
1509  "dimDomain and number of children has to be the same");
1510  }
1511 
1512  inline void evaluate(const typename Traits::ElementType& e,
1513  const typename Traits::DomainType& x,
1514  typename Traits::RangeType& y) const
1515  {
1516  // get and bind local functions space
1517  lfs.bind(e);
1518  lfs_cache.update();
1519  x_view.bind(lfs_cache);
1520  x_view.read(xl);
1521  x_view.unbind();
1522 
1523  // get Jacobian of geometry
1524  const typename Traits::ElementType::Geometry::JacobianInverseTransposed
1525  JgeoIT = e.geometry().jacobianInverseTransposed(x);
1526 
1527  const typename Traits::ElementType::Geometry::JacobianInverseTransposed::size_type N =
1528  Traits::ElementType::Geometry::JacobianInverseTransposed::rows;
1529 
1530  y = 0.0;
1531 
1532  // some handy variables for the curl computation in 2D
1533  RF sign = -1.0;
1534  int i2;
1535 
1536  // loop over childs of VectorGFS
1537  for(unsigned int k=0; k != TypeTree::degree(lfs); ++k) {
1538 
1539  // get local Jacobians/gradients of the shape functions
1540  std::vector<typename LBTraits::JacobianType> J(lfs.child(k).size());
1541  lfs.child(k).finiteElement().localBasis().evaluateJacobian(x,J);
1542 
1543  RF d_k_phi;
1544  // pick out the right derivative
1545  i2 = 1-k;
1546  for(typename LFS::Traits::SizeType i=0; i<lfs.child(k).size(); i++) {
1547  // compute i2-th derivative of k-th child
1548  d_k_phi =
1549  SingleDerivativeComputationHelper<
1550  typename Traits::ElementType::Geometry::JacobianInverseTransposed,
1551  typename Traits::ElementType::Geometry::JacobianInverseTransposed::field_type,
1552  N>::template compute_derivative<typename LBTraits::JacobianType::row_type>
1553  (JgeoIT,J[i][0],i2);
1554 
1555  y += sign * xl[lfs.child(k).localIndex(i)] * d_k_phi;
1556  }
1557  sign *= -1.0;
1558  }
1559  }
1560 
1562  inline const typename Traits::GridViewType& getGridView() const
1563  {
1564  return pgfs->gridView();
1565  }
1566 
1567  private :
1568  typedef Dune::PDELab::LocalFunctionSpace<GFS> LFS;
1569  typedef Dune::PDELab::LFSIndexCache<LFS> LFSCache;
1570  typedef typename X::template ConstLocalView<LFSCache> XView;
1571 
1572  std::shared_ptr<GFS const> pgfs;
1573  mutable LFS lfs;
1574  mutable LFSCache lfs_cache;
1575  mutable XView x_view;
1576  mutable std::vector<RF> xl;
1577  mutable std::vector<JT> J;
1578  std::shared_ptr<const X> px;
1579  }; // end class VectorDiscreteGridFunctionCurl (2D)
1580 
1582  } // namespace PDELab
1583 } // namespace Dune
1584 
1585 #endif // DUNE_PDELAB_GRIDFUNCTIONSPACE_GRIDFUNCTIONSPACEUTILITIES_HH
Dune::DenseVector::axpy
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:576
Dune::DiagonalMatrix
A diagonal matrix of static size.
Definition: diagonalmatrix.hh:51
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:69
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:95
Dune::PDELab::DiscreteGridFunctionCurl
convert a grid function space and a coefficient vector into a grid function of the curl
Definition: gridfunctionspaceutilities.hh:193
Dune::PDELab::DiscreteGridFunctionCurl::getGridView
const Traits::GridViewType & getGridView() const
get a reference to the GridView
Definition: gridfunctionspaceutilities.hh:260
Dune::PDELab::DiscreteGridFunctionCurl::DiscreteGridFunctionCurl
DiscreteGridFunctionCurl(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_)
Construct a DiscreteGridFunctionCurl.
Definition: gridfunctionspaceutilities.hh:226
Dune::PDELab::DiscreteGridFunctionCurl::DiscreteGridFunctionCurl
DiscreteGridFunctionCurl(const GFS &gfs, const X &x_)
Construct a DiscreteGridFunctionCurl.
Definition: gridfunctionspaceutilities.hh:210
Dune::PDELab::DiscreteGridFunctionGlobalCurl
convert a single component function space with experimental global finite elements into a grid functi...
Definition: gridfunctionspaceutilities.hh:375
Dune::PDELab::DiscreteGridFunctionGlobalCurl::DiscreteGridFunctionGlobalCurl
DiscreteGridFunctionGlobalCurl(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_)
Construct a DiscreteGridFunctionGlobalCurl.
Definition: gridfunctionspaceutilities.hh:413
Dune::PDELab::DiscreteGridFunctionGlobalCurl::DiscreteGridFunctionGlobalCurl
DiscreteGridFunctionGlobalCurl(const GFS &gfs, const X &x_)
Construct a DiscreteGridFunctionGlobalCurl.
Definition: gridfunctionspaceutilities.hh:398
Dune::PDELab::DiscreteGridFunctionGlobalCurl::getGridView
const Traits::GridViewType & getGridView() const
get a reference to the GridView
Definition: gridfunctionspaceutilities.hh:461
Dune::PDELab::DiscreteGridFunctionGradient
convert a single component function space with a grid function representing the gradient
Definition: gridfunctionspaceutilities.hh:505
Dune::PDELab::DiscreteGridFunctionGradient::getGridView
const Traits::GridViewType & getGridView() const
get a reference to the GridView
Definition: gridfunctionspaceutilities.hh:588
Dune::PDELab::DiscreteGridFunctionGradient::DiscreteGridFunctionGradient
DiscreteGridFunctionGradient(const GFS &gfs, const X &x_)
Construct a DiscreteGridFunctionGradient.
Definition: gridfunctionspaceutilities.hh:530
Dune::PDELab::DiscreteGridFunctionGradient::DiscreteGridFunctionGradient
DiscreteGridFunctionGradient(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_)
Construct a DiscreteGridFunctionGradient.
Definition: gridfunctionspaceutilities.hh:543
Dune::PDELab::DiscreteGridFunctionPiola
DiscreteGridFunction with Piola transformation.
Definition: gridfunctionspaceutilities.hh:620
Dune::PDELab::DiscreteGridFunctionPiola::DiscreteGridFunctionPiola
DiscreteGridFunctionPiola(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_)
Construct a DiscreteGridFunctionPiola.
Definition: gridfunctionspaceutilities.hh:655
Dune::PDELab::DiscreteGridFunctionPiola::getGridView
const Traits::GridViewType & getGridView() const
get a reference to the GridView
Definition: gridfunctionspaceutilities.hh:692
Dune::PDELab::DiscreteGridFunctionPiola::DiscreteGridFunctionPiola
DiscreteGridFunctionPiola(const GFS &gfs, const X &x_)
Construct a DiscreteGridFunctionPiola.
Definition: gridfunctionspaceutilities.hh:640
Dune::PDELab::DiscreteGridFunction
convert a grid function space and a coefficient vector into a grid function
Definition: gridfunctionspaceutilities.hh:76
Dune::PDELab::DiscreteGridFunction::DiscreteGridFunction
DiscreteGridFunction(const GFS &gfs, const X &x_)
Construct a DiscreteGridFunction.
Definition: gridfunctionspaceutilities.hh:102
Dune::PDELab::DiscreteGridFunction::DiscreteGridFunction
DiscreteGridFunction(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_)
Construct a DiscreteGridFunction.
Definition: gridfunctionspaceutilities.hh:117
Dune::PDELab::DiscreteGridFunction::getGridView
const Traits::GridViewType & getGridView() const
get a reference to the GridView
Definition: gridfunctionspaceutilities.hh:150
Dune::PDELab::GridFunctionBase
leaf of a function tree
Definition: function.hh:302
Dune::PDELab::GridFunctionInterface::Traits
T Traits
Export type traits.
Definition: function.hh:193
Dune::PDELab::JacobianToCurl
extract the curl of a function from the jacobian of that function
Definition: jacobiantocurl.hh:27
Dune::PDELab::LocalFunctionSpace
Create a local function space from a global function space.
Definition: localfunctionspace.hh:754
Dune::PDELab::VectorDiscreteGridFunctionCurl< T, X, 2 >::getGridView
const Traits::GridViewType & getGridView() const
get a reference to the GridView
Definition: gridfunctionspaceutilities.hh:1562
Dune::PDELab::VectorDiscreteGridFunctionCurl< T, X, 2 >::VectorDiscreteGridFunctionCurl
VectorDiscreteGridFunctionCurl(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_)
Construct a VectorDiscreteGridFunctionCurl.
Definition: gridfunctionspaceutilities.hh:1500
Dune::PDELab::VectorDiscreteGridFunctionCurl< T, X, 2 >::VectorDiscreteGridFunctionCurl
VectorDiscreteGridFunctionCurl(const GFS &gfs, const X &x_)
Construct a VectorDiscreteGridFunctionCurl.
Definition: gridfunctionspaceutilities.hh:1483
Dune::PDELab::VectorDiscreteGridFunctionCurl< T, X, 3 >::VectorDiscreteGridFunctionCurl
VectorDiscreteGridFunctionCurl(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_)
Construct a VectorDiscreteGridFunctionCurl.
Definition: gridfunctionspaceutilities.hh:1350
Dune::PDELab::VectorDiscreteGridFunctionCurl< T, X, 3 >::getGridView
const Traits::GridViewType & getGridView() const
get a reference to the GridView
Definition: gridfunctionspaceutilities.hh:1422
Dune::PDELab::VectorDiscreteGridFunctionCurl< T, X, 3 >::VectorDiscreteGridFunctionCurl
VectorDiscreteGridFunctionCurl(const GFS &gfs, const X &x_)
Construct a VectorDiscreteGridFunctionCurl.
Definition: gridfunctionspaceutilities.hh:1333
Dune::PDELab::VectorDiscreteGridFunctionCurl
Compute curl of vector-valued functions.
Definition: gridfunctionspaceutilities.hh:1257
Dune::PDELab::VectorDiscreteGridFunctionDiv
Compute divergence of vector-valued functions.
Definition: gridfunctionspaceutilities.hh:1129
Dune::PDELab::VectorDiscreteGridFunctionDiv::VectorDiscreteGridFunctionDiv
VectorDiscreteGridFunctionDiv(const GFS &gfs, const X &x_)
Construct a VectorDiscreteGridFunctionDiv.
Definition: gridfunctionspaceutilities.hh:1152
Dune::PDELab::VectorDiscreteGridFunctionDiv::VectorDiscreteGridFunctionDiv
VectorDiscreteGridFunctionDiv(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_)
Construct a VectorDiscreteGridFunctionDiv.
Definition: gridfunctionspaceutilities.hh:1169
Dune::PDELab::VectorDiscreteGridFunctionDiv::getGridView
const Traits::GridViewType & getGridView() const
get a reference to the GridView
Definition: gridfunctionspaceutilities.hh:1224
Dune::PDELab::VectorDiscreteGridFunctionGradient
Equivalent of DiscreteGridFunctionGradient for vector-valued functions.
Definition: gridfunctionspaceutilities.hh:923
Dune::PDELab::VectorDiscreteGridFunctionGradient::getGridView
const Traits::GridViewType & getGridView() const
get a reference to the GridView
Definition: gridfunctionspaceutilities.hh:1015
Dune::PDELab::VectorDiscreteGridFunctionGradient::VectorDiscreteGridFunctionGradient
VectorDiscreteGridFunctionGradient(const GFS &gfs, const X &x_)
Construct a VectorDiscreteGridFunctionGradient.
Definition: gridfunctionspaceutilities.hh:953
Dune::PDELab::VectorDiscreteGridFunctionGradient::VectorDiscreteGridFunctionGradient
VectorDiscreteGridFunctionGradient(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_)
Construct a VectorDiscreteGridFunctionGradient.
Definition: gridfunctionspaceutilities.hh:968
Dune::PDELab::VectorDiscreteGridFunction
DiscreteGridFunction for vector-valued functions.
Definition: gridfunctionspaceutilities.hh:744
Dune::PDELab::VectorDiscreteGridFunction::VectorDiscreteGridFunction
VectorDiscreteGridFunction(const GFS &gfs, const X &x_, const Remap &remap_)
construct
Definition: gridfunctionspaceutilities.hh:821
Dune::PDELab::VectorDiscreteGridFunction::getGridView
const Traits::GridViewType & getGridView() const
get a reference to the GridView
Definition: gridfunctionspaceutilities.hh:882
Dune::PDELab::VectorDiscreteGridFunction::VectorDiscreteGridFunction
VectorDiscreteGridFunction(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_, const Remap &remap_)
Construct a VectorDiscreteGridFunction.
Definition: gridfunctionspaceutilities.hh:848
Dune::PDELab::VectorDiscreteGridFunction::VectorDiscreteGridFunction
VectorDiscreteGridFunction(std::shared_ptr< const GFS > gfs, std::shared_ptr< const X > x_, std::size_t start=0)
Construct a VectorDiscreteGridFunction.
Definition: gridfunctionspaceutilities.hh:795
Dune::PDELab::VectorDiscreteGridFunction::VectorDiscreteGridFunction
VectorDiscreteGridFunction(const GFS &gfs, const X &x_, std::size_t start=0)
construct
Definition: gridfunctionspaceutilities.hh:776
exceptions.hh
A few common exception classes.
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune::TypeTree::degree
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:76
Dune::TypeTree::StaticDegree
decltype(Node::degree()) StaticDegree
Returns the statically known degree of the given Node type as a std::integral_constant.
Definition: nodeinterface.hh:104
Dune
Dune namespace.
Definition: alignedallocator.hh:11
Dune::stackobject_to_shared_ptr
std::shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:70
Dune::sign
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:177
shared_ptr.hh
This file implements several utilities related to std::shared_ptr.
Dune::AlwaysFalse
template which always yields a false value
Definition: typetraits.hh:124
Dune::BasisInterfaceSwitch
Switch for uniform treatment of local and global basis classes.
Definition: interfaceswitch.hh:152
Dune::BasisInterfaceSwitch::dimRange
static const std::size_t dimRange
export dimension of the values
Definition: interfaceswitch.hh:163
Dune::BasisInterfaceSwitch::RangeField
Basis::Traits::RangeField RangeField
export field type of the values
Definition: interfaceswitch.hh:161
Dune::BasisInterfaceSwitch::Range
Basis::Traits::Range Range
export vector type of the values
Definition: interfaceswitch.hh:165
Dune::FiniteElementInterfaceSwitch
Switch for uniform treatment of finite element with either the local or the global interface.
Definition: interfaceswitch.hh:28
Dune::PDELab::DiscreteGridFunctionCurlTraits
Helper class to calculate the Traits of DiscreteGridFunctionCurl.
Definition: gridfunctionspaceutilities.hh:294
Dune::PDELab::FunctionTraits< GV::Grid::ctype, GV::dimension, Dune::FieldVector< GV::Grid::ctype, GV::dimension >, RF, m, R >::DomainType
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:50
Dune::PDELab::FunctionTraits::RangeFieldType
RF RangeFieldType
Export type for range field.
Definition: function.hh:53
Dune::PDELab::GridFunctionTraits
traits class holding the function signature, same as in local function
Definition: function.hh:183
Dune::PDELab::PowerCompositeGridFunctionTraits::ElementType
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:119
Dune::PDELab::PowerCompositeGridFunctionTraits::GridViewType
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:116
Dune::PDELab::SingleDerivativeComputationHelper< Dune::DiagonalMatrix< RF, size >, RF, size >::compute_derivative
static RF compute_derivative(const Dune::DiagonalMatrix< RF, size > &mat, const T &t, const unsigned int k)
Definition: gridfunctionspaceutilities.hh:1104
Dune::PDELab::SingleDerivativeComputationHelper< Dune::FieldMatrix< RF, size, size >, RF, size >::compute_derivative
static RF compute_derivative(const Dune::FieldMatrix< RF, size, size > &mat, const T &t, const unsigned int k)
Definition: gridfunctionspaceutilities.hh:1080
Dune::PDELab::SingleDerivativeComputationHelper
Helper class to compute a single derivative of scalar basis functions.
Definition: gridfunctionspaceutilities.hh:1048
Dune::PDELab::SingleDerivativeComputationHelper::compute_derivative
static RF compute_derivative(const Mat &mat, const T &t, const unsigned int k)
Definition: gridfunctionspaceutilities.hh:1054
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (May 16, 22:29, 2024)