DUNE PDELab (git)

globalvaluedlocalfiniteelement.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
5// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
6
7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
9
10#include <array>
11#include <numeric>
12
15#include <dune/common/math.hh>
17
18#include <dune/geometry/referenceelements.hh>
19
20#include <dune/localfunctions/common/localbasis.hh>
21#include <dune/localfunctions/common/localfiniteelementtraits.hh>
22
23namespace Dune::Functions::Impl
24{
25
39 struct ContravariantPiolaTransformator
40 {
45 template<typename Values, typename LocalCoordinate, typename Geometry>
46 static auto apply(Values& values,
47 const LocalCoordinate& xi,
48 const Geometry& geometry)
49 {
50 auto jacobianTransposed = geometry.jacobianTransposed(xi);
51 auto integrationElement = geometry.integrationElement(xi);
52
53 for (auto& value : values)
54 {
55 auto tmp = value;
56 jacobianTransposed.mtv(tmp, value);
57 value /= integrationElement;
58 }
59 }
60
70 template<typename Gradients, typename LocalCoordinate, typename Geometry>
71 static auto applyJacobian(Gradients& gradients,
72 const LocalCoordinate& xi,
73 const Geometry& geometry)
74 {
75 auto jacobianTransposed = geometry.jacobianTransposed(xi);
76 auto integrationElement = geometry.integrationElement(xi);
77 for (auto& gradient : gradients)
78 {
79 auto tmp = gradient;
80 gradient = 0;
81 for (size_t k=0; k<gradient.M(); k++)
82 for (size_t l=0; l<tmp.N(); l++)
83 // Use sparseRange because jacobianTransposed may be a sparse DiagonalMatrix
84 for(auto&& [jacobianTransposed_l_j, j] : sparseRange(jacobianTransposed[l]))
85 gradient[j][k] += jacobianTransposed_l_j * tmp[l][k];
86 gradient /= integrationElement;
87 }
88 }
89
97 template<class Function, class LocalCoordinate, class Element>
98 class LocalValuedFunction
99 {
100 const Function& f_;
101 const Element& element_;
102
103 using LocalValue = LocalCoordinate;
104
105 public:
106
107 LocalValuedFunction(const Function& f, const Element& element)
108 : f_(f), element_(element)
109 {}
110
111 auto operator()(const LocalCoordinate& xi) const
112 {
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 = LocalValue{};
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 gradient = 0;
181 for (size_t j=0; j<gradient.N(); j++)
182 for (size_t k=0; k<gradient.M(); k++)
183 // Use sparseRange because jacobianTransposed may be a sparse DiagonalMatrix
184 for(auto&& [jacobianInverseTransposed_j_l, l] : sparseRange(jacobianInverseTransposed[j]))
185 gradient[j][k] += jacobianInverseTransposed_j_l * tmp[l][k];
186 }
187 }
188
196 template<class Function, class LocalCoordinate, class Element>
197 class LocalValuedFunction
198 {
199 const Function& f_;
200 const Element& element_;
201
202 public:
203
204 LocalValuedFunction(const Function& f, const Element& element)
205 : f_(f), element_(element)
206 {}
207
208 auto operator()(const LocalCoordinate& xi) const
209 {
210 auto globalValue = f_(xi);
211
212 // Apply the inverse Piola transform
213 auto jacobianTransposed = element_.geometry().jacobianTransposed(xi);
214
215 auto localValue = globalValue;
216 jacobianTransposed.mv(globalValue, localValue);
217
218 return localValue;
219 }
220 };
221 };
222
229 template<class Transformator, class LocalValuedLocalBasis, class Element>
230 class GlobalValuedLocalBasis
231 {
232 public:
233 using Traits = typename LocalValuedLocalBasis::Traits;
234
237 void bind(const LocalValuedLocalBasis& localValuedLocalBasis, const Element& element)
238 {
239 localValuedLocalBasis_ = &localValuedLocalBasis;
240 element_ = &element;
241 }
242
245 auto size() const
246 {
247 return localValuedLocalBasis_->size();
248 }
249
251 void evaluateFunction(const typename Traits::DomainType& x,
252 std::vector<typename Traits::RangeType>& out) const
253 {
254 localValuedLocalBasis_->evaluateFunction(x,out);
255
256 Transformator::apply(out, x, element_->geometry());
257 }
258
264 void evaluateJacobian(const typename Traits::DomainType& x,
265 std::vector<typename Traits::JacobianType>& out) const
266 {
267 localValuedLocalBasis_->evaluateJacobian(x,out);
268
269 Transformator::applyJacobian(out, x, element_->geometry());
270 }
271
278 void partial(const std::array<unsigned int,2>& order,
279 const typename Traits::DomainType& x,
280 std::vector<typename Traits::RangeType>& out) const
281 {
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));
287 out.resize(size());
288
289 // TODO: The following is wasteful: We compute the full Jacobian and then return
290 // only a part of it. While we need the full Jacobian of the underlying local-valued LFE,
291 // it should be possible to compute only a partial Piola transform for the requested
292 // partial derivatives.
293 std::vector<typename Traits::JacobianType> fullJacobian;
294 localValuedLocalBasis_->evaluateJacobian(x,fullJacobian);
295
296 Transformator::applyJacobian(fullJacobian, x, element_->geometry());
297
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];
301
302 } else
303 DUNE_THROW(NotImplemented, "Partial derivatives of order 2 or higher");
304 }
305
307 auto order() const
308 {
309 return localValuedLocalBasis_->order();
310 }
311
312 const LocalValuedLocalBasis* localValuedLocalBasis_;
313 const Element* element_;
314 };
315
324 template<class Transformator, class LocalValuedLocalInterpolation, class Element>
325 class GlobalValuedLocalInterpolation
326 {
327 public:
330 void bind(const LocalValuedLocalInterpolation& localValuedLocalInterpolation, const Element& element)
331 {
332 localValuedLocalInterpolation_ = &localValuedLocalInterpolation;
333 element_ = &element;
334 }
335
336 template<typename F, typename C>
337 void interpolate (const F& f, std::vector<C>& out) const
338 {
339 using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
340 typename Transformator::template LocalValuedFunction<F,LocalCoordinate,Element> localValuedFunction(f, *element_);
341 localValuedLocalInterpolation_->interpolate(localValuedFunction, out);
342 }
343
344 private:
345 const LocalValuedLocalInterpolation* localValuedLocalInterpolation_;
346 const Element* element_;
347 };
348
349
356 template<class Transformator, class LocalValuedLFE, class Element>
357 class GlobalValuedLocalFiniteElement
358 {
359 using LocalBasis = GlobalValuedLocalBasis<Transformator,
360 typename LocalValuedLFE::Traits::LocalBasisType,
361 Element>;
362 using LocalInterpolation = GlobalValuedLocalInterpolation<Transformator,
363 typename LocalValuedLFE::Traits::LocalInterpolationType,
364 Element>;
365
366 public:
369 using Traits = LocalFiniteElementTraits<LocalBasis,
370 typename LocalValuedLFE::Traits::LocalCoefficientsType,
371 LocalInterpolation>;
372
373 GlobalValuedLocalFiniteElement() {}
374
375 void bind(const LocalValuedLFE& localValuedLFE, const Element& element)
376 {
377 globalValuedLocalBasis_.bind(localValuedLFE.localBasis(), element);
378 globalValuedLocalInterpolation_.bind(localValuedLFE.localInterpolation(), element);
379 localValuedLFE_ = &localValuedLFE;
380 }
381
384 const typename Traits::LocalBasisType& localBasis() const
385 {
386 return globalValuedLocalBasis_;
387 }
388
391 const typename Traits::LocalCoefficientsType& localCoefficients() const
392 {
393 return localValuedLFE_->localCoefficients();
394 }
395
398 const typename Traits::LocalInterpolationType& localInterpolation() const
399 {
400 return globalValuedLocalInterpolation_;
401 }
402
404 std::size_t size() const
405 {
406 return localValuedLFE_->size();
407 }
408
411 GeometryType type() const
412 {
413 return localValuedLFE_->type();
414 }
415
416 private:
417
418 typename Traits::LocalBasisType globalValuedLocalBasis_;
419 typename Traits::LocalInterpolationType globalValuedLocalInterpolation_;
420 const LocalValuedLFE* localValuedLFE_;
421 };
422
423} // namespace Dune::Functions::Impl
424
425#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)