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