dune-functions  2.9.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/fmatrix.hh>
10 #include <dune/common/fvector.hh>
11 #include <dune/common/math.hh>
12 #include <dune/common/rangeutilities.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  for (auto& gradient : gradients)
75  {
76  auto tmp = gradient;
77  gradient = 0;
78  for (size_t k=0; k<gradient.M(); k++)
79  for (size_t l=0; l<tmp.N(); l++)
80  // Use sparseRange because jacobianTransposed may be a sparse DiagonalMatrix
81  for(auto&& [jacobianTransposed_l_j, j] : sparseRange(jacobianTransposed[l]))
82  gradient[j][k] += jacobianTransposed_l_j * tmp[l][k];
83  gradient /= integrationElement;
84  }
85  }
86 
94  template<class Function, class LocalCoordinate, class Element>
95  class LocalValuedFunction
96  {
97  const Function& f_;
98  const Element& element_;
99 
100  public:
101 
102  LocalValuedFunction(const Function& f, const Element& element)
103  : f_(f), element_(element)
104  {}
105 
106  auto operator()(const LocalCoordinate& xi) const
107  {
108  auto&& f = Dune::Impl::makeFunctionWithCallOperator<LocalCoordinate>(f_);
109  auto globalValue = f(xi);
110 
111  // Apply the inverse Piola transform
112  auto jacobianInverseTransposed = element_.geometry().jacobianInverseTransposed(xi);
113  auto integrationElement = element_.geometry().integrationElement(xi);
114 
115  auto localValue = globalValue;
116  jacobianInverseTransposed.mtv(globalValue, localValue);
117  localValue *= integrationElement;
118 
119  return localValue;
120  }
121  };
122  };
123 
137  struct CovariantPiolaTransformator
138  {
143  template<typename Values, typename LocalCoordinate, typename Geometry>
144  static auto apply(Values& values,
145  const LocalCoordinate& xi,
146  const Geometry& geometry)
147  {
148  auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
149 
150  for (auto& value : values)
151  {
152  auto tmp = value;
153  jacobianInverseTransposed.mv(tmp, value);
154  }
155  }
156 
166  template<typename Gradients, typename LocalCoordinate, typename Geometry>
167  static auto applyJacobian(Gradients& gradients,
168  const LocalCoordinate& xi,
169  const Geometry& geometry)
170  {
171  auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
172 
173  for (auto& gradient : gradients)
174  {
175  auto tmp = gradient;
176  gradient = 0;
177  for (size_t j=0; j<gradient.N(); j++)
178  for (size_t k=0; k<gradient.M(); k++)
179  // Use sparseRange because jacobianTransposed may be a sparse DiagonalMatrix
180  for(auto&& [jacobianInverseTransposed_j_l, l] : sparseRange(jacobianInverseTransposed[j]))
181  gradient[j][k] += jacobianInverseTransposed_j_l * tmp[l][k];
182  }
183  }
184 
192  template<class Function, class LocalCoordinate, class Element>
193  class LocalValuedFunction
194  {
195  const Function& f_;
196  const Element& element_;
197 
198  public:
199 
200  LocalValuedFunction(const Function& f, const Element& element)
201  : f_(f), element_(element)
202  {}
203 
204  auto operator()(const LocalCoordinate& xi) const
205  {
206  auto&& f = Dune::Impl::makeFunctionWithCallOperator<LocalCoordinate>(f_);
207  auto globalValue = f(xi);
208 
209  // Apply the inverse Piola transform
210  auto jacobianTransposed = element_.geometry().jacobianTransposed(xi);
211 
212  auto localValue = globalValue;
213  jacobianTransposed.mv(globalValue, localValue);
214 
215  return localValue;
216  }
217  };
218  };
219 
226  template<class Transformator, class LocalValuedLocalBasis, class Element>
227  class GlobalValuedLocalBasis
228  {
229  public:
230  using Traits = typename LocalValuedLocalBasis::Traits;
231 
234  void bind(const LocalValuedLocalBasis& localValuedLocalBasis, const Element& element)
235  {
236  localValuedLocalBasis_ = &localValuedLocalBasis;
237  element_ = &element;
238  }
239 
242  auto size() const
243  {
244  return localValuedLocalBasis_->size();
245  }
246 
248  void evaluateFunction(const typename Traits::DomainType& x,
249  std::vector<typename Traits::RangeType>& out) const
250  {
251  localValuedLocalBasis_->evaluateFunction(x,out);
252 
253  Transformator::apply(out, x, element_->geometry());
254  }
255 
261  void evaluateJacobian(const typename Traits::DomainType& x,
262  std::vector<typename Traits::JacobianType>& out) const
263  {
264  localValuedLocalBasis_->evaluateJacobian(x,out);
265 
266  Transformator::applyJacobian(out, x, element_->geometry());
267  }
268 
275  void partial(const std::array<unsigned int,2>& order,
276  const typename Traits::DomainType& x,
277  std::vector<typename Traits::RangeType>& out) const
278  {
279  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
280  if (totalOrder == 0) {
281  evaluateFunction(x, out);
282  } else if (totalOrder == 1) {
283  auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
284  out.resize(size());
285 
286  // TODO: The following is wasteful: We compute the full Jacobian and then return
287  // only a part of it. While we need the full Jacobian of the underlying local-valued LFE,
288  // it should be possible to compute only a partial Piola transform for the requested
289  // partial derivatives.
290  std::vector<typename Traits::JacobianType> fullJacobian;
291  localValuedLocalBasis_->evaluateJacobian(x,fullJacobian);
292 
293  Transformator::applyJacobian(fullJacobian, x, element_->geometry());
294 
295  for (std::size_t i=0; i<out.size(); i++)
296  for (std::size_t j=0; j<out[i].size(); j++)
297  out[i][j] = fullJacobian[i][j][direction];
298 
299  } else
300  DUNE_THROW(NotImplemented, "Partial derivatives of order 2 or higher");
301  }
302 
304  auto order() const
305  {
306  return localValuedLocalBasis_->order();
307  }
308 
309  const LocalValuedLocalBasis* localValuedLocalBasis_;
310  const Element* element_;
311  };
312 
321  template<class Transformator, class LocalValuedLocalInterpolation, class Element>
322  class GlobalValuedLocalInterpolation
323  {
324  public:
327  void bind(const LocalValuedLocalInterpolation& localValuedLocalInterpolation, const Element& element)
328  {
329  localValuedLocalInterpolation_ = &localValuedLocalInterpolation;
330  element_ = &element;
331  }
332 
333  template<typename F, typename C>
334  void interpolate (const F& f, std::vector<C>& out) const
335  {
336  using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
337  typename Transformator::template LocalValuedFunction<F,LocalCoordinate,Element> localValuedFunction(f, *element_);
338  localValuedLocalInterpolation_->interpolate(localValuedFunction, out);
339  }
340 
341  private:
342  const LocalValuedLocalInterpolation* localValuedLocalInterpolation_;
343  const Element* element_;
344  };
345 
346 
353  template<class Transformator, class LocalValuedLFE, class Element>
354  class GlobalValuedLocalFiniteElement
355  {
356  using LocalBasis = GlobalValuedLocalBasis<Transformator,
357  typename LocalValuedLFE::Traits::LocalBasisType,
358  Element>;
359  using LocalInterpolation = GlobalValuedLocalInterpolation<Transformator,
360  typename LocalValuedLFE::Traits::LocalInterpolationType,
361  Element>;
362 
363  public:
366  using Traits = LocalFiniteElementTraits<LocalBasis,
367  typename LocalValuedLFE::Traits::LocalCoefficientsType,
368  LocalInterpolation>;
369 
370  GlobalValuedLocalFiniteElement() {}
371 
372  void bind(const LocalValuedLFE& localValuedLFE, const Element& element)
373  {
374  globalValuedLocalBasis_.bind(localValuedLFE.localBasis(), element);
375  globalValuedLocalInterpolation_.bind(localValuedLFE.localInterpolation(), element);
376  localValuedLFE_ = &localValuedLFE;
377  }
378 
381  const typename Traits::LocalBasisType& localBasis() const
382  {
383  return globalValuedLocalBasis_;
384  }
385 
388  const typename Traits::LocalCoefficientsType& localCoefficients() const
389  {
390  return localValuedLFE_->localCoefficients();
391  }
392 
395  const typename Traits::LocalInterpolationType& localInterpolation() const
396  {
397  return globalValuedLocalInterpolation_;
398  }
399 
401  std::size_t size() const
402  {
403  return localValuedLFE_->size();
404  }
405 
408  GeometryType type() const
409  {
410  return localValuedLFE_->type();
411  }
412 
413  private:
414 
415  typename Traits::LocalBasisType globalValuedLocalBasis_;
416  typename Traits::LocalInterpolationType globalValuedLocalInterpolation_;
417  const LocalValuedLFE* localValuedLFE_;
418  };
419 
420 } // namespace Dune::Functions::Impl
421 
422 #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