AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
DiscreteLocalFunction.inc.hpp
1#pragma once
2
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>
14
15#include <dune/common/ftraits.hh>
16
17namespace AMDiS {
18namespace Impl {
19
20// specialization for containers with gather method
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))>
24{
25 coeff.gather(localView, localCoeff);
26}
27
28// implementation for containers with multi-index access
29template <class Coeff, class LV, class LC,
30 class MI = typename LV::MultiIndex>
31auto gather(Coeff const& coeff, LV const& localView, LC& localCoeff, Dune::PriorityTag<1>)
32 -> std::void_t<decltype(coeff[std::declval<MI>()])>
33{
34 localCoeff.resize(localView.size());
35 auto it = localCoeff.begin();
36 for (auto const& idx : nodeIndices(localView))
37 *it++ = coeff[idx];
38}
39
40// fallback implementation
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>])>
44{
45 auto backend = Dune::Functions::istlVectorBackend(coeff);
46 localCoeff.resize(localView.size());
47 auto it = localCoeff.begin();
48 for (auto const& idx : nodeIndices(localView))
49 *it++ = backend[idx];
50}
51
52} // end namespace Impl
53
54
55template <class C, class GB, class TP, class R>
56 template <class Type>
57class DiscreteFunction<C const,GB,TP,R>::LocalFunction
58{
59 using DomainType = typename DiscreteFunction::Domain;
60 using RangeType = typename DiscreteFunction::Range;
61
62 template <class T>
63 using DerivativeRange = typename DerivativeTraits<RangeType(DomainType), T>::Range;
64
65public:
66 using Domain = typename EntitySet::LocalCoordinate;
67 using Range = DerivativeRange<Type>;
68
69 enum { hasDerivative = std::is_same<Type, tag::value>::value };
70
71private:
72 using LocalView = typename GlobalBasis::LocalView;
73 using Element = typename EntitySet::Element;
74 using Geometry = typename Element::Geometry;
75
76public:
78 template <class LV,
79 REQUIRES(Concepts::Similar<LV,LocalView>)>
80 LocalFunction(LV&& localView, TP const& treePath, std::shared_ptr<C const> coefficients, Type type)
81 : localView_(FWD(localView))
82 , treePath_(treePath)
83 , coefficients_(std::move(coefficients))
84 , type_(type)
85 {}
86
88 void bind(Element const& element)
89 {
90 localView_.bind(element);
91 geometry_.emplace(element.geometry());
92 Impl::gather(*coefficients_, localView_, localCoefficients_, Dune::PriorityTag<4>{});
93 }
94
96 void unbind()
97 {
98 localView_.unbind();
99 geometry_.reset();
100 }
101
103 bool bound() const
104 {
105 return !!geometry_;
106 }
107
109 Range operator() (const Domain& local) const
110 {
111 assert(bound());
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);
124 }
125
126 return Range{};
127 }
128
130 template <class DerType>
132 {
133 static_assert(std::is_same_v<Type, tag::value>);
134 LocalFunction<DerType> df {localView_, treePath_, coefficients_, type};
135 if (bound())
136 df.bind(localContext());
137 return df;
138 }
139
141 auto order() const
142 -> decltype(AMDiS::order(std::declval<SubTree>()))
143 {
144 assert(bound());
145 return AMDiS::order(subTree());
146 }
147
149 Element const& localContext() const
150 {
151 assert(bound());
152 return localView_.element();
153 }
154
155private:
156
157 template <class LocalCoordinate>
158 auto evaluateFunction (const LocalCoordinate& local) const;
159
160 template <class LocalCoordinate>
161 auto evaluateJacobian (const LocalCoordinate& local) const;
162
163 template <class LocalCoordinate>
164 auto evaluateDivergence (const LocalCoordinate& local) const;
165
166 template <class LocalCoordinate>
167 auto evaluatePartialDerivative (const LocalCoordinate& local) const;
168
169 // get the geometry of the bound element, that is constructed in the bind() method
170 Geometry const& geometry() const
171 {
172 assert(!!geometry_);
173 return *geometry_;
174 }
175
176 // return the sub-tree of the basis-tree associated to the tree-path
177 SubTree const& subTree() const
178 {
179 return Dune::TypeTree::child(localView_.tree(), treePath_);
180 }
181
182 // return the sub-tree of the basis-tree-cache associated to the tree-path
183 // NOTE: this is only return in case the local-view does provide a treeCache() function
184 template <class LV,
185 class = decltype(std::declval<LV>().treeCache())>
186 auto const& subTreeCache(Dune::PriorityTag<1>) const
187 {
188 return Dune::TypeTree::child(localView_.treeCache(), treePath_);
189 }
190
191 // construct a new tree-cache that stores a reference to the sub-tree
192 template <class>
193 auto subTreeCache(Dune::PriorityTag<0>) const
194 {
195 return makeNodeCache(subTree());
196 }
197
198 decltype(auto) subTreeCache() const
199 {
200 return subTreeCache<LocalView>(Dune::PriorityTag<2>());
201 }
202
203private:
204 LocalView localView_;
205 TP treePath_;
206 std::shared_ptr<Coefficients const> coefficients_;
207 Type type_;
208
209 std::optional<Geometry> geometry_;
210 std::vector<Coefficient> localCoefficients_;
211};
212
213
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
219{
220 Range y;
221 Recursive::forEach(y, [](auto& v) { v = 0; });
222
223 Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
224 {
225 auto const& shapeFunctionValues = node.localBasisValuesAt(local);
226 std::size_t size = node.finiteElement().size();
227
228 // Get range entry associated to this node
229 auto re = Dune::Functions::flatVectorView(hierarchicNodeToRangeMap(tp, y));
230
231 for (std::size_t i = 0; i < size; ++i) {
232 // Get coefficient associated to i-th shape function
233 auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
234
235 // Get value of i-th shape function
236 auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
237
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) {
242 auto&& c_j = c[j];
243 for(std::size_t k = 0; k < dimV; ++k)
244 re[j*dimV + k] += c_j*v[k];
245 }
246 }
247 });
248
249 return y;
250}
251
252
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
258{
259 Range dy;
260 Recursive::forEach(dy, [](auto& v) { v = 0; });
261
262 Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
263 {
264 LocalToGlobalBasisAdapter localBasis(node, this->geometry());
265 auto const& gradients = localBasis.gradientsAt(local);
266
267 // Get range entry associated to this node
268 auto re = Dune::Functions::flatVectorView(hierarchicNodeToRangeMap(tp, dy));
269
270 for (std::size_t i = 0; i < localBasis.size(); ++i) {
271 // Get coefficient associated to i-th shape function
272 auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
273
274 // Get value of i-th transformed reference gradient
275 auto grad = Dune::Functions::flatVectorView(gradients[i]);
276
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) {
281 auto&& c_j = c[j];
282 for(std::size_t k = 0; k < dimV; ++k)
283 re[j*dimV + k] += c_j*grad[k];
284 }
285 }
286 });
287
288 return dy;
289}
290
291
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
297{
298 static_assert(static_size_v<Range> == 1, "Divergence implemented for scalar coefficients only.");
299
300 Range dy;
301 Recursive::forEach(dy, [](auto& v) { v = 0; });
302
303 auto&& node = this->subTreeCache();
304
305 LocalToGlobalBasisAdapter localBasis(node.child(0), this->geometry());
306 auto const& gradients = localBasis.gradientsAt(local);
307
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]);
312
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];
316 }
317
318 return dy;
319}
320
321
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
327{
328 std::size_t comp = this->type_.comp;
329
330 Range dy;
331 Recursive::forEach(dy, [](auto& v) { v = 0; });
332
333 Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
334 {
335 LocalToGlobalBasisAdapter localBasis(node, this->geometry());
336 auto const& partial = localBasis.partialsAt(local, comp);
337
338 // Get range entry associated to this node
339 auto re = Dune::Functions::flatVectorView(hierarchicNodeToRangeMap(tp, dy));
340
341 for (std::size_t i = 0; i < localBasis.size(); ++i) {
342 // Get coefficient associated to i-th shape function
343 auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
344
345 // Get value of i-th transformed reference partial_derivative
346 auto d_comp = Dune::Functions::flatVectorView(partial[i]);
347
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) {
352 auto&& c_j = c[j];
353 for(std::size_t k = 0; k < dimV; ++k)
354 re[j*dimV + k] += c_j*d_comp[k];
355 }
356 }
357 });
358
359 return dy;
360}
361
362
363} // end namespace AMDiS
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