7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
13#include <dune/common/fmatrix.hh>
14#include <dune/common/fvector.hh>
15#include <dune/common/math.hh>
16#include <dune/common/rangeutilities.hh>
18#include <dune/geometry/referenceelements.hh>
20#include <dune/localfunctions/common/localbasis.hh>
21#include <dune/localfunctions/common/localfiniteelementtraits.hh>
23namespace Dune::Functions::Impl
39 struct ContravariantPiolaTransformator
45 template<
typename Values,
typename LocalCoordinate,
typename Geometry>
46 static auto apply(Values& values,
47 const LocalCoordinate& xi,
48 const Geometry& geometry)
50 auto jacobianTransposed = geometry.jacobianTransposed(xi);
51 auto integrationElement = geometry.integrationElement(xi);
53 for (
auto& value : values)
56 jacobianTransposed.mtv(tmp, value);
57 value /= integrationElement;
70 template<
typename Gradients,
typename LocalCoordinate,
typename Geometry>
71 static auto applyJacobian(Gradients& gradients,
72 const LocalCoordinate& xi,
73 const Geometry& geometry)
75 auto jacobianTransposed = geometry.jacobianTransposed(xi);
76 auto integrationElement = geometry.integrationElement(xi);
77 for (
auto& gradient : gradients)
81 for (
size_t k=0; k<gradient.M(); k++)
82 for (
size_t l=0; l<tmp.N(); l++)
84 for(
auto&& [jacobianTransposed_l_j, j] : sparseRange(jacobianTransposed[l]))
85 gradient[j][k] += jacobianTransposed_l_j * tmp[l][k];
86 gradient /= integrationElement;
97 template<
class Function,
class LocalCoordinate,
class Element>
98 class LocalValuedFunction
101 const Element& element_;
103 using LocalValue = LocalCoordinate;
107 LocalValuedFunction(
const Function& f,
const Element& element)
108 : f_(f), element_(element)
111 auto operator()(
const LocalCoordinate& xi)
const
113 auto globalValue = f_(xi);
116 auto jacobianInverseTransposed = element_.geometry().jacobianInverseTransposed(xi);
117 auto integrationElement = element_.geometry().integrationElement(xi);
119 auto localValue = LocalValue{};
120 jacobianInverseTransposed.mtv(globalValue, localValue);
121 localValue *= integrationElement;
141 struct CovariantPiolaTransformator
147 template<
typename Values,
typename LocalCoordinate,
typename Geometry>
148 static auto apply(Values& values,
149 const LocalCoordinate& xi,
150 const Geometry& geometry)
152 auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
154 for (
auto& value : values)
157 jacobianInverseTransposed.mv(tmp, value);
170 template<
typename Gradients,
typename LocalCoordinate,
typename Geometry>
171 static auto applyJacobian(Gradients& gradients,
172 const LocalCoordinate& xi,
173 const Geometry& geometry)
175 auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
177 for (
auto& gradient : gradients)
181 for (
size_t j=0; j<gradient.N(); j++)
182 for (
size_t k=0; k<gradient.M(); k++)
184 for(
auto&& [jacobianInverseTransposed_j_l, l] : sparseRange(jacobianInverseTransposed[j]))
185 gradient[j][k] += jacobianInverseTransposed_j_l * tmp[l][k];
196 template<
class Function,
class LocalCoordinate,
class Element>
197 class LocalValuedFunction
200 const Element& element_;
204 LocalValuedFunction(
const Function& f,
const Element& element)
205 : f_(f), element_(element)
208 auto operator()(
const LocalCoordinate& xi)
const
210 auto globalValue = f_(xi);
213 auto jacobianTransposed = element_.geometry().jacobianTransposed(xi);
215 auto localValue = globalValue;
216 jacobianTransposed.mv(globalValue, localValue);
229 template<
class Transformator,
class LocalValuedLocalBasis,
class Element>
230 class GlobalValuedLocalBasis
233 using Traits =
typename LocalValuedLocalBasis::Traits;
237 void bind(
const LocalValuedLocalBasis& localValuedLocalBasis,
const Element& element)
239 localValuedLocalBasis_ = &localValuedLocalBasis;
247 return localValuedLocalBasis_->size();
251 void evaluateFunction(
const typename Traits::DomainType& x,
252 std::vector<typename Traits::RangeType>& out)
const
254 localValuedLocalBasis_->evaluateFunction(x,out);
256 Transformator::apply(out, x, element_->geometry());
264 void evaluateJacobian(
const typename Traits::DomainType& x,
265 std::vector<typename Traits::JacobianType>& out)
const
267 localValuedLocalBasis_->evaluateJacobian(x,out);
269 Transformator::applyJacobian(out, x, element_->geometry());
278 void partial(
const std::array<unsigned int,2>& order,
279 const typename Traits::DomainType& x,
280 std::vector<typename Traits::RangeType>& out)
const
282 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
283 if (totalOrder == 0) {
284 evaluateFunction(x, out);
285 }
else if (totalOrder == 1) {
286 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
293 std::vector<typename Traits::JacobianType> fullJacobian;
294 localValuedLocalBasis_->evaluateJacobian(x,fullJacobian);
296 Transformator::applyJacobian(fullJacobian, x, element_->geometry());
298 for (std::size_t i=0; i<out.size(); i++)
299 for (std::size_t j=0; j<out[i].size(); j++)
300 out[i][j] = fullJacobian[i][j][direction];
303 DUNE_THROW(NotImplemented,
"Partial derivatives of order 2 or higher");
309 return localValuedLocalBasis_->order();
312 const LocalValuedLocalBasis* localValuedLocalBasis_;
313 const Element* element_;
324 template<
class Transformator,
class LocalValuedLocalInterpolation,
class Element>
325 class GlobalValuedLocalInterpolation
330 void bind(
const LocalValuedLocalInterpolation& localValuedLocalInterpolation,
const Element& element)
332 localValuedLocalInterpolation_ = &localValuedLocalInterpolation;
336 template<
typename F,
typename C>
337 void interpolate (
const F& f, std::vector<C>& out)
const
339 using LocalCoordinate =
typename Element::Geometry::LocalCoordinate;
340 typename Transformator::template LocalValuedFunction<F,LocalCoordinate,Element> localValuedFunction(f, *element_);
341 localValuedLocalInterpolation_->interpolate(localValuedFunction, out);
345 const LocalValuedLocalInterpolation* localValuedLocalInterpolation_;
346 const Element* element_;
356 template<
class Transformator,
class LocalValuedLFE,
class Element>
357 class GlobalValuedLocalFiniteElement
359 using LocalBasis = GlobalValuedLocalBasis<Transformator,
360 typename LocalValuedLFE::Traits::LocalBasisType,
362 using LocalInterpolation = GlobalValuedLocalInterpolation<Transformator,
363 typename LocalValuedLFE::Traits::LocalInterpolationType,
369 using Traits = LocalFiniteElementTraits<LocalBasis,
370 typename LocalValuedLFE::Traits::LocalCoefficientsType,
373 GlobalValuedLocalFiniteElement() {}
375 void bind(
const LocalValuedLFE& localValuedLFE,
const Element& element)
377 globalValuedLocalBasis_.bind(localValuedLFE.localBasis(), element);
378 globalValuedLocalInterpolation_.bind(localValuedLFE.localInterpolation(), element);
379 localValuedLFE_ = &localValuedLFE;
384 const typename Traits::LocalBasisType& localBasis()
const
386 return globalValuedLocalBasis_;
391 const typename Traits::LocalCoefficientsType& localCoefficients()
const
393 return localValuedLFE_->localCoefficients();
398 const typename Traits::LocalInterpolationType& localInterpolation()
const
400 return globalValuedLocalInterpolation_;
404 std::size_t size()
const
406 return localValuedLFE_->size();
411 GeometryType type()
const
413 return localValuedLFE_->type();
418 typename Traits::LocalBasisType globalValuedLocalBasis_;
419 typename Traits::LocalInterpolationType globalValuedLocalInterpolation_;
420 const LocalValuedLFE* localValuedLFE_;