3#include <amdis/Output.hpp>
4#include <amdis/algorithm/ForEach.hpp>
5#include <amdis/common/DerivativeTraits.hpp>
6#include <amdis/common/FieldMatVec.hpp>
7#include <amdis/common/Index.hpp>
8#include <amdis/common/Math.hpp>
9#include <amdis/functions/HierarchicNodeToRangeMap.hpp>
10#include <amdis/functions/NodeIndices.hpp>
11#include <amdis/functions/Order.hpp>
12#include <amdis/typetree/FiniteElementType.hpp>
13#include <amdis/utility/LocalToGlobalAdapter.hpp>
15#include <dune/common/ftraits.hh>
21template <
class Coeff,
class LV,
class LC>
22auto gather(Coeff
const& coeff, LV
const& localView, LC& localCoeff, Dune::PriorityTag<2>)
23 -> std::void_t<
decltype(coeff.gather(localView, localCoeff))>
25 coeff.gather(localView, localCoeff);
29template <
class Coeff,
class LV,
class LC,
31auto gather(Coeff
const& coeff, LV
const& localView, LC& localCoeff, Dune::PriorityTag<1>)
32 -> std::void_t<decltype(coeff[std::declval<MI>()])>
34 localCoeff.resize(localView.size());
35 auto it = localCoeff.begin();
36 for (
auto const& idx : nodeIndices(localView))
41template <
class Coeff,
class LV,
class LC>
42auto gather(Coeff
const& coeff, LV
const& localView, LC& localCoeff, Dune::PriorityTag<0>)
43 -> std::void_t<
decltype(coeff[index_<0>])>
45 auto backend = Dune::Functions::istlVectorBackend(coeff);
46 localCoeff.resize(localView.size());
47 auto it = localCoeff.begin();
48 for (
auto const& idx : nodeIndices(localView))
55template <
class C,
class GB,
class TP,
class R>
59 using DomainType =
typename DiscreteFunction::Domain;
60 using RangeType =
typename DiscreteFunction::Range;
63 using DerivativeRange =
typename DerivativeTraits<RangeType(DomainType), T>::Range;
66 using Domain =
typename EntitySet::LocalCoordinate;
67 using Range = DerivativeRange<Type>;
69 enum { hasDerivative = std::is_same<Type, tag::value>::value };
73 using Element =
typename EntitySet::Element;
74 using Geometry =
typename Element::Geometry;
79 REQUIRES(Concepts::Similar<LV,LocalView>)>
81 : localView_(FWD(localView))
88 void bind(Element
const& element)
90 localView_.bind(element);
91 geometry_.emplace(element.geometry());
92 Impl::gather(*coefficients_, localView_, localCoefficients_, Dune::PriorityTag<4>{});
109 Range operator() (
const Domain& local)
const
112 if constexpr (std::is_same_v<Type, tag::value>)
113 return evaluateFunction(local);
114 else if constexpr (std::is_same_v<Type, tag::jacobian>)
115 return evaluateJacobian(local);
116 else if constexpr (std::is_same_v<Type, tag::partial>)
117 return evaluatePartialDerivative(local);
118 else if constexpr (std::is_same_v<Type, tag::divergence>) {
119 static constexpr bool isVectorNode
120 = SubTree::isPower &&
int(SubTree::degree()) == GridView::dimensionworld;
121 static_assert(isVectorNode);
122 if constexpr (isVectorNode)
123 return evaluateDivergence(local);
130 template <
class DerType>
133 static_assert(std::is_same_v<Type, tag::value>);
136 df.
bind(localContext());
142 -> decltype(AMDiS::order(std::declval<SubTree>()))
145 return AMDiS::order(subTree());
152 return localView_.element();
157 template <
class LocalCoordinate>
158 auto evaluateFunction (
const LocalCoordinate& local)
const;
160 template <
class LocalCoordinate>
161 auto evaluateJacobian (
const LocalCoordinate& local)
const;
163 template <
class LocalCoordinate>
164 auto evaluateDivergence (
const LocalCoordinate& local)
const;
166 template <
class LocalCoordinate>
167 auto evaluatePartialDerivative (
const LocalCoordinate& local)
const;
170 Geometry
const& geometry()
const
177 SubTree
const& subTree()
const
179 return Dune::TypeTree::child(localView_.tree(), treePath_);
185 class =
decltype(std::declval<LV>().treeCache())>
186 auto const& subTreeCache(Dune::PriorityTag<1>)
const
188 return Dune::TypeTree::child(localView_.treeCache(), treePath_);
193 auto subTreeCache(Dune::PriorityTag<0>)
const
195 return makeNodeCache(subTree());
198 decltype(
auto) subTreeCache()
const
200 return subTreeCache<LocalView>(Dune::PriorityTag<2>());
206 std::shared_ptr<Coefficients const> coefficients_;
209 std::optional<Geometry> geometry_;
210 std::vector<Coefficient> localCoefficients_;
214template <
class C,
class GB,
class TP,
class R>
215 template <
class Type>
216 template <
class LocalCoordinate>
217auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
218::evaluateFunction(LocalCoordinate
const& local)
const
221 Recursive::forEach(y, [](
auto& v) { v = 0; });
223 Traversal::forEachLeafNode(this->subTreeCache(), [&](
auto const& node,
auto const& tp)
225 auto const& shapeFunctionValues = node.localBasisValuesAt(local);
226 std::size_t size = node.finiteElement().size();
229 auto re = Dune::Functions::flatVectorView(hierarchicNodeToRangeMap(tp, y));
231 for (std::size_t i = 0; i < size; ++i) {
233 auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
236 auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
238 std::size_t dimC = c.size();
239 std::size_t dimV = v.size();
240 assert(dimC*dimV == std::size_t(re.size()));
241 for(std::size_t j = 0; j < dimC; ++j) {
243 for(std::size_t k = 0; k < dimV; ++k)
244 re[j*dimV + k] += c_j*v[k];
253template <
class C,
class GB,
class TP,
class R>
254 template <
class Type>
255 template <
class LocalCoordinate>
256auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
257::evaluateJacobian(LocalCoordinate
const& local)
const
260 Recursive::forEach(dy, [](
auto& v) { v = 0; });
262 Traversal::forEachLeafNode(this->subTreeCache(), [&](
auto const& node,
auto const& tp)
264 LocalToGlobalBasisAdapter localBasis(node, this->geometry());
265 auto const& gradients = localBasis.gradientsAt(local);
268 auto re = Dune::Functions::flatVectorView(hierarchicNodeToRangeMap(tp, dy));
270 for (std::size_t i = 0; i < localBasis.size(); ++i) {
272 auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
275 auto grad = Dune::Functions::flatVectorView(gradients[i]);
277 std::size_t dimC = c.size();
278 std::size_t dimV = grad.size();
279 assert(dimC*dimV == std::size_t(re.size()));
280 for(std::size_t j = 0; j < dimC; ++j) {
282 for(std::size_t k = 0; k < dimV; ++k)
283 re[j*dimV + k] += c_j*grad[k];
292template <
class C,
class GB,
class TP,
class R>
293 template <
class Type>
294 template <
class LocalCoordinate>
295auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
296::evaluateDivergence(LocalCoordinate
const& local)
const
298 static_assert(static_size_v<Range> == 1,
"Divergence implemented for scalar coefficients only.");
301 Recursive::forEach(dy, [](
auto& v) { v = 0; });
303 auto&& node = this->subTreeCache();
305 LocalToGlobalBasisAdapter localBasis(node.child(0), this->geometry());
306 auto const& gradients = localBasis.gradientsAt(local);
308 auto re = Dune::Functions::flatVectorView(dy);
309 assert(
int(re.size()) == 1);
310 for (std::size_t i = 0; i < localBasis.size(); ++i) {
311 auto grad = Dune::Functions::flatVectorView(gradients[i]);
313 assert(
int(grad.size()) == GridView::dimensionworld);
314 for (std::size_t j = 0; j < GridView::dimensionworld; ++j)
315 re[0] += localCoefficients_[node.child(j).localIndex(i)] * grad[j];
322template <
class C,
class GB,
class TP,
class R>
323 template <
class Type>
324 template <
class LocalCoordinate>
325auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
326::evaluatePartialDerivative(LocalCoordinate
const& local)
const
328 std::size_t comp = this->type_.comp;
331 Recursive::forEach(dy, [](
auto& v) { v = 0; });
333 Traversal::forEachLeafNode(this->subTreeCache(), [&](
auto const& node,
auto const& tp)
335 LocalToGlobalBasisAdapter localBasis(node, this->geometry());
336 auto const& partial = localBasis.partialsAt(local, comp);
339 auto re = Dune::Functions::flatVectorView(hierarchicNodeToRangeMap(tp, dy));
341 for (std::size_t i = 0; i < localBasis.size(); ++i) {
343 auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
346 auto d_comp = Dune::Functions::flatVectorView(partial[i]);
348 std::size_t dimC = c.size();
349 std::size_t dimV = d_comp.size();
350 assert(dimC*dimV == std::size_t(re.size()));
351 for(std::size_t j = 0; j < dimC; ++j) {
353 for(std::size_t k = 0; k < dimV; ++k)
354 re[j*dimV + k] += c_j*d_comp[k];
Definition: DiscreteLocalFunction.inc.hpp:58
bool bound() const
Check whether this localfunction is bound to an element.
Definition: DiscreteLocalFunction.inc.hpp:103
auto order() const -> decltype(AMDiS::order(std::declval< SubTree >()))
The polynomial degree of the LocalFunctions basis-functions.
Definition: DiscreteLocalFunction.inc.hpp:141
LocalFunction(LV &&localView, TP const &treePath, std::shared_ptr< C const > coefficients, Type type)
Constructor. Stores a copy of the DiscreteFunction.
Definition: DiscreteLocalFunction.inc.hpp:80
void bind(Element const &element)
Bind the LocalView to the element.
Definition: DiscreteLocalFunction.inc.hpp:88
Element const & localContext() const
Return the bound element.
Definition: DiscreteLocalFunction.inc.hpp:149
void unbind()
Unbind the LocalView from the element.
Definition: DiscreteLocalFunction.inc.hpp:96
LocalFunction< DerType > makeDerivative(DerType type) const
Create a LocalFunction representing the derivative.
Definition: DiscreteLocalFunction.inc.hpp:131
A mutable view on the subspace of a DOFVector,.
Definition: DiscreteFunction.hpp:45
Coefficients & coefficients()
Return the mutable DOFVector.
Definition: DiscreteFunction.hpp:113
AMDiS::LocalView< Self > LocalView
Type of the local view on the restriction of the basis to a single element.
Definition: GlobalBasis.hpp:61
constexpr bool MultiIndex
A multi-index type.
Definition: Concepts.hpp:149
constexpr bool LocalView
A Dune::Functions::LocalView type.
Definition: Concepts.hpp:181
Definition: DerivativeTraits.hpp:29