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
20namespace 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