dune-functions  2.8.0
globalvaluedlocalfiniteelement.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_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
5 
6 #include <array>
7 #include <numeric>
8 
9 #include <dune/common/deprecated.hh>
10 #include <dune/common/fmatrix.hh>
11 #include <dune/common/fvector.hh>
12 #include <dune/common/math.hh>
13 
14 #include <dune/geometry/referenceelements.hh>
15 
16 #include <dune/localfunctions/common/localbasis.hh>
17 #include <dune/localfunctions/common/localfiniteelementtraits.hh>
18 #include <dune/localfunctions/common/localinterpolation.hh>
19 
20 namespace Dune::Functions::Impl
21 {
22 
36  struct ContravariantPiolaTransformator
37  {
42  template<typename Values, typename LocalCoordinate, typename Geometry>
43  static auto apply(Values& values,
44  const LocalCoordinate& xi,
45  const Geometry& geometry)
46  {
47  auto jacobianTransposed = geometry.jacobianTransposed(xi);
48  auto integrationElement = geometry.integrationElement(xi);
49 
50  for (auto& value : values)
51  {
52  auto tmp = value;
53  jacobianTransposed.mtv(tmp, value);
54  value /= integrationElement;
55  }
56  }
57 
67  template<typename Gradients, typename LocalCoordinate, typename Geometry>
68  static auto applyJacobian(Gradients& gradients,
69  const LocalCoordinate& xi,
70  const Geometry& geometry)
71  {
72  auto jacobianTransposed = geometry.jacobianTransposed(xi);
73  auto integrationElement = geometry.integrationElement(xi);
74 
75  for (auto& gradient : gradients)
76  {
77  auto tmp = gradient;
78 
79  for (size_t j=0; j<gradient.N(); j++)
80  for (size_t k=0; k<gradient.M(); k++)
81  {
82  gradient[j][k] = 0;
83  for (size_t l=0; l<tmp.N(); l++)
84  gradient[j][k] += jacobianTransposed[l][j] * tmp[l][k];
85  gradient[j][k] /= integrationElement;
86  }
87 
88  }
89  }
90 
98  template<class Function, class LocalCoordinate, class Element>
99  class LocalValuedFunction
100  {
101  const Function& f_;
102  const Element& element_;
103 
104  public:
105 
106  LocalValuedFunction(const Function& f, const Element& element)
107  : f_(f), element_(element)
108  {}
109 
110  auto operator()(const LocalCoordinate& xi) const
111  {
112  auto&& f = Dune::Impl::makeFunctionWithCallOperator<LocalCoordinate>(f_);
113  auto globalValue = f(xi);
114 
115  // Apply the inverse Piola transform
116  auto jacobianInverseTransposed = element_.geometry().jacobianInverseTransposed(xi);
117  auto integrationElement = element_.geometry().integrationElement(xi);
118 
119  auto localValue = globalValue;
120  jacobianInverseTransposed.mtv(globalValue, localValue);
121  localValue *= integrationElement;
122 
123  return localValue;
124  }
125  };
126  };
127 
141  struct CovariantPiolaTransformator
142  {
147  template<typename Values, typename LocalCoordinate, typename Geometry>
148  static auto apply(Values& values,
149  const LocalCoordinate& xi,
150  const Geometry& geometry)
151  {
152  auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
153 
154  for (auto& value : values)
155  {
156  auto tmp = value;
157  jacobianInverseTransposed.mv(tmp, value);
158  }
159  }
160 
170  template<typename Gradients, typename LocalCoordinate, typename Geometry>
171  static auto applyJacobian(Gradients& gradients,
172  const LocalCoordinate& xi,
173  const Geometry& geometry)
174  {
175  auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
176 
177  for (auto& gradient : gradients)
178  {
179  auto tmp = gradient;
180 
181  for (size_t j=0; j<gradient.N(); j++)
182  for (size_t k=0; k<gradient.M(); k++)
183  {
184  gradient[j][k] = 0;
185  for (size_t l=0; l<tmp.N(); l++)
186  gradient[j][k] += jacobianInverseTransposed[j][l] * tmp[l][k];
187  }
188 
189  }
190  }
191 
199  template<class Function, class LocalCoordinate, class Element>
200  class LocalValuedFunction
201  {
202  const Function& f_;
203  const Element& element_;
204 
205  public:
206 
207  LocalValuedFunction(const Function& f, const Element& element)
208  : f_(f), element_(element)
209  {}
210 
211  auto operator()(const LocalCoordinate& xi) const
212  {
213  auto&& f = Dune::Impl::makeFunctionWithCallOperator<LocalCoordinate>(f_);
214  auto globalValue = f(xi);
215 
216  // Apply the inverse Piola transform
217  auto jacobianTransposed = element_.geometry().jacobianTransposed(xi);
218 
219  auto localValue = globalValue;
220  jacobianTransposed.mv(globalValue, localValue);
221 
222  return localValue;
223  }
224  };
225  };
226 
233  template<class Transformator, class LocalValuedLocalBasis, class Element>
234  class GlobalValuedLocalBasis
235  {
236  public:
237  using Traits = typename LocalValuedLocalBasis::Traits;
238 
241  void bind(const LocalValuedLocalBasis& localValuedLocalBasis, const Element& element)
242  {
243  localValuedLocalBasis_ = &localValuedLocalBasis;
244  element_ = &element;
245  }
246 
249  auto size() const
250  {
251  return localValuedLocalBasis_->size();
252  }
253 
255  void evaluateFunction(const typename Traits::DomainType& x,
256  std::vector<typename Traits::RangeType>& out) const
257  {
258  localValuedLocalBasis_->evaluateFunction(x,out);
259 
260  Transformator::apply(out, x, element_->geometry());
261  }
262 
268  void evaluateJacobian(const typename Traits::DomainType& x,
269  std::vector<typename Traits::JacobianType>& out) const
270  {
271  localValuedLocalBasis_->evaluateJacobian(x,out);
272 
273  Transformator::applyJacobian(out, x, element_->geometry());
274  }
275 
282  void partial(const std::array<unsigned int,2>& order,
283  const typename Traits::DomainType& x,
284  std::vector<typename Traits::RangeType>& out) const
285  {
286  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
287  if (totalOrder == 0) {
288  evaluateFunction(x, out);
289  } else if (totalOrder == 1) {
290  auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
291  out.resize(size());
292 
293  // TODO: The following is wasteful: We compute the full Jacobian and then return
294  // only a part of it. While we need the full Jacobian of the underlying local-valued LFE,
295  // it should be possible to compute only a partial Piola transform for the requested
296  // partial derivatives.
297  std::vector<typename Traits::JacobianType> fullJacobian;
298  localValuedLocalBasis_->evaluateJacobian(x,fullJacobian);
299 
300  Transformator::applyJacobian(fullJacobian, x, element_->geometry());
301 
302  for (std::size_t i=0; i<out.size(); i++)
303  for (std::size_t j=0; j<out[i].size(); j++)
304  out[i][j] = fullJacobian[i][j][direction];
305 
306  } else
307  DUNE_THROW(NotImplemented, "Partial derivatives of order 2 or higher");
308  }
309 
311  auto order() const
312  {
313  return localValuedLocalBasis_->order();
314  }
315 
316  const LocalValuedLocalBasis* localValuedLocalBasis_;
317  const Element* element_;
318  };
319 
328  template<class Transformator, class LocalValuedLocalInterpolation, class Element>
329  class GlobalValuedLocalInterpolation
330  {
331  public:
334  void bind(const LocalValuedLocalInterpolation& localValuedLocalInterpolation, const Element& element)
335  {
336  localValuedLocalInterpolation_ = &localValuedLocalInterpolation;
337  element_ = &element;
338  }
339 
340  template<typename F, typename C>
341  void interpolate (const F& f, std::vector<C>& out) const
342  {
343  using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
344  typename Transformator::template LocalValuedFunction<F,LocalCoordinate,Element> localValuedFunction(f, *element_);
345  localValuedLocalInterpolation_->interpolate(localValuedFunction, out);
346  }
347 
348  private:
349  const LocalValuedLocalInterpolation* localValuedLocalInterpolation_;
350  const Element* element_;
351  };
352 
353 
360  template<class Transformator, class LocalValuedLFE, class Element>
361  class GlobalValuedLocalFiniteElement
362  {
363  using LocalBasis = GlobalValuedLocalBasis<Transformator,
364  typename LocalValuedLFE::Traits::LocalBasisType,
365  Element>;
366  using LocalInterpolation = GlobalValuedLocalInterpolation<Transformator,
367  typename LocalValuedLFE::Traits::LocalInterpolationType,
368  Element>;
369 
370  public:
373  using Traits = LocalFiniteElementTraits<LocalBasis,
374  typename LocalValuedLFE::Traits::LocalCoefficientsType,
375  LocalInterpolation>;
376 
377  GlobalValuedLocalFiniteElement() {}
378 
379  void bind(const LocalValuedLFE& localValuedLFE, const Element& element)
380  {
381  globalValuedLocalBasis_.bind(localValuedLFE.localBasis(), element);
382  globalValuedLocalInterpolation_.bind(localValuedLFE.localInterpolation(), element);
383  localValuedLFE_ = &localValuedLFE;
384  }
385 
388  const typename Traits::LocalBasisType& localBasis() const
389  {
390  return globalValuedLocalBasis_;
391  }
392 
395  const typename Traits::LocalCoefficientsType& localCoefficients() const
396  {
397  return localValuedLFE_->localCoefficients();
398  }
399 
402  const typename Traits::LocalInterpolationType& localInterpolation() const
403  {
404  return globalValuedLocalInterpolation_;
405  }
406 
408  std::size_t size() const
409  {
410  return localValuedLFE_->size();
411  }
412 
415  GeometryType type() const
416  {
417  return localValuedLFE_->type();
418  }
419 
420  private:
421 
422  typename Traits::LocalBasisType globalValuedLocalBasis_;
423  typename Traits::LocalInterpolationType globalValuedLocalInterpolation_;
424  const LocalValuedLFE* localValuedLFE_;
425  };
426 
427 } // namespace Dune::Functions::Impl
428 
429 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
void interpolate(const B &basis, C &&coeff, const F &f, const BV &bv, const NTRE &nodeToRangeEntry)
Interpolate given function in discrete function space.
Definition: interpolate.hh:202