1#ifndef __DUNE_ACFEM_FUNCTIONS_BINDABLETENSORFUNCTION_HH__
2#define __DUNE_ACFEM_FUNCTIONS_BINDABLETENSORFUNCTION_HH__
4#include <dune/fem/function/localfunction/bindable.hh>
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"
15#include "placeholders/placeholdertraits.hh"
16#include "placeholders/localfunctionderivatives.hh"
22 namespace GridFunction {
28 using namespace Literals;
30 using Tensor::indeterminate;
31 using Tensor::TensorTraits;
32 using Tensor::derivative;
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;
47 template<
class Expr,
class Gr
idPart, std::
size_t AutoDiffOrder, std::
size_t IndeterminateId>
48 class BindableTensorFunction;
51 template<
class Expr,
class Gr
idPart, 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> >
58 using TensorSignature =
typename TensorTraits<Expr>::Signature;
59 using RangeType =
typename Tensor::FieldVectorVectorType<FieldType, TensorSignature>;
61 using BaseType = Fem::BindableGridFunction<GridPart, RangeType>;
63 using DomainType =
typename BaseType::DomainType;
64 using JacobianRangeType =
typename BaseType::JacobianRangeType;
65 using HessianRangeType =
typename BaseType::HessianRangeType;
66 using EntityType =
typename BaseType::EntityType;
68 using BaseType::global;
69 using BaseType::entity;
70 using BaseType::gridPart;
74 static_assert(!ExamineOr<Expr, Tensor::IsAutoDiffPlaceholder>::value,
75 "Expr should not contain AD placeholders.");
79 static constexpr std::size_t autoDiffOrder_ =
80 (hasConstantPolynomialDegreeAtMost<0, Expr, IndeterminateId>()
82 : (hasConstantPolynomialDegreeAtMost<1, Expr, IndeterminateId>()
83 ? std::min(AutoDiffOrder, 1UL)
86 using IndeterminateType =
87 std::decay_t<decltype(asExpression(Tensor::indeterminate<IndeterminateId>(std::declval<DomainType>())))>;
89 using ValuesType =
decltype(
90 replaceIndeterminate<IndeterminateId>(std::declval<Expr>(), std::declval<IndeterminateType&>())
93 static_assert(!std::is_rvalue_reference<ValuesType>::value,
94 "Expr must not be an rvalue reference.");
96 using ValuesCSECascade = std::decay_t<
decltype(
97 Tensor::autoDiffCSECascade<autoDiffOrder_>(std::declval<ValuesType&&>(), std::declval<IndeterminateType>())
100 using ValuesCSEType = std::decay_t<
decltype(
101 injectSubExpressions(std::declval<ValuesType&&>(), std::declval<ValuesCSECascade>())
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())
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_))
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_))
158 void bind(
const EntityType& entity)
160 BaseType::bind(entity);
161 forEach(bindablePlaceholders_, [&](
auto&& p) {
168 forEach(unbindablePlaceholders_, [&](
auto&& p) {
174 BaseType::bind(EntityType());
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
183 initIndeterminate(point, ExplicitlyDependsOn<indeterminateId_, ValuesType>{});
184 initPlaceholders(point, valueSetValuePlaceholders_);
185 initPlaceholders(valueIndeterminateSetValuePlaceholders_);
187 valuesCSECascade_.evaluate(valuesCSECascade_.sExpTraits().upTo(0_c));
193 template<
class Point,
195 void evaluate(
const Point& point, RangeType& result)
const
197 result = (*this)(point);
201 template<
class Point,
203 void jacobian(
const Point& point, JacobianRangeType& result)
const
205 if constexpr (autoDiffOrder_ < 1) {
208 initIndeterminate(point, ExplicitlyDependsOn<indeterminateId_, ValuesType>{});
209 initPlaceholders(point, valueSetValuePlaceholders_);
210 initPlaceholders(valueIndeterminateSetValuePlaceholders_);
212 const auto& evalTraits = valuesCSECascade_.sExpTraits().upTo(1_c);
213 valuesCSECascade_.evaluate(evalTraits);
216 assign(result, derivative(evalTraits(values_), indeterminate_));
221 template<
class Point,
223 void hessian(
const Point& point, HessianRangeType& result)
const
225 if constexpr (autoDiffOrder_ < 2) {
228 static_assert(autoDiffOrder_ >= 2,
"We were constructed without 2nd derivatives ...");
230 initIndeterminate(point, ExplicitlyDependsOn<indeterminateId_, ValuesType>{});
231 initPlaceholders(point, valueSetValuePlaceholders_);
232 initPlaceholders(valueIndeterminateSetValuePlaceholders_);
234 const auto& evalTraits = valuesCSECascade_.sExpTraits().upTo(2_c);
235 valuesCSECascade_.evaluate(evalTraits);
238 assign(result, derivative<2>(evalTraits(values_), indeterminate_));
242 ValuesType expression()&&
247 decltype(
auto) expression()&
252 decltype(
auto) expression()
const&
262 template<std::
size_t I>
265 return gridFunction(gridPart(), restriction<0>(values_, Seq<I>{}), requestedADOrder(), indeterminateId());
272 template<std::
size_t I>
275 return gridFunction(gridPart(), restriction<0>(values_, Seq<I>{}), requestedADOrder(), indeterminateId());
282 template<std::
size_t I>
287 restriction<0>(std::move(values_), Seq<I>{}),
293 constexpr unsigned int order()
const
295 return Expressions::polynomialDegree<IndeterminateId>(values_);
298 std::string expressionName()
const
300 return expression_.name();
303 std::string name()
const
308 void setName(
const std::string & name)
313 static constexpr auto indeterminateId()
315 return IndexConstant<indeterminateId_>{};
318 static constexpr auto autoDiffOrder()
320 return IndexConstant<(hasConstantPolynomialDegreeOf<0, Expr, IndeterminateId>()
321 ? Policy::autoDiffLimit()
322 : (hasConstantPolynomialDegreeOf<1, Expr, IndeterminateId>() && autoDiffOrder_ >= 1
323 ? Policy::autoDiffLimit()
324 : autoDiffOrder_))>{};
327 static constexpr auto requestedADOrder()
329 return IndexConstant<requestedADOrder_>{};
332 decltype(
auto) indeterminate()
const& {
333 return indeterminate_;
336 decltype(
auto) indeterminate()& {
337 return indeterminate_;
340 IndeterminateType indeterminate()&& {
341 return indeterminate_;
344 static constexpr auto signature()
346 return typename TensorTraits<ValuesType>::Signature{};
349 void printInfo(std::ostream& out = std::clog)
const
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;
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) {
385 std::clog << t.name() << std::endl;
388 std::clog <<
" #CSE expressions: " << count << std::endl;
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)>
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)
407 tensor(to[0]) = from;
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)
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
425 indeterminate_.operand(0_c) = global(point);
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
437 struct IsRelevantPlaceholder
439 && !Tensor::IsAutoDiffPlaceholder<T>::value
440 && !RefersConst<T>::value
444 template<
class T,
class Cascade>
445 static auto extractPlaceholders(T&& expr, Cascade&& cascade)
448 extract<IsRelevantPlaceholder>(
449 std::forward<T>(expr),
450 extract<IsRelevantPlaceholder>(
451 std::forward<Cascade>(cascade),
454 PredicateWrapper<AreRuntimeEqual>{}
457 PredicateWrapper<AreRuntimeEqual>{}
461 template<
class T,
class Cascade>
462 using PlaceholderTuple = std::decay_t<decltype(extractPlaceholders(std::declval<T>(), std::declval<Cascade>()))>;
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)
473 forLoop<size<Placeholders>()>([&](
auto i) {
474 get<i.value>(placeholders).setValue(point);
481 template<
class Placeholders>
482 void initPlaceholders(
const Placeholders& placeholders)
const
484 forLoop<size<Placeholders>()>([&](
auto i) {
485 get<i.value>(placeholders).setValue(indeterminate_);
494 mutable IndeterminateType indeterminate_;
501 ValuesType expression_;
511 mutable ValuesCSECascade valuesCSECascade_;
512 ValuesCSEType values_;
531 PlaceholderTuple<ValuesCSEType&, ValuesCSECascade&> valuePlaceholders_;
532 FilteredTuple<HasBind,
decltype(valuePlaceholders_)> bindablePlaceholders_;
533 FilteredTuple<HasUnbind,
decltype(valuePlaceholders_)> unbindablePlaceholders_;
535 FilteredTuple<HasSetValue,
decltype(valuePlaceholders_)> valueSetValuePlaceholders_;
537 FilteredTuple<HasIndeterminateSetValue,
decltype(valuePlaceholders_)> valueIndeterminateSetValuePlaceholders_;
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
Terminals may derive from this class to express that they are expressions.
Definition: terminal.hh:25
Definition: quadraturepoint.hh:29