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