diff --git a/src/amdis/FileWriter.hpp b/src/amdis/FileWriter.hpp index a8d95f2552064d8c01eba82c56863b2dfda3796c..6f6488af74291444aa7693d4a0935db3961f05ea 100644 --- a/src/amdis/FileWriter.hpp +++ b/src/amdis/FileWriter.hpp @@ -43,14 +43,17 @@ namespace AMDiS constexpr std::size_t VTKFieldSize = Size<Range>; - template <class GlobalBasis, class RangeType, class TreePath> + template <class GB, class VT, class TP> class FileWriter : public FileWriterInterface { - private: // typedefs and static constants + using GlobalBasis = GB; + using ValueType = VT; + using TreePath = TP; + + private: using GridView = typename GlobalBasis::GridView; - using DiscreteFunction = AMDiS::DiscreteFunction<GlobalBasis,RangeType,TreePath>; - using Range = typename DiscreteFunction::Range; + using Range = typename DiscreteFunction<GB,VT,TP>::Range; /// Dimension of the mesh static constexpr int dim = GridView::dimension; @@ -60,26 +63,25 @@ namespace AMDiS public: /// Constructor. - FileWriter(std::string const& baseName, - DiscreteFunction const& discreteFct) - : FileWriterInterface(baseName) + FileWriter(std::string const& name, DiscreteFunction<GB,VT,TP> const& discreteFct) + : FileWriterInterface(name) , discreteFct_(discreteFct) , animation_(false) { - Parameters::get(baseName + "->ParaView animation", animation_); + Parameters::get(name + "->ParaView animation", animation_); - int m = Parameters::get<int>(baseName + "->ParaView mode").value_or(0); + int m = Parameters::get<int>(name + "->ParaView mode").value_or(0); mode_ = (m == 0 ? Dune::VTK::ascii : Dune::VTK::appendedraw); - init(baseName, ValueCategory_t<Range>{}); + init(name, ValueCategory_t<Range>{}); } template <class ValueCat> void init(std::string const& baseName, ValueCat) { - int subSampling = Parameters::get<int>(baseName + "->subsampling").value_or(0); - if (subSampling > 0) - vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(gridView(), subSampling); + auto subSampling = Parameters::get<int>(baseName + "->subsampling"); + if (subSampling) + vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(gridView(), subSampling.value()); else vtkWriter_ = std::make_shared<Dune::VTKWriter<GridView>>(gridView()); @@ -109,7 +111,7 @@ namespace AMDiS } private: - DiscreteFunction discreteFct_; + DiscreteFunction<GB,VT,TP> discreteFct_; std::shared_ptr<Dune::VTKWriter<GridView>> vtkWriter_; std::shared_ptr<Dune::VTKSequenceWriter<GridView>> vtkSeqWriter_; @@ -122,12 +124,12 @@ namespace AMDiS }; - template <class GlobalBasis, class Range, class TreePath> - std::shared_ptr<FileWriter<GlobalBasis,Range,TreePath>> - makeFileWriterPtr(std::string baseName, - DiscreteFunction<GlobalBasis,Range,TreePath> const& discreteFct) + template <class GlobalBasis, class ValueType, class TreePath> + std::shared_ptr<FileWriter<GlobalBasis,ValueType,TreePath>> + makeFileWriterPtr(std::string const& name, + DiscreteFunction<GlobalBasis,ValueType,TreePath> const& discreteFct) { - return std::make_shared<FileWriter<GlobalBasis,Range,TreePath>>(baseName, discreteFct); + return std::make_shared<FileWriter<GlobalBasis,ValueType,TreePath>>(name, discreteFct); } } // end namespace AMDiS diff --git a/src/amdis/gridfunctions/DOFVectorView.hpp b/src/amdis/gridfunctions/DOFVectorView.hpp index c39cf8df1d37a3735215caf8367caaa0f083cfb1..734c5e1ebc35037d2a13f252b35f9dae9b972e04 100644 --- a/src/amdis/gridfunctions/DOFVectorView.hpp +++ b/src/amdis/gridfunctions/DOFVectorView.hpp @@ -6,19 +6,24 @@ namespace AMDiS { /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction - template <class GlobalBasisType, class RangeType, class TreePathType> + template <class GB, class VT, class TP> class DOFVectorView - : public DiscreteFunction<GlobalBasisType, RangeType, TreePathType> + : public DiscreteFunction<GB, VT, TP> { using Self = DOFVectorView; - using Super = DiscreteFunction<GlobalBasisType, RangeType, TreePathType>; + using Super = DiscreteFunction<GB, VT, TP>; - using GlobalBasis = GlobalBasisType; - using TreePath = TreePathType; + using GlobalBasis = GB; + using TreePath = TP; public: + /// Constructor forwards to the treePath constructor, with empty TreePath + DOFVectorView(DOFVector<GB,VT>& dofVector) + : DOFVectorView{dofVector, Dune::TypeTree::hybridTreePath()} + {} + /// Constructor. Stores a pointer to the mutable `dofvector`. - DOFVectorView(DOFVector<GlobalBasis,RangeType>& dofVector, TreePath const& treePath) + DOFVectorView(DOFVector<GB,VT>& dofVector, TP const& treePath) : Super(dofVector, treePath) , mutableDofVector_(&dofVector) {} @@ -26,6 +31,13 @@ namespace AMDiS public: /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no /// reference to this DOFVector in the expression. + /** + * **Example:** + * ``` + * auto v = makeDOFVectorView(prob.getSolutionVector(),0); + * v.interpolate([](auto const& x) { return x[0]; }); + * ``` + **/ template <class Expr> void interpolate_noalias(Expr&& expr) { @@ -34,15 +46,25 @@ namespace AMDiS auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), basis.gridView()); - Dune::Functions::interpolate(basis, treePath, coefficients(), std::forward<decltype(gridFct)>(gridFct)); + Dune::Functions::interpolate(basis, treePath, coefficients(), + std::forward<decltype(gridFct)>(gridFct)); } - /// Interpolation of GridFunction to DOFVector + /// \brief Interpolation of GridFunction to DOFVector + /** + * **Example:** + * ``` + * auto v = makeDOFVectorView(prob.getSolutionVector(),0); + * v.interpolate(v + [](auto const& x) { return x[0]; }); + * ``` + * Allows to have a reference to the DOFVector in the expression, e.g. as + * \ref DiscreteFunction or \ref gradientAtQP() of a DiscreteFunction. + **/ template <class Expr> void interpolate(Expr&& expr) { // create temporary copy of data - DOFVector<GlobalBasis,RangeType> tmp(coefficients()); + DOFVector<GB,VT> tmp(coefficients()); Self tmpView{tmp, this->treePath()}; tmpView.interpolate_noalias(std::forward<Expr>(expr)); @@ -58,6 +80,7 @@ namespace AMDiS return *this; } + /// \brief interpolate `(*this) + expr` to DOFVector template <class Expr> DOFVectorView& operator+=(Expr&& expr) { @@ -65,6 +88,7 @@ namespace AMDiS return *this; } + /// \brief interpolate `(*this) - expr` to DOFVector template <class Expr> DOFVectorView& operator-=(Expr&& expr) { @@ -72,31 +96,37 @@ namespace AMDiS return *this; } - /// Return the mutable DOFVector - DOFVector<GlobalBasis,RangeType>& coefficients() { return *mutableDofVector_; } + DOFVector<GB,VT>& coefficients() { return *mutableDofVector_; } /// Return the const DOFVector using Super::coefficients; protected: - DOFVector<GlobalBasis,RangeType>* mutableDofVector_; + DOFVector<GB,VT>* mutableDofVector_; }; +#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION + // Deduction guide for DOFVectorView class + template <class GlobalBasis, class ValueType> + DOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector) + -> DOFVectorView<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>; +#endif + /// A Generator for a mutable \ref DOFVectorView - template <class GlobalBasis, class RangeType, class TreePath> - auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType>& dofVector, TreePath const& treePath) + template <class GlobalBasis, class ValueType, class TreePath> + auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector, TreePath const& treePath) { - return DOFVectorView<GlobalBasis, RangeType, TreePath>{dofVector, treePath}; + return DOFVectorView<GlobalBasis, ValueType, TreePath>{dofVector, treePath}; } /// A Generator for a mutable \ref DOFVectorView - template <class GlobalBasis, class RangeType> - auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType>& dofVector) + template <class GlobalBasis, class ValueType> + auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector) { auto treePath = Dune::TypeTree::hybridTreePath(); - return DOFVectorView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath}; + return DOFVectorView<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath}; } } // end namespace AMDiS diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp index 396fe04c842e7b1d35958247bbe6db63606fb81a..3f8dc51055328c8cd3479a9c0aca71866a62f4b5 100644 --- a/src/amdis/gridfunctions/DiscreteFunction.hpp +++ b/src/amdis/gridfunctions/DiscreteFunction.hpp @@ -17,13 +17,19 @@ namespace AMDiS /// \brief A view on a subspace of a \ref DOFVector /** * \ingroup GridFunctions + * + * \tparam GB GlobalBasis type that models \ref Dune::Functions::Concept::GlobalBasis + * \tparam VT Coefficient type of the DOFVector + * \tparam TP TreePath type a realization of \ref Dune::TypeTree::HybridTreePath **/ - template <class GlobalBasisType, class RangeType, class TreePathType> + template <class GB, class VT, class TP> class DiscreteFunction { - public: - using GlobalBasis = GlobalBasisType; - using TreePath = TreePathType; + static_assert(std::is_arithmetic<VT>::value, ""); + + private: + using GlobalBasis = GB; + using TreePath = TP; using Tree = typename GlobalBasis::LocalView::Tree; using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>; @@ -32,159 +38,39 @@ namespace AMDiS using GridView = typename GlobalBasis::GridView; using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>; + public: + /// Global coordinates of the EntitySet using Domain = typename EntitySet::GlobalCoordinate; - using Range = RangeType_t<SubTree>; - static_assert(std::is_arithmetic<RangeType>::value, ""); - // Don't know how to determine Range with non-trivial RangeType - - using RawSignature = typename Dune::Functions::SignatureTraits<Range(Domain)>::RawSignature; - using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>; - using DerivativeRange = typename DerivativeTraits::Range; - using LocalDomain = typename EntitySet::LocalCoordinate; - using Element = typename EntitySet::Element; - using Geometry = typename Element::Geometry; + /// Range type of this DiscreteFunction + using Range = RangeType_t<SubTree,VT>; + /// \brief This GridFunction has no derivative function, it can be created + /// by \ref DiscreteGridFunction. enum { hasDerivative = false }; - public: // a local view on the gradients - - /// A LocalFunction representing the derivative of the DOFVector - class GradientLocalFunction - { - public: - using Domain = LocalDomain; - using Range = DerivativeRange; - - enum { hasDerivative = false }; - - private: - using LocalView = typename GlobalBasis::LocalView; - - public: - GradientLocalFunction(DiscreteFunction const& globalFunction) - : globalFunction_(&globalFunction) - , localView_(globalFunction_->basis().localView()) - , subTree_(&child(localView_.tree(), globalFunction_->treePath())) - {} - - void bind(Element const& element) - { - localView_.bind(element); - geometry_.emplace(element.geometry()); - bound_ = true; - } - - void unbind() - { - localView_.unbind(); - geometry_.reset(); - bound_ = false; - } - - /// Evaluate Gradient at bound element in local coordinates - Range operator()(Domain const& x) const; - - friend int order(GradientLocalFunction const& self) - { - assert( self.bound_ ); - return std::max(0, polynomialDegree(*self.subTree_)-1); - } - - /// Return the bound element - Element const& localContext() const - { - assert( bound_ ); - return localView_.element(); - } - - private: - DiscreteFunction const* globalFunction_; - LocalView localView_; - - SubTree const* subTree_; - - Dune::Std::optional<Geometry> geometry_; - bool bound_ = false; - }; - - - public: // a local view on the values - - /// A LocalFunction, i.e., an element local view on the DOFVector - class LocalFunction - { - public: - using Domain = typename DiscreteFunction::LocalDomain; - using Range = typename DiscreteFunction::Range; - - enum { hasDerivative = true }; - - private: - using LocalView = typename GlobalBasis::LocalView; - - public: - LocalFunction(DiscreteFunction const& globalFunction) - : globalFunction_(&globalFunction) - , localView_(globalFunction_->basis().localView()) - , subTree_(&child(localView_.tree(), globalFunction_->treePath())) - {} - - void bind(Element const& element) - { - localView_.bind(element); - bound_ = true; - } - - void unbind() - { - localView_.unbind(); - bound_ = false; - } - - /// Evaluate LocalFunction at bound element in local coordinates - Range operator()(Domain const& x) const; - - /// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction - friend GradientLocalFunction derivative(LocalFunction const& localFunction) - { - static_assert(isValidRange<DerivativeTraits>(), - "Derivative of DOFVector not defined."); - return GradientLocalFunction{*localFunction.globalFunction_}; - } - - friend int order(LocalFunction const& self) - { - assert( self.bound_ ); - return polynomialDegree(*self.subTree_); - } - - /// Return the bound element - Element const& localContext() const - { - assert( bound_ ); - return localView_.element(); - } - - private: - DiscreteFunction const* globalFunction_; - LocalView localView_; - SubTree const* subTree_; - - bool bound_ = false; - }; + public: + /// A LocalFunction representing the derivative of the DOFVector on a bound element + class GradientLocalFunction; + /// A LocalFunction representign the value the DOFVector on a bound element + class LocalFunction; public: + /// Constructor forwards to the treePath constructor, with empty TreePath + DiscreteFunction(DOFVector<GB,VT> const& dofVector) + : DiscreteFunction{dofVector, Dune::TypeTree::hybridTreePath()} + {} + /// Constructor. Stores a pointer to the dofVector and a copy of the treePath. - DiscreteFunction(DOFVector<GlobalBasis,RangeType> const& dofVector, TreePath const& treePath) + DiscreteFunction(DOFVector<GB,VT> const& dofVector, TP const& treePath) : dofVector_(&dofVector) , treePath_(treePath) , entitySet_(dofVector.basis().gridView()) , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(dofVector.basis(), treePath)) {} - /// Evaluate DiscreteFunction in global coordinates. NOTE: expensive + /// \brief Evaluate DiscreteFunction in global coordinates. NOTE: expensive Range operator()(Domain const& x) const; /// \brief Create a local function for this view on the DOFVector. \relates LocalFunction @@ -193,51 +79,59 @@ namespace AMDiS return LocalFunction{self}; } + /// \brief Return a \ref Dune::Functions::GridViewEntitySet EntitySet const& entitySet() const { return entitySet_; } public: - /// Return global basis + /// \brief Return global basis bound to the DOFVector GlobalBasis const& basis() const { return dofVector_->basis(); } - /// Return treePath associated with this view + /// \brief Return treePath associated with this discrete function TreePath const& treePath() const { return treePath_; } - /// Return const coefficient vector - DOFVector<GlobalBasis,RangeType> const& coefficients() const + /// \brief Return const coefficient vector + DOFVector<GB,VT> const& coefficients() const { return *dofVector_; } protected: - DOFVector<GlobalBasis,RangeType> const* dofVector_; - TreePath const treePath_; - + DOFVector<GB,VT> const* dofVector_; + TreePath treePath_; EntitySet entitySet_; NodeToRangeEntry nodeToRangeEntry_; }; + +#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION + // Deduction guide for DiscreteFunction class + template <class GlobalBasis, class ValueType> + DiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector) + -> DiscreteFunction<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>; +#endif + /// A Generator for a \ref DiscreteFunction - template <class GlobalBasis, class RangeType, class TreePath> - auto makeDiscreteFunction(DOFVector<GlobalBasis, RangeType> const& dofVector, TreePath const& treePath) + template <class GlobalBasis, class ValueType, class TreePath> + auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector, TreePath const& treePath) { - return DiscreteFunction<GlobalBasis, RangeType, TreePath>{dofVector, treePath}; + return DiscreteFunction<GlobalBasis, ValueType, TreePath>{dofVector, treePath}; } /// A Generator for a \ref DiscreteFunction - template <class GlobalBasis, class RangeType> - auto makeDiscreteFunction(DOFVector<GlobalBasis, RangeType> const& dofVector) + template <class GlobalBasis, class ValueType> + auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector) { auto treePath = Dune::TypeTree::hybridTreePath(); - return DiscreteFunction<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath}; + return DiscreteFunction<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath}; } } // end namespace AMDiS diff --git a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp index fc6faa734e72cd37685941a895faf11f7cd42c95..730e4b82d30dd55733eb438efa376077ff4a83c4 100644 --- a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp +++ b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp @@ -7,9 +7,140 @@ namespace AMDiS { -template <class GB, class RT, class TP> -typename DiscreteFunction<GB,RT,TP>::Range DiscreteFunction<GB,RT,TP>:: -LocalFunction::operator()(LocalDomain const& x) const +template <class GB, class VT, class TP> +class DiscreteFunction<GB,VT,TP>::LocalFunction +{ +public: + using Domain = typename EntitySet::LocalCoordinate; + using Range = typename DiscreteFunction::Range; + + enum { hasDerivative = true }; + +private: + using LocalView = typename GlobalBasis::LocalView; + using Element = typename EntitySet::Element; + using Geometry = typename Element::Geometry; + +public: + LocalFunction(DiscreteFunction const& globalFunction) + : globalFunction_(&globalFunction) + , localView_(globalFunction_->basis().localView()) + , subTree_(&child(localView_.tree(), globalFunction_->treePath())) + {} + + /// \brief Bind the LocalView to the element + void bind(Element const& element) + { + localView_.bind(element); + bound_ = true; + } + + /// \brief Unbind the LocalView from the element + void unbind() + { + localView_.unbind(); + bound_ = false; + } + + /// \brief Evaluate LocalFunction at bound element in local coordinates + Range operator()(Domain const& x) const; + + /// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction + friend GradientLocalFunction derivative(LocalFunction const& localFunction) + { + return GradientLocalFunction{*localFunction.globalFunction_}; + } + + /// \brief The \ref polynomialDegree() of the LocalFunctions + friend int order(LocalFunction const& self) + { + assert( self.bound_ ); + return polynomialDegree(*self.subTree_); + } + + /// \brief Return the bound element + Element const& localContext() const + { + assert( bound_ ); + return localView_.element(); + } + +private: + DiscreteFunction const* globalFunction_; + LocalView localView_; + SubTree const* subTree_; + bool bound_ = false; +}; + + +template <class GB, class VT, class TP> +class DiscreteFunction<GB,VT,TP>::GradientLocalFunction +{ + using R = typename DiscreteFunction::Range; + using D = typename DiscreteFunction::Domain; + using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature; + using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>; + +public: + using Domain = typename EntitySet::LocalCoordinate; + using Range = typename DerivativeTraits::Range; + + enum { hasDerivative = false }; + +private: + using LocalView = typename GlobalBasis::LocalView; + using Element = typename EntitySet::Element; + using Geometry = typename Element::Geometry; + +public: + GradientLocalFunction(DiscreteFunction const& globalFunction) + : globalFunction_(&globalFunction) + , localView_(globalFunction_->basis().localView()) + , subTree_(&child(localView_.tree(), globalFunction_->treePath())) + {} + + void bind(Element const& element) + { + localView_.bind(element); + geometry_.emplace(element.geometry()); + bound_ = true; + } + + void unbind() + { + localView_.unbind(); + geometry_.reset(); + bound_ = false; + } + + /// Evaluate Gradient at bound element in local coordinates + Range operator()(Domain const& x) const; + + friend int order(GradientLocalFunction const& self) + { + assert( self.bound_ ); + return std::max(0, polynomialDegree(*self.subTree_)-1); + } + + /// Return the bound element + Element const& localContext() const + { + assert( bound_ ); + return localView_.element(); + } + +private: + DiscreteFunction const* globalFunction_; + LocalView localView_; + SubTree const* subTree_; + Dune::Std::optional<Geometry> geometry_; + bool bound_ = false; +}; + + +template <class GB, class VT, class TP> +typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>:: +LocalFunction::operator()(Domain const& x) const { assert( bound_ ); auto y = Range(0); @@ -51,9 +182,9 @@ LocalFunction::operator()(LocalDomain const& x) const } -template <class GB, class RT, class TP> -typename DiscreteFunction<GB,RT,TP>::DerivativeRange DiscreteFunction<GB,RT,TP>:: -GradientLocalFunction::operator()(LocalDomain const& x) const +template <class GB, class VT, class TP> +typename DiscreteFunction<GB,VT,TP>::GradientLocalFunction::Range DiscreteFunction<GB,VT,TP>:: +GradientLocalFunction::operator()(Domain const& x) const { assert( bound_ ); Range dy; @@ -109,12 +240,12 @@ GradientLocalFunction::operator()(LocalDomain const& x) const } -template <class GB, class RT, class TP> -typename DiscreteFunction<GB,RT,TP>::Range DiscreteFunction<GB,RT,TP>:: +template <class GB, class VT, class TP> +typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>:: operator()(Domain const& x) const { - using Grid = typename GB::GridView::Grid; - using IS = typename GB::GridView::IndexSet; + using Grid = typename GlobalBasis::GridView::Grid; + using IS = typename GlobalBasis::GridView::IndexSet; auto const& gv = this->basis().gridView(); Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()}; diff --git a/src/amdis/utility/RangeType.hpp b/src/amdis/utility/RangeType.hpp index cee8d959e970b86769290983aa38c2797ec0d86d..27116104982776775ae39e9369c5d117623c9536 100644 --- a/src/amdis/utility/RangeType.hpp +++ b/src/amdis/utility/RangeType.hpp @@ -11,31 +11,43 @@ namespace AMDiS { namespace Impl { - template <class Node, bool isLeaf, bool isPower, bool isComposite> + template <class Node, class R, bool isLeaf, bool isPower, bool isComposite> struct RangeTypeImpl { static_assert( isLeaf || isPower || isComposite, "Unknown node-type for range definition" ); }; } - template <class Node> + /// \brief Range type of a node in the basis tree, composed of the leaf basis + /// range types. + /** + * Generate the range type by recursively combining leaf range types to a + * hybrid node range type. Range types for PowerNodes are thereby constructed + * as Dune::FieldVector of child range types. CompositeNodes produce a + * MultiTypeVector of the difference child range types. + * + * \tparam Node Type of a basis-tree node + * \tparam R Coefficient type [double] + **/ + template <class Node, class R = double> using RangeType_t = - typename Impl::RangeTypeImpl<Node, Node::isLeaf, Node::isPower, Node::isComposite>::type; + typename Impl::RangeTypeImpl<Node, R, Node::isLeaf, Node::isPower, Node::isComposite>::type; namespace Impl { // Leaf node - template <class Node> - struct RangeTypeImpl<Node, true, false, false> + template <class Node, class R> + struct RangeTypeImpl<Node, R, true, false, false> { using LocalBasis = typename Node::FiniteElement::Traits::LocalBasisType; - using type = typename LocalBasis::Traits::RangeType; + using T = typename LocalBasis::Traits::RangeType; + using type = std::decay_t<decltype(std::declval<R>() * std::declval<T>())>; }; // Power node - template <class Node> - struct RangeTypeImpl<Node, false, true, false> + template <class Node, class R> + struct RangeTypeImpl<Node, R, false, true, false> { using ChildNode = typename Node::template Child<0>::type; @@ -49,12 +61,12 @@ namespace AMDiS using type = Dune::FieldVector<T, int(Node::CHILDREN)>; }; - using type = typename FlatType<ChildNode::isLeaf, RangeType_t<ChildNode>>::type; + using type = typename FlatType<ChildNode::isLeaf, RangeType_t<ChildNode,R>>::type; }; // Composite node - template <class Node> - struct RangeTypeImpl<Node, false, false, true> + template <class Node, class R> + struct RangeTypeImpl<Node, R, false, false, true> { template <class Idx> struct RangeTypeGenerator; template <std::size_t... I> @@ -62,7 +74,7 @@ namespace AMDiS { template <std::size_t J> using ChildNode = typename Node::template Child<J>::type; - using type = MultiTypeVector<RangeType_t<ChildNode<I>>...>; + using type = MultiTypeVector<RangeType_t<ChildNode<I>,R>...>; }; using type = typename RangeTypeGenerator<MakeSeq_t<Node::CHILDREN>>::type; diff --git a/test/DiscreteFunctionTest.cpp b/test/DiscreteFunctionTest.cpp index bd3edbf45d235c9a79bcab458469f793600b8324..c2cf3e4f733c26302e798f0be67c06ba76a7e83b 100644 --- a/test/DiscreteFunctionTest.cpp +++ b/test/DiscreteFunctionTest.cpp @@ -2,6 +2,15 @@ // vi: set et ts=4 sw=2 sts=2: #include <iostream> +#include <type_traits> + +#include <quadmath.h> + +namespace std +{ + template <> + struct is_arithmetic<__float128> : true_type {}; +} #include <amdis/AMDiS.hpp> #include <amdis/ProblemStat.hpp> @@ -97,6 +106,10 @@ int main(int argc, char** argv) dlocalFct.unbind(); } + auto V0 = makeDOFVector<__float128>(prob.globalBasis()); + auto v0 = makeDOFVectorView(V0); + v0 << expr; + AMDiS::finalize(); return 0; }