3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
14#include <dune/geometry/referenceelements.hh>
16#include <dune/localfunctions/common/localbasis.hh>
17#include <dune/localfunctions/common/localfiniteelementtraits.hh>
19namespace Dune::Functions::Impl
35 struct ContravariantPiolaTransformator
41 template<
typename Values,
typename LocalCoordinate,
typename Geometry>
42 static auto apply(Values& values,
43 const LocalCoordinate& xi,
44 const Geometry& geometry)
46 auto jacobianTransposed = geometry.jacobianTransposed(xi);
47 auto integrationElement = geometry.integrationElement(xi);
49 for (
auto& value : values)
52 jacobianTransposed.mtv(tmp, value);
53 value /= integrationElement;
66 template<
typename Gradients,
typename LocalCoordinate,
typename Geometry>
67 static auto applyJacobian(Gradients& gradients,
68 const LocalCoordinate& xi,
69 const Geometry& geometry)
71 auto jacobianTransposed = geometry.jacobianTransposed(xi);
72 auto integrationElement = geometry.integrationElement(xi);
73 for (
auto& gradient : gradients)
77 for (
size_t k=0; k<gradient.M(); k++)
78 for (
size_t l=0; l<tmp.N(); l++)
80 for(
auto&& [jacobianTransposed_l_j, j] :
sparseRange(jacobianTransposed[l]))
81 gradient[j][k] += jacobianTransposed_l_j * tmp[l][k];
82 gradient /= integrationElement;
93 template<
class Function,
class LocalCoordinate,
class Element>
94 class LocalValuedFunction
97 const Element& element_;
99 using LocalValue = LocalCoordinate;
103 LocalValuedFunction(
const Function& f,
const Element& element)
104 : f_(f), element_(element)
107 auto operator()(
const LocalCoordinate& xi)
const
109 auto globalValue = f_(xi);
112 auto jacobianInverseTransposed = element_.geometry().jacobianInverseTransposed(xi);
113 auto integrationElement = element_.geometry().integrationElement(xi);
115 auto localValue = LocalValue{};
116 jacobianInverseTransposed.mtv(globalValue, localValue);
117 localValue *= integrationElement;
137 struct CovariantPiolaTransformator
143 template<
typename Values,
typename LocalCoordinate,
typename Geometry>
144 static auto apply(Values& values,
145 const LocalCoordinate& xi,
146 const Geometry& geometry)
148 auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
150 for (
auto& value : values)
153 jacobianInverseTransposed.mv(tmp, value);
166 template<
typename Gradients,
typename LocalCoordinate,
typename Geometry>
167 static auto applyJacobian(Gradients& gradients,
168 const LocalCoordinate& xi,
169 const Geometry& geometry)
171 auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
173 for (
auto& gradient : gradients)
177 for (
size_t j=0; j<gradient.N(); j++)
178 for (
size_t k=0; k<gradient.M(); k++)
180 for(
auto&& [jacobianInverseTransposed_j_l, l] :
sparseRange(jacobianInverseTransposed[j]))
181 gradient[j][k] += jacobianInverseTransposed_j_l * tmp[l][k];
192 template<
class Function,
class LocalCoordinate,
class Element>
193 class LocalValuedFunction
196 const Element& element_;
200 LocalValuedFunction(
const Function& f,
const Element& element)
201 : f_(f), element_(element)
204 auto operator()(
const LocalCoordinate& xi)
const
206 auto globalValue = f_(xi);
209 auto jacobianTransposed = element_.geometry().jacobianTransposed(xi);
211 auto localValue = globalValue;
212 jacobianTransposed.mv(globalValue, localValue);
225 template<
class Transformator,
class LocalValuedLocalBasis,
class Element>
226 class GlobalValuedLocalBasis
229 using Traits =
typename LocalValuedLocalBasis::Traits;
233 void bind(
const LocalValuedLocalBasis& localValuedLocalBasis,
const Element& element)
235 localValuedLocalBasis_ = &localValuedLocalBasis;
243 return localValuedLocalBasis_->size();
247 void evaluateFunction(
const typename Traits::DomainType& x,
248 std::vector<typename Traits::RangeType>& out)
const
250 localValuedLocalBasis_->evaluateFunction(x,out);
252 Transformator::apply(out, x, element_->geometry());
260 void evaluateJacobian(
const typename Traits::DomainType& x,
261 std::vector<typename Traits::JacobianType>& out)
const
263 localValuedLocalBasis_->evaluateJacobian(x,out);
265 Transformator::applyJacobian(out, x, element_->geometry());
274 void partial(
const std::array<unsigned int,2>& order,
275 const typename Traits::DomainType& x,
276 std::vector<typename Traits::RangeType>& out)
const
279 if (totalOrder == 0) {
280 evaluateFunction(x, out);
281 }
else if (totalOrder == 1) {
282 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
289 std::vector<typename Traits::JacobianType> fullJacobian;
290 localValuedLocalBasis_->evaluateJacobian(x,fullJacobian);
292 Transformator::applyJacobian(fullJacobian, x, element_->geometry());
294 for (std::size_t i=0; i<out.size(); i++)
295 for (std::size_t j=0; j<out[i].size(); j++)
296 out[i][j] = fullJacobian[i][j][direction];
299 DUNE_THROW(NotImplemented,
"Partial derivatives of order 2 or higher");
305 return localValuedLocalBasis_->order();
308 const LocalValuedLocalBasis* localValuedLocalBasis_;
309 const Element* element_;
320 template<
class Transformator,
class LocalValuedLocalInterpolation,
class Element>
321 class GlobalValuedLocalInterpolation
326 void bind(
const LocalValuedLocalInterpolation& localValuedLocalInterpolation,
const Element& element)
328 localValuedLocalInterpolation_ = &localValuedLocalInterpolation;
332 template<
typename F,
typename C>
333 void interpolate (
const F& f, std::vector<C>& out)
const
335 using LocalCoordinate =
typename Element::Geometry::LocalCoordinate;
336 typename Transformator::template LocalValuedFunction<F,LocalCoordinate,Element> localValuedFunction(f, *element_);
337 localValuedLocalInterpolation_->interpolate(localValuedFunction, out);
341 const LocalValuedLocalInterpolation* localValuedLocalInterpolation_;
342 const Element* element_;
352 template<
class Transformator,
class LocalValuedLFE,
class Element>
353 class GlobalValuedLocalFiniteElement
355 using LocalBasis = GlobalValuedLocalBasis<Transformator,
356 typename LocalValuedLFE::Traits::LocalBasisType,
358 using LocalInterpolation = GlobalValuedLocalInterpolation<Transformator,
359 typename LocalValuedLFE::Traits::LocalInterpolationType,
365 using Traits = LocalFiniteElementTraits<LocalBasis,
366 typename LocalValuedLFE::Traits::LocalCoefficientsType,
369 GlobalValuedLocalFiniteElement() {}
371 void bind(
const LocalValuedLFE& localValuedLFE,
const Element& element)
373 globalValuedLocalBasis_.bind(localValuedLFE.localBasis(), element);
374 globalValuedLocalInterpolation_.bind(localValuedLFE.localInterpolation(), element);
375 localValuedLFE_ = &localValuedLFE;
382 return globalValuedLocalBasis_;
389 return localValuedLFE_->localCoefficients();
396 return globalValuedLocalInterpolation_;
400 std::size_t
size()
const
402 return localValuedLFE_->size();
409 return localValuedLFE_->type();
416 const LocalValuedLFE* localValuedLFE_;
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
auto sparseRange(Range &&range)
Allow structured-binding for-loops for sparse iterators.
Definition: rangeutilities.hh:722
Some useful basic math stuff.
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Utilities for reduction like operations on ranges.
LB LocalBasisType
Definition: localfiniteelementtraits.hh:16
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24