dune-pdelab  2.5-dev
discretegridviewfunction.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 #ifndef DUNE_PDELAB_FUNCTION_DISCRETEGRIDVIEWFUNCTION_HH
4 #define DUNE_PDELAB_FUNCTION_DISCRETEGRIDVIEWFUNCTION_HH
5 
6 #include <cstdlib>
7 #include <vector>
8 #include <memory>
9 #include <type_traits>
10 
11 #include <dune/common/exceptions.hh>
12 #include <dune/common/fvector.hh>
13 
14 #include <dune/localfunctions/common/interfaceswitch.hh>
15 
20 
21 #include <dune/functions/gridfunctions/gridviewfunction.hh>
22 
23 namespace std {
25  template<typename T, int N, typename R2>
26  struct common_type<Dune::FieldVector<T,N>, R2>
27  {
28  using type = Dune::FieldVector<typename std::common_type<T,R2>::type,N>;
29  };
31  template<typename T, int N, typename R2>
32  struct common_type<Dune::FieldVector<T,N>, Dune::FieldVector<R2,N>>
33  {
34  using type = Dune::FieldVector<typename std::common_type<T,R2>::type,N>;
35  };
36 }
37 
38 namespace Dune {
39 namespace PDELab {
40 
41 template<typename Signature, typename E, template<class> class D, int B,
42  int diffOrder>
44  Functions::Imp::GridFunctionTraits<
45  typename DiscreteGridViewFunctionTraits<Signature,E,D,B,diffOrder-1>::DerivativeSignature
46  ,E,D,B>
47 {};
48 
49 template<typename Signature, typename E, template<class> class D, int B>
50 struct DiscreteGridViewFunctionTraits<Signature,E,D,B,0> :
51  Functions::Imp::GridFunctionTraits<Signature,E,D,B>
52 {};
53 
68 template<typename GFS, typename V, int diffOrder = 0>
70 {
71 public:
72  using GridView = typename GFS::Traits::GridView;
73  using EntitySet = Functions::GridViewEntitySet<GridView, 0>;
74 
75  using Domain = typename EntitySet::GlobalCoordinate;
76  using LocalBasisTraits = typename GFS::Traits::FiniteElementMap::Traits::FiniteElement::Traits::LocalBasisType::Traits;
77  using LocalBasisRange = typename LocalBasisTraits::RangeType;
78  using VectorRange = typename V::ElementType;
79  using ElementaryRange = typename std::common_type<LocalBasisRange, VectorRange>::type;
80 
81  using LocalDomain = typename EntitySet::LocalCoordinate;
82  using Element = typename EntitySet::Element;
83 
84  using Traits =
86  Functions::DefaultDerivativeTraits, 16, diffOrder>;
87 
88  using Range = typename Traits::Range; // this is actually either the Range, the Jacobian or Hessian
89 
90  using Basis = GFS;
91  using GridFunctionSpace = GFS;
92  using Vector = V;
93  enum { maxDiffOrder = LocalBasisTraits::diffOrder - diffOrder };
94 
96  {
99  using XView = typename Vector::template ConstLocalView<LFSCache>;
100 
101  public:
102 
107  using size_type = std::size_t;
108 
109  enum { maxDiffOrder = LocalBasisTraits::diffOrder - diffOrder };
110 
111  LocalFunction(const shared_ptr<const GridFunctionSpace> gfs, const shared_ptr<const Vector> v)
112  : pgfs_(gfs)
113  , v_(v)
114  , lfs_(*pgfs_)
115  , lfs_cache_(lfs_)
116  , x_view_(*v_)
117  , xl_(pgfs_->maxLocalSize())
118  , yb_(pgfs_->maxLocalSize())
119  , element_(nullptr)
120  {}
121 
128  void bind(const Element& element)
129  {
130  element_ = &element;
131  lfs_.bind(element);
132  lfs_cache_.update();
133  x_view_.bind(lfs_cache_);
134  x_view_.read(xl_);
135  x_view_.unbind();
136  }
137 
138  void unbind()
139  {
140  element_ = nullptr;
141  }
142 
143  const Element& localContext() const
144  {
145 #ifndef NDEBUG
146  if (!element_)
147  DUNE_THROW(InvalidStateException,"can't call localContext on unbound DiscreteGridViewFunction::LocalFunction");
148 #endif
149  return *element_;
150  }
151 
168  {
170  diff(t.pgfs_, t.v_);
171  // TODO: do we really want this?
172  if (t.element_) diff.bind(*t.element_);
173  return diff;
174  }
175 
185  Range
186  operator()(const Domain& coord)
187  {
188  return evaluate<LocalBasisTraits::diffOrder, diffOrder>(coord);
189  };
190 
191  private:
192  template<int maxDiffOrder, int dOrder>
193  typename std::enable_if<(dOrder > 2 or dOrder > maxDiffOrder),
194  Range>::type
195  evaluate(const Domain& coord) const
196  {
197  if (diffOrder > 2) DUNE_THROW(NotImplemented,
198  "Derivatives are only implemented up to degree 2");
199  if (diffOrder > maxDiffOrder) DUNE_THROW(NotImplemented,
200  "Derivative of degree " << diffOrder << " is not provided by the local basis");
201  };
202 
203  template<int maxDiffOrder, int dOrder>
204  typename std::enable_if<dOrder == 0,
205  Range>::type
206  evaluate(const Domain& coord) const
207  {
208  Range r(0);
209  auto& basis = lfs_.finiteElement().localBasis();
210  basis.evaluateFunction(coord,yb_);
211  for (size_type i = 0; i < yb_.size(); ++i)
212  {
213  r.axpy(xl_[i],yb_[i]);
214  }
215  return r;
216  }
217 
218  template<int maxDiffOrder, int dOrder>
219  typename std::enable_if<maxDiffOrder >= 1 and dOrder == 1,
220  Range>::type
221  evaluate(const Domain& coord) const
222  {
223  Range r(0);
224  // get Jacobian of geometry
225  const typename Element::Geometry::JacobianInverseTransposed
226  JgeoIT = element_->geometry().jacobianInverseTransposed(coord);
227 
228  // get local Jacobians/gradients of the shape functions
229  lfs_.finiteElement().localBasis().evaluateJacobian(coord,yb_);
230 
231  Range gradphi;
232  r = 0;
233  for(std::size_t i = 0; i < yb_.size(); ++i) {
234  assert(gradphi.size() == yb_[i].size());
235  for(std::size_t j = 0; j < gradphi.size(); ++j) {
236  // compute global gradient of shape function i
237  // graphi += {J^{-1}}^T * yb_i0
238  JgeoIT.mv(yb_[i][j], gradphi[j]);
239 
240  // sum up global gradients, weighting them with the appropriate coeff
241  // r \in R^{1,dim}
242  // r_0 += xl_i * grad \phi
243  r[j].axpy(xl_[i], gradphi[j]);
244  }
245  }
246  return r;
247  }
248 
249  template<int maxDiffOrder, int dOrder>
250  typename std::enable_if<maxDiffOrder >= 2 and dOrder == 2,
251  Range>::type
252  evaluate(const Domain& coord) const
253  {
254  Range r(0);
255  // TODO: we currently require affine geometries.
256  if (! element_->geometry().affine())
257  DUNE_THROW(NotImplemented, "Due to missing features in the Geometry interface, "
258  "the computation of higher derivatives (>=2) works only for affine transformations.");
259  // get Jacobian of geometry
260  const typename Element::Geometry::JacobianInverseTransposed
261  JgeoIT = element_->geometry().jacobianInverseTransposed(coord);
262 
263  // TODO: we currently only implement the hessian...
264  // a proper implementation will require TMP magic.
265  static const unsigned int dim = GridView::dimensionworld;
266  // static_assert(
267  // isHessian<Range>::value,
268  // "Currently the only higher order derivative we support is the Hessian of scalar functions");
269 
270  // get local hessian of the shape functions
271  array<std::size_t, dim> directions;
272  for(std::size_t i = 0; i < dim; ++i) {
273  // evaluate diagonal entries
274  directions[i] = 2;
275  // evaluate offdiagonal entries
276  directions[i] = 1;
277  for(std::size_t j = i+1; j < dim; ++j) {
278  directions[j] = 1;
279  lfs_.finiteElement().localBasis().partial(directions,coord,yb_);
280  assert( yb_.size() == 1); // TODO: we only implement the hessian of scalar functions
281  for(std::size_t n = 0; n < yb_.size(); ++n) {
282  // sum up derivatives, weighting them with the appropriate coeff
283  r[i][j] += xl_[i] * yb_[j];
284  }
285  // use symmetry of the hessian
286  r[j][i] = r[i][j];
287  directions[j] = 0;
288  }
289  directions[i] = 0;
290  }
291  // transform back to global coordinates
292  for(std::size_t i = 0; i < dim; ++i)
293  for(std::size_t j = i; j < dim; ++j)
294  r[i][j] *= JgeoIT[i][j] * JgeoIT[i][j];
295  return r;
296  }
297 
298  protected:
299 
300  const shared_ptr<const GridFunctionSpace> pgfs_;
301  const shared_ptr<const Vector> v_;
304  XView x_view_;
305  mutable std::vector<ElementaryRange> xl_;
306  mutable std::vector<Range> yb_;
308  };
309 
311  : pgfs_(stackobject_to_shared_ptr(gfs)),v_(stackobject_to_shared_ptr(v))
312  {}
313 
314  DiscreteGridViewFunction(std::shared_ptr<const GridFunctionSpace> pgfs, std::shared_ptr<const Vector> v)
315  : pgfs_(pgfs),v_(v)
316  {}
317 
318  // this is part of the interface in dune-functions
319  const Basis& basis() const
320  {
321  return *pgfs_;
322  }
324  {
325  return *pgfs_;
326  }
327 
328  const V& dofs() const
329  {
330  return *v_;
331  }
332 
333  // TODO: Implement this using hierarchic search
334  Range operator() (const Domain& x) const
335  {
336  DUNE_THROW(NotImplemented,"not implemented");
337  }
338 
340  {
341  return DiscreteGridViewFunction<GFS,V,diffOrder+1>(t.pgfs_, t.v_);
342  }
343 
354  {
355  return LocalFunction(t.pgfs_, t.v_);
356  }
357 
361  EntitySet entitySet() const
362  {
363  return pgfs_->gridView();
364  }
365 
366 private:
367 
368  const shared_ptr<const GridFunctionSpace> pgfs_;
369  const shared_ptr<const Vector> v_;
370 
371 };
372 
373 } // end of namespace Dune::PDELab
374 } // end of namespace Dune
375 
376 #endif // DUNE_PDELAB_FUNCTION_DISCRETEGRIDVIEWFUNCTION_HH
static const int dim
Definition: adaptivity.hh:83
LFS lfs_
Definition: discretegridviewfunction.hh:302
GlobalFunction::Element Element
Definition: discretegridviewfunction.hh:106
GFS GridFunctionSpace
Definition: discretegridviewfunction.hh:91
typename EntitySet::LocalCoordinate LocalDomain
Definition: discretegridviewfunction.hh:81
typename std::common_type< LocalBasisRange, VectorRange >::type ElementaryRange
Definition: discretegridviewfunction.hh:79
const Element * element_
Definition: discretegridviewfunction.hh:307
std::size_t size_type
Definition: discretegridviewfunction.hh:107
const V & dofs() const
Definition: discretegridviewfunction.hh:328
void bind(const Element &element)
Bind LocalFunction to grid element.
Definition: discretegridviewfunction.hh:128
friend LocalFunction localFunction(const DiscreteGridViewFunction &t)
Get local function of wrapped function.
Definition: discretegridviewfunction.hh:353
LocalFunction(const shared_ptr< const GridFunctionSpace > gfs, const shared_ptr< const Vector > v)
Definition: discretegridviewfunction.hh:111
typename Traits::Range Range
Definition: discretegridviewfunction.hh:88
std::vector< ElementaryRange > xl_
Definition: discretegridviewfunction.hh:305
const GridFunctionSpace & gridFunctionSpace() const
Definition: discretegridviewfunction.hh:323
DiscreteGridViewFunction(std::shared_ptr< const GridFunctionSpace > pgfs, std::shared_ptr< const Vector > v)
Definition: discretegridviewfunction.hh:314
GlobalFunction::Range Range
Definition: discretegridviewfunction.hh:105
const Basis & basis() const
Definition: discretegridviewfunction.hh:319
Range operator()(const Domain &coord)
Evaluate LocalFunction at bound element.
Definition: discretegridviewfunction.hh:186
LFSCache lfs_cache_
Definition: discretegridviewfunction.hh:303
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
typename GFS::Traits::GridView GridView
Definition: discretegridviewfunction.hh:72
EntitySet entitySet() const
Get associated EntitySet.
Definition: discretegridviewfunction.hh:361
V Vector
Definition: discretegridviewfunction.hh:92
Dune::FieldVector< typename std::common_type< T, R2 >::type, N > type
Definition: discretegridviewfunction.hh:28
friend DiscreteGridViewFunction< GFS, V, diffOrder+1 >::LocalFunction derivative(const LocalFunction &t)
free function to obtain the derivative of a LocalFunction
Definition: discretegridviewfunction.hh:167
typename EntitySet::GlobalCoordinate Domain
Definition: discretegridviewfunction.hh:75
Dune::FieldVector< typename std::common_type< T, R2 >::type, N > type
Definition: discretegridviewfunction.hh:34
XView x_view_
Definition: discretegridviewfunction.hh:304
DiscreteGridViewFunction(const GridFunctionSpace &gfs, const Vector &v)
Definition: discretegridviewfunction.hh:310
void unbind()
Definition: discretegridviewfunction.hh:138
LocalDomain Domain
Definition: discretegridviewfunction.hh:104
STL namespace.
const shared_ptr< const GridFunctionSpace > pgfs_
Definition: discretegridviewfunction.hh:300
Definition: discretegridviewfunction.hh:95
typename V::ElementType VectorRange
Definition: discretegridviewfunction.hh:78
Definition: discretegridviewfunction.hh:43
const shared_ptr< const Vector > v_
Definition: discretegridviewfunction.hh:301
typename EntitySet::Element Element
Definition: discretegridviewfunction.hh:82
typename LocalBasisTraits::RangeType LocalBasisRange
Definition: discretegridviewfunction.hh:77
const Element & localContext() const
Definition: discretegridviewfunction.hh:143
typename GFS::Traits::FiniteElementMap::Traits::FiniteElement::Traits::LocalBasisType::Traits LocalBasisTraits
Definition: discretegridviewfunction.hh:76
std::vector< Range > yb_
Definition: discretegridviewfunction.hh:306
A discrete function defined over a GridFunctionSpace.
Definition: discretegridviewfunction.hh:69
Functions::GridViewEntitySet< GridView, 0 > EntitySet
Definition: discretegridviewfunction.hh:73
GFS Basis
Definition: discretegridviewfunction.hh:90
friend DiscreteGridViewFunction< GFS, V, diffOrder+1 > derivative(const DiscreteGridViewFunction &t)
Definition: discretegridviewfunction.hh:339