DUNE-ACFEM (unstable)

bindabletensorfunction.hh
1#ifndef __DUNE_ACFEM_FUNCTIONS_BINDABLETENSORFUNCTION_HH__
2#define __DUNE_ACFEM_FUNCTIONS_BINDABLETENSORFUNCTION_HH__
3
4#include <dune/fem/function/localfunction/bindable.hh>
5
6#include "../common/ostream.hh"
7#include "../common/quadraturepoint.hh"
8#include "../expressions/examine.hh"
9#include "../expressions/polynomialtraits.hh"
10#include "../tensors/optimization/subexpression.hh"
11#include "../tensors/bindings/dune/fieldvectortraits.hh"
12#include "../tensors/ostream.hh"
13
14#include "policy.hh"
15#include "placeholders/placeholdertraits.hh"
16#include "placeholders/localfunctionderivatives.hh"
17
18namespace Dune {
19
20 namespace ACFem {
21
22 namespace GridFunction {
23
28 using namespace Literals;
29
30 using Tensor::indeterminate;
31 using Tensor::TensorTraits;
32 using Tensor::derivative;
33 using Tensor::Seq;
34 using Expressions::replaceIndeterminate;
35 using Expressions::IsClosure;
36 using Expressions::extract;
37 using Expressions::extractPlaceholders;
38 using Expressions::extractIndeterminates;
39 using Expressions::PlaceholderTuple;
40 using Expressions::ExplicitlyDependsOn;
42 using Expressions::finalize;
44 using Expressions::AreRuntimeEqual;
45 using Expressions::ExamineOr;
46
47 template<class Expr, class GridPart, std::size_t AutoDiffOrder, std::size_t IndeterminateId>
48 class BindableTensorFunction;
49
51 template<class Expr, class GridPart, std::size_t AutoDiffOrder = Policy::autoDiffLimit(), std::size_t IndeterminateId = Policy::indeterminateId()>
53 : public Fem::BindableGridFunction<GridPart, Tensor::FieldVectorVectorType<FloatingPointClosure<typename TensorTraits<Expr>::FieldType>, typename TensorTraits<Expr>::Signature> >
54 , public Expressions::SelfExpression<BindableTensorFunction<Expr, GridPart, AutoDiffOrder, IndeterminateId> >
55 {
56 public:
58 using TensorSignature = typename TensorTraits<Expr>::Signature;
59 using RangeType = typename Tensor::FieldVectorVectorType<FieldType, TensorSignature>;
60 private:
61 using BaseType = Fem::BindableGridFunction<GridPart, RangeType>;
62 public:
63 using DomainType = typename BaseType::DomainType;
64 using JacobianRangeType = typename BaseType::JacobianRangeType;
65 using HessianRangeType = typename BaseType::HessianRangeType;
66 using EntityType = typename BaseType::EntityType;
67
68 using BaseType::global;
69 using BaseType::entity;
70 using BaseType::gridPart;
71
72 static constexpr std::size_t indeterminateId_ = IndeterminateId;
73
74 static_assert(!ExamineOr<Expr, Tensor::IsAutoDiffPlaceholder>::value,
75 "Expr should not contain AD placeholders.");
76
78 static constexpr std::size_t requestedADOrder_ = AutoDiffOrder;
79 static constexpr std::size_t autoDiffOrder_ =
80 (hasConstantPolynomialDegreeAtMost<0, Expr, IndeterminateId>()
81 ? 0
82 : (hasConstantPolynomialDegreeAtMost<1, Expr, IndeterminateId>()
83 ? std::min(AutoDiffOrder, 1UL)
84 : AutoDiffOrder));
85
86 using IndeterminateType =
87 std::decay_t<decltype(asExpression(Tensor::indeterminate<IndeterminateId>(std::declval<DomainType>())))>;
88
89 using ValuesType = decltype(
90 replaceIndeterminate<IndeterminateId>(std::declval<Expr>(), std::declval<IndeterminateType&>())
91 );
92
93 static_assert(!std::is_rvalue_reference<ValuesType>::value,
94 "Expr must not be an rvalue reference.");
95
96 using ValuesCSECascade = std::decay_t<decltype(
97 Tensor::autoDiffCSECascade<autoDiffOrder_>(std::declval<ValuesType&&>(), std::declval<IndeterminateType>())
98 )>;
99
100 using ValuesCSEType = std::decay_t<decltype(
101 injectSubExpressions(std::declval<ValuesType&&>(), std::declval<ValuesCSECascade>())
102 )>;
103
105 BindableTensorFunction(Expr&& t, const GridPart& gridPart)
106 : BaseType(gridPart)
107 , indeterminate_(asExpression(Tensor::indeterminate<IndeterminateId>(DomainType(0))))
108 , expression_(replaceIndeterminate<indeterminateId_>(std::forward<Expr>(t), indeterminate_))
109 , valuesCSECascade_(Tensor::autoDiffCSECascade<autoDiffOrder_>(static_cast<ValuesType&&>(expression_), indeterminate_))
110 , values_(injectSubExpressions(static_cast<ValuesType>(expression_), valuesCSECascade_))
111 , valuePlaceholders_(extractPlaceholders(values_, valuesCSECascade_))
112 , bindablePlaceholders_(filteredTuple<HasBind>(valuePlaceholders_))
113 , unbindablePlaceholders_(filteredTuple<HasUnbind>(valuePlaceholders_))
114 , valueSetValuePlaceholders_(filteredTuple<HasSetValue>(valuePlaceholders_))
115 , valueIndeterminateSetValuePlaceholders_(filteredTuple<HasIndeterminateSetValue>(valuePlaceholders_))
116 , name_(expression_.name())
117 {}
118
125 : BaseType(other)
126 , indeterminate_(asExpression(Tensor::indeterminate<IndeterminateId>(DomainType())))
127 , expression_(replaceIndeterminate<indeterminateId_>(const_cast<ValuesType&&>(other.expression_), indeterminate_))
128 , valuesCSECascade_(Tensor::autoDiffCSECascade<autoDiffOrder_>(static_cast<ValuesType&&>(expression_), indeterminate_))
129 , values_(injectSubExpressions(static_cast<ValuesType>(expression_), valuesCSECascade_))
130 , valuePlaceholders_(extractPlaceholders(values_, valuesCSECascade_))
131 , bindablePlaceholders_(filteredTuple<HasBind>(valuePlaceholders_))
132 , unbindablePlaceholders_(filteredTuple<HasUnbind>(valuePlaceholders_))
133 , valueSetValuePlaceholders_(filteredTuple<HasSetValue>(valuePlaceholders_))
134 , valueIndeterminateSetValuePlaceholders_(filteredTuple<HasIndeterminateSetValue>(valuePlaceholders_))
135 , name_(other.name_)
136 {}
137
140 : BaseType(other)
141 , indeterminate_(asExpression(Tensor::indeterminate<IndeterminateId>(DomainType())))
142 , expression_(replaceIndeterminate<indeterminateId_>(const_cast<ValuesType&&>(other.expression_), indeterminate_))
143 , valuesCSECascade_(Tensor::autoDiffCSECascade<autoDiffOrder_>(static_cast<ValuesType&&>(expression_), indeterminate_))
144 , values_(injectSubExpressions(static_cast<ValuesType>(expression_), valuesCSECascade_))
145 , valuePlaceholders_(extractPlaceholders(values_, valuesCSECascade_))
146 , bindablePlaceholders_(filteredTuple<HasBind>(valuePlaceholders_))
147 , unbindablePlaceholders_(filteredTuple<HasUnbind>(valuePlaceholders_))
148 , valueSetValuePlaceholders_(filteredTuple<HasSetValue>(valuePlaceholders_))
149 , valueIndeterminateSetValuePlaceholders_(filteredTuple<HasIndeterminateSetValue>(valuePlaceholders_))
150 , name_(std::move(other.name_))
151 {}
152
158 void bind(const EntityType& entity)
159 {
160 BaseType::bind(entity);
161 forEach(bindablePlaceholders_, [&](auto&& p) {
162 p.bind(entity);
163 });
164 }
165
166 void unbind()
167 {
168 forEach(unbindablePlaceholders_, [&](auto&& p) {
169 p.unbind();
170 });
171 BaseType::unbind();
172#ifndef NDEBUG
173 // in order to catch bogus bind cascades
174 BaseType::bind(EntityType());
175#endif
176 }
177
179 template<class Point,
180 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
181 const auto& operator()(const Point& point) const
182 {
183 initIndeterminate(point, ExplicitlyDependsOn<indeterminateId_, ValuesType>{});
184 initPlaceholders(point, valueSetValuePlaceholders_);
185 initPlaceholders(valueIndeterminateSetValuePlaceholders_);
186
187 valuesCSECascade_.evaluate(valuesCSECascade_.sExpTraits().upTo(0_c));
188
189 return values_;
190 }
191
193 template<class Point,
194 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
195 void evaluate(const Point& point, RangeType& result) const
196 {
197 result = (*this)(point);
198 }
199
201 template<class Point,
202 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
203 void jacobian(const Point& point, JacobianRangeType& result) const
204 {
205 if constexpr (autoDiffOrder_ < 1) {
206 result = 0;
207 } else {
208 initIndeterminate(point, ExplicitlyDependsOn<indeterminateId_, ValuesType>{});
209 initPlaceholders(point, valueSetValuePlaceholders_);
210 initPlaceholders(valueIndeterminateSetValuePlaceholders_);
211
212 const auto& evalTraits = valuesCSECascade_.sExpTraits().upTo(1_c);
213 valuesCSECascade_.evaluate(evalTraits);
214
215 // The call to derivative() will just peek out the pre-evaluated derivative
216 assign(result, derivative(evalTraits(values_), indeterminate_));
217 }
218 }
219
221 template<class Point,
222 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
223 void hessian(const Point& point, HessianRangeType& result) const
224 {
225 if constexpr (autoDiffOrder_ < 2) {
226 result = 0;
227 } else {
228 static_assert(autoDiffOrder_ >= 2, "We were constructed without 2nd derivatives ...");
229
230 initIndeterminate(point, ExplicitlyDependsOn<indeterminateId_, ValuesType>{});
231 initPlaceholders(point, valueSetValuePlaceholders_);
232 initPlaceholders(valueIndeterminateSetValuePlaceholders_);
233
234 const auto& evalTraits = valuesCSECascade_.sExpTraits().upTo(2_c);
235 valuesCSECascade_.evaluate(evalTraits);
236
237 // The call to derivative() will just peek out the pre-evaluated derivative
238 assign(result, derivative<2>(evalTraits(values_), indeterminate_));
239 }
240 }
241
242 ValuesType expression()&&
243 {
244 return expression_;
245 }
246
247 decltype(auto) expression()&
248 {
249 return expression_;
250 }
251
252 decltype(auto) expression() const&
253 {
254 return expression_;
255 }
256
262 template<std::size_t I>
264 {
265 return gridFunction(gridPart(), restriction<0>(values_, Seq<I>{}), requestedADOrder(), indeterminateId());
266 }
267
272 template<std::size_t I>
273 auto operator[](const IndexConstant<I>)const&
274 {
275 return gridFunction(gridPart(), restriction<0>(values_, Seq<I>{}), requestedADOrder(), indeterminateId());
276 }
277
282 template<std::size_t I>
283 auto operator[](const IndexConstant<I>)const&&
284 {
285 return gridFunction(
286 gridPart(),
287 restriction<0>(std::move(values_), Seq<I>{}),
288 requestedADOrder(),
289 indeterminateId()
290 );
291 }
292
293 constexpr unsigned int order() const
294 {
295 return Expressions::polynomialDegree<IndeterminateId>(values_);
296 }
297
298 std::string expressionName() const
299 {
300 return expression_.name();
301 }
302
303 std::string name() const
304 {
305 return name_;
306 }
307
308 void setName(const std::string & name)
309 {
310 name_ = name;
311 }
312
313 static constexpr auto indeterminateId()
314 {
315 return IndexConstant<indeterminateId_>{};
316 }
317
318 static constexpr auto autoDiffOrder()
319 {
320 return IndexConstant<(hasConstantPolynomialDegreeOf<0, Expr, IndeterminateId>()
321 ? Policy::autoDiffLimit()
322 : (hasConstantPolynomialDegreeOf<1, Expr, IndeterminateId>() && autoDiffOrder_ >= 1
323 ? Policy::autoDiffLimit()
324 : autoDiffOrder_))>{};
325 }
326
327 static constexpr auto requestedADOrder()
328 {
329 return IndexConstant<requestedADOrder_>{};
330 }
331
332 decltype(auto) indeterminate() const& {
333 return indeterminate_;
334 }
335
336 decltype(auto) indeterminate()& {
337 return indeterminate_;
338 }
339
340 IndeterminateType indeterminate()&& {
341 return indeterminate_;
342 }
343
344 static constexpr auto signature()
345 {
346 return typename TensorTraits<ValuesType>::Signature{};
347 }
348
349 void printInfo(std::ostream& out = std::clog) const
350 {
351 //std::clog << "Memory Footprint: " << sizeof(*this) << std::endl;
352 std::clog << "IndeterminateType: " << typeString<IndeterminateType>() << std::endl;
353 forLoop<3>([](auto i) {
354 using I = decltype(i);
355 std::clog << "Constant Degree of " << I::value << ": " << hasConstantPolynomialDegreeOf<I::value, ValuesType, IndeterminateId>() << std::endl;
356 });
357 std::clog << "Estimated order: " << order() << std::endl;
358 std::clog << "Value Placeholders: "
359 << " #" << size(valuePlaceholders_) << ": "
360 << typeString<decltype(valuePlaceholders_)>() << std::endl;
361 std::clog << " bindable: "
362 << " #" << size(bindablePlaceholders_) << ": "
363 << typeString<decltype(bindablePlaceholders_)>() << std::endl;
364 std::clog << " unbindable: "
365 << " #" << size(unbindablePlaceholders_) << ": "
366 << typeString<decltype(unbindablePlaceholders_)>() << std::endl;
367 std::clog << " set value: "
368 << " #" << size(valueSetValuePlaceholders_) << ": "
369 << typeString<decltype(valueSetValuePlaceholders_)>() << std::endl;
370 std::clog << " set indet: "
371 << " #" << size(valueIndeterminateSetValuePlaceholders_) << ": "
372 << typeString<decltype(valueIndeterminateSetValuePlaceholders_)>() << std::endl;
373 std::clog << "ValuesType:" << std::endl
374 << " " << values_.name() << std::endl
375 << " " << typeString<decltype(values_)>() << std::endl;
376 std::clog << " depends explicitly on " << indeterminate_.name() << ": " << ExplicitlyDependsOn<indeterminateId_, ValuesType>::value << std::endl;
377 std::clog << " is zero: " << ExpressionTraits<ValuesType>::isZero << std::endl;
378 std::clog << " tensor signature: " << TensorTraits<ValuesType>::signature() << std::endl;
379 std::clog << " fem signature: " << TensorTraits<RangeType>::signature() << std::endl;
380 std::clog << " CSE expressions: " << std::endl;
381 std::size_t count = 0;
382 forLoop<ValuesCSECascade::depth()>([this,&count](auto i) {
383 forEach(valuesCSECascade_.template get<i.value>(),[&count](const auto& t) {
384 ++count;
385 std::clog << t.name() << std::endl;
386 });
387 });
388 std::clog << " #CSE expressions: " << count << std::endl;
389 }
390
391 private:
392 template<class FemResult, class TensorResult>
393 struct NeedsScalarDisambiguation
394 : BoolConstant<((TensorTraits<FemResult>::rank - TensorTraits<TensorResult>::rank) == 1
395 && TensorTraits<FemResult>::template dim<0>() == 1)>
396 {};
397
403 template<class FemResult, class TensorResult,
404 std::enable_if_t<(NeedsScalarDisambiguation<FemResult, TensorResult>::value), int> = 0>
405 static void assign(FemResult& to, TensorResult&& from)
406 {
407 tensor(to[0]) = from;
408 }
409
410 template<class FemResult, class TensorResult,
411 std::enable_if_t<(!NeedsScalarDisambiguation<FemResult, TensorResult>::value), int> = 0>
412 static void assign(FemResult& to, TensorResult&& from)
413 {
414 tensor(to) = from;
415 }
416
418 template<
419 class Point, class IsDependent,
420 std::enable_if_t<((IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value)
421 && IsDependent::value), int> = 0>
422 void initIndeterminate(const Point& point, IsDependent) const
423 {
424 //assert(entity() != EntityType());
425 indeterminate_.operand(0_c) = global(point);
426 }
427
429 template<
430 class Point, class IsDependent,
431 std::enable_if_t<((IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value)
432 && !IsDependent::value), int> = 0>
433 void initIndeterminate(const Point& point, IsDependent) const
434 {}
435
436 template<class T>
437 struct IsRelevantPlaceholder
438 : BoolConstant<(IsPlaceholderExpression<T>::value
439 && !Tensor::IsAutoDiffPlaceholder<T>::value
440 && !RefersConst<T>::value
441 )>
442 {};
443
444 template<class T, class Cascade>
445 static auto extractPlaceholders(T&& expr, Cascade&& cascade)
446 {
447 return
448 extract<IsRelevantPlaceholder>(
449 std::forward<T>(expr),
450 extract<IsRelevantPlaceholder>(
451 std::forward<Cascade>(cascade),
452 std::tuple<>{},
453 FalseType{}, // recursive
454 PredicateWrapper<AreRuntimeEqual>{}
455 ),
456 FalseType{}, // recursive
457 PredicateWrapper<AreRuntimeEqual>{}
458 );
459 }
460
461 template<class T, class Cascade>
462 using PlaceholderTuple = std::decay_t<decltype(extractPlaceholders(std::declval<T>(), std::declval<Cascade>()))>;
463
469 template<class Point, class Placeholders,
470 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
471 static void initPlaceholders(const Point& point, const Placeholders& placeholders)
472 {
473 forLoop<size<Placeholders>()>([&](auto i) {
474 get<i.value>(placeholders).setValue(point);
475 });
476 }
477
481 template<class Placeholders>
482 void initPlaceholders(const Placeholders& placeholders) const
483 {
484 forLoop<size<Placeholders>()>([&](auto i) {
485 get<i.value>(placeholders).setValue(indeterminate_);
486 });
487 }
488
494 mutable IndeterminateType indeterminate_;
495
501 ValuesType expression_;
502
511 mutable ValuesCSECascade valuesCSECascade_;
512 ValuesCSEType values_;
514
531 PlaceholderTuple<ValuesCSEType&, ValuesCSECascade&> valuePlaceholders_;
532 FilteredTuple<HasBind, decltype(valuePlaceholders_)> bindablePlaceholders_;
533 FilteredTuple<HasUnbind, decltype(valuePlaceholders_)> unbindablePlaceholders_;
534
535 FilteredTuple<HasSetValue, decltype(valuePlaceholders_)> valueSetValuePlaceholders_;
536
537 FilteredTuple<HasIndeterminateSetValue, decltype(valuePlaceholders_)> valueIndeterminateSetValuePlaceholders_;
538 std::string name_;
540 };
541
543
544 } // NS GridFunction
545
546 } // NS ACFem
547
548} // NS Dune
549
550#endif // __DUNE_ACFEM_FUNCTIONS_BINDABLETENSORFUNCTION_HH__
Wrap a tensor into a Fem::BindableGridFunction.
Definition: bindabletensorfunction.hh:55
void bind(const EntityType &entity)
Bind ourselves and all placeholders contained in the value expression.
Definition: bindabletensorfunction.hh:158
static constexpr std::size_t indeterminateId_
Id of the indeterminate.
Definition: bindabletensorfunction.hh:72
BindableTensorFunction(BindableTensorFunction &&other)
Move CTOR.
Definition: bindabletensorfunction.hh:139
auto operator[](const IndexConstant< I >) &
Compile-time constant access to the given component.
Definition: bindabletensorfunction.hh:263
auto operator[](const IndexConstant< I >) const &&
Compile-time constant access to the given access.
Definition: bindabletensorfunction.hh:283
BindableTensorFunction(const BindableTensorFunction &other)
Copy CTOR.
Definition: bindabletensorfunction.hh:124
BindableTensorFunction(Expr &&t, const GridPart &gridPart)
Standard constructor from tensor expression and GridPart.
Definition: bindabletensorfunction.hh:105
auto operator[](const IndexConstant< I >) const &
Compile-time constant access to the given access.
Definition: bindabletensorfunction.hh:273
void hessian(const Point &point, HessianRangeType &result) const
Dune::Fem::BindableGridFunction::hessian()
Definition: bindabletensorfunction.hh:223
static constexpr std::size_t requestedADOrder_
Max. order of AD differentiation.
Definition: bindabletensorfunction.hh:78
void jacobian(const Point &point, JacobianRangeType &result) const
Dune::Fem::BindableGridFunction::jacobian()
Definition: bindabletensorfunction.hh:203
void evaluate(const Point &point, RangeType &result) const
Dune::Fem::BindableGridFunction::evaluate()
Definition: bindabletensorfunction.hh:195
constexpr decltype(auto) asExpression(T &&t)
Return a non-closure expression as is.
Definition: interface.hh:122
typename ClosureTraits< T >::ExpressionType EnclosedType
Type alias shortcut to get hold of the type contained in an "expression closure".
Definition: interface.hh:77
typename FloatingPointClosureHelper< T >::Type FloatingPointClosure
Template alias.
Definition: fieldpromotion.hh:74
auto gridFunction(const GridPart &gridPart, T &&t, IndexConstant< MaxOrder > maxOrder=IndexConstant< MaxOrder >{}, IndexConstant< IndeterminateId > id=IndexConstant< IndeterminateId >{})
Generate a BindableGridFunction which wraps a copy of T.
Definition: gridfunction.hh:23
constexpr decltype(auto) get(T &&t, IndexConstant< I >=IndexConstant< I >{})
Access to the i-the element.
Definition: access.hh:129
auto & assign(T1 &t1, T2 &&t2)
Assign one tuple-alike to another by looping over the elements.
Definition: assign.hh:40
constexpr std::size_t size()
Gives the number of elements in tuple-likes and std::integer_sequence.
Definition: size.hh:73
auto filteredTuple(Tuple &&t)
Generate a sub-tuple of the elements where F<TupleElement<N, Tuple> > is derived from std::true_type,...
Definition: subtuple.hh:224
std::decay_t< decltype(filteredTuple< F >(std::declval< Tuple >()))> FilteredTuple
Generate the type of the return value of filteredTuple().
Definition: subtuple.hh:233
Constant< bool, V > BoolConstant
Short-cut for integral constant of type bool.
Definition: types.hh:48
Constant< std::size_t, V > IndexConstant
Short-cut for integral constant of type std::size_t.
Definition: types.hh:44
BoolConstant< false > FalseType
Alias for std::false_type.
Definition: types.hh:110
STL namespace.
Terminals may derive from this class to express that they are expressions.
Definition: terminal.hh:25
Definition: quadraturepoint.hh:29
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jan 3, 23:40, 2025)