From 0e947c5c64607a35b6bc398d6f6d3c8538a93873 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sat, 7 Dec 2019 23:56:26 +0100 Subject: [PATCH 01/10] renamed DOFVectorView into DiscreteFunction --- src/amdis/DOFVector.hpp | 8 +- src/amdis/ProblemStat.hpp | 1 - src/amdis/gridfunctions/CMakeLists.txt | 2 +- src/amdis/gridfunctions/DOFVectorView.hpp | 150 ------- src/amdis/gridfunctions/DiscreteFunction.hpp | 164 +++++-- .../gridfunctions/DiscreteFunction.inc.hpp | 416 ++---------------- .../DiscreteLocalFunction.inc.hpp | 383 ++++++++++++++++ src/amdis/io/FileWriterCreator.hpp | 4 +- src/amdis/io/VTKWriter.hpp | 8 +- 9 files changed, 562 insertions(+), 574 deletions(-) delete mode 100644 src/amdis/gridfunctions/DOFVectorView.hpp create mode 100644 src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp diff --git a/src/amdis/DOFVector.hpp b/src/amdis/DOFVector.hpp index 1230f216..4549edf9 100644 --- a/src/amdis/DOFVector.hpp +++ b/src/amdis/DOFVector.hpp @@ -88,7 +88,7 @@ namespace AMDiS auto child(TreePath const& path = {}) { auto&& tp = makeTreePath(path); - return makeDOFVectorView(*this, tp); + return makeDiscreteFunction(*this, tp); } template @@ -101,7 +101,7 @@ namespace AMDiS /// Interpolation of GridFunction to DOFVector, assuming that there is no /// reference to this DOFVector in the expression. - /// See \ref DOFVectorView::interpolate_noalias + /// See \ref DiscreteFunction::interpolate_noalias template void interpolate_noalias(Expr&& expr, Tag strategy) { @@ -109,7 +109,7 @@ namespace AMDiS } /// Interpolation of GridFunction to DOFVector. - /// See \ref DOFVectorView::interpolate + /// See \ref DiscreteFunction::interpolate template void interpolate(Expr&& expr, Tag strategy) { @@ -117,7 +117,7 @@ namespace AMDiS } /// Interpolation of GridFunction to DOFVector. - /// See \ref DOFVectorView::interpolate + /// See \ref DiscreteFunction::interpolate template DOFVector& operator<<(Expr&& expr) { diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp index b50009b1..bea802b2 100644 --- a/src/amdis/ProblemStat.hpp +++ b/src/amdis/ProblemStat.hpp @@ -40,7 +40,6 @@ #include #include -#include #include diff --git a/src/amdis/gridfunctions/CMakeLists.txt b/src/amdis/gridfunctions/CMakeLists.txt index ead5dfa0..e160da3c 100644 --- a/src/amdis/gridfunctions/CMakeLists.txt +++ b/src/amdis/gridfunctions/CMakeLists.txt @@ -8,7 +8,7 @@ install(FILES DerivativeGridFunction.hpp DiscreteFunction.hpp DiscreteFunction.inc.hpp - DOFVectorView.hpp + DiscreteLocalFunction.inc.hpp FunctorGridFunction.hpp GridFunction.hpp OperationsGridFunction.hpp diff --git a/src/amdis/gridfunctions/DOFVectorView.hpp b/src/amdis/gridfunctions/DOFVectorView.hpp deleted file mode 100644 index 99741d16..00000000 --- a/src/amdis/gridfunctions/DOFVectorView.hpp +++ /dev/null @@ -1,150 +0,0 @@ -#pragma once - -#include -#include -#include - -namespace AMDiS -{ - /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction - template - class DOFVectorView - : public DiscreteFunction - { - using Self = DOFVectorView; - using Super = DiscreteFunction; - - using GlobalBasis = GB; - using TreePath = TP; - - public: - /// Constructor. Stores a pointer to the mutable `dofvector`. - template - DOFVectorView(DOFVector& dofVector, PreTreePath const& preTreePath) - : Super(dofVector, makeTreePath(preTreePath)) - , mutableDofVector_(&dofVector) - {} - - /// Constructor forwards to the treePath constructor, with empty TreePath - DOFVectorView(DOFVector& dofVector) - : DOFVectorView(dofVector, Dune::TypeTree::hybridTreePath()) - {} - - public: - /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no - /// reference to this DOFVector in the expression. - /** - * **Example:** - * ``` - * auto v = makeDOFVectorView(prob.solutionVector(),0); - * v.interpolate_noalias([](auto const& x) { return x[0]; }); - * ``` - **/ - template - void interpolate_noalias(Expr&& expr, Tag strategy = {}) - { - auto const& basis = *this->basis(); - auto const& treePath = this->treePath(); - - auto&& gf = makeGridFunction(FWD(expr), basis.gridView()); - - if (std::is_same::value) { - auto counter = coefficients(); - AMDiS::interpolate(basis, coefficients(), gf, treePath, counter); - - coefficients().forEach([&counter](std::size_t dof, auto& coeff) - { - coeff /= std::max(double(counter.at(dof)), 1.0); - }); - } else { - AMDiS::interpolate(basis, coefficients(), gf, treePath); - } - } - - /// \brief Interpolation of GridFunction to DOFVector - /** - * **Example:** - * ``` - * auto v = makeDOFVectorView(prob.solutionVector(),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 - void interpolate(Expr&& expr, Tag strategy = {}) - { - // create temporary copy of data - DOFVector tmp(coefficients()); - - Self tmpView{tmp, this->treePath()}; - tmpView.interpolate_noalias(FWD(expr), strategy); - - // move data from temporary vector into stored DOFVector - coefficients().backend() = std::move(tmp.backend()); - } - - /// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate() - template - DOFVectorView& operator<<(Expr&& expr) - { - interpolate(FWD(expr)); - return *this; - } - - /// \brief interpolate `(*this) + expr` to DOFVector - template - DOFVectorView& operator+=(Expr&& expr) - { - interpolate((*this) + expr); - return *this; - } - - /// \brief interpolate `(*this) - expr` to DOFVector - template - DOFVectorView& operator-=(Expr&& expr) - { - interpolate((*this) - expr); - return *this; - } - - /// Return the mutable DOFVector - DOFVector& coefficients() { return *mutableDofVector_; } - - /// Return the const DOFVector - using Super::coefficients; - - protected: - DOFVector* mutableDofVector_; - }; - - -#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION - // Deduction guide for DOFVectorView class - template - DOFVectorView(DOFVector& dofVector, PreTreePath const& preTreePath) - -> DOFVectorView; - - // Deduction guide for DOFVectorView class - template - DOFVectorView(DOFVector& dofVector) - -> DOFVectorView>; -#endif - - /// A Generator for a mutable \ref DOFVectorView - template - auto makeDOFVectorView(DOFVector& dofVector, PreTreePath const& preTreePath) - { - auto treePath = makeTreePath(preTreePath); - return DOFVectorView{dofVector, treePath}; - } - - /// A Generator for a mutable \ref DOFVectorView - template - auto makeDOFVectorView(DOFVector& dofVector) - { - auto treePath = Dune::TypeTree::hybridTreePath(); - return DOFVectorView>{dofVector, treePath}; - } - -} // end namespace AMDiS diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp index 8a9ba447..f9de7bd9 100644 --- a/src/amdis/gridfunctions/DiscreteFunction.hpp +++ b/src/amdis/gridfunctions/DiscreteFunction.hpp @@ -19,20 +19,57 @@ namespace AMDiS /// \class DiscreteFunction /// \brief A view on a subspace of a \ref DOFVector /** - * \ingroup GridFunctions - * - * \tparam GB Type of the global basis - * \tparam VT Coefficient type of the DOFVector - * \tparam TP A realization of \ref Dune::TypeTree::HybridTreePath - * - * **Requirements:** - * - GB models \ref Dune::Functions::Concept::GlobalBasis - **/ - template - class DiscreteFunction + * \ingroup GridFunctions + * + * \tparam GB Type of the global basis + * \tparam VT Coefficient type of the DOFVector + * \tparam TP A realization of \ref Dune::TypeTree::HybridTreePath + * \tparam is_const Specifies whether a const or mutable view is implemented. + * + * **Requirements:** + * - GB models \ref Dune::Functions::Concept::GlobalBasis + **/ + template + class DiscreteFunction; + +#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION + // Deduction guide for DiscreteFunction class + template + DiscreteFunction(DOFVector const& dofVector) + -> DiscreteFunction, true>; + + template + DiscreteFunction(DOFVector& dofVector) + -> DiscreteFunction, false>; +#endif + + + /// A Generator for a const \ref DiscreteFunction + template > + auto makeDiscreteFunction(DOFVector const& dofVector, + PreTreePath const& preTreePath = {}) { - // static_assert(Dune::IsNumber::value, ""); + auto treePath = makeTreePath(preTreePath); + return DiscreteFunction{dofVector, treePath}; + } + /// A Generator for a mutable \ref DiscreteFunction + template > + auto makeDiscreteFunction(DOFVector& dofVector, + PreTreePath const& preTreePath = {}) + { + auto treePath = makeTreePath(preTreePath); + return DiscreteFunction{dofVector, treePath}; + } + + + + /// A Const DiscreteFunction + template + class DiscreteFunction + { private: using GlobalBasis = GB; using TreePath = TP; @@ -71,18 +108,13 @@ namespace AMDiS public: /// Constructor. Stores a pointer to the dofVector and a copy of the treePath. - DiscreteFunction(DOFVector const& dofVector, TP const& treePath) + DiscreteFunction(DOFVector const& dofVector, TP const& treePath = {}) : dofVector_(&dofVector) , treePath_(treePath) , entitySet_(dofVector.basis()->gridView()) , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(*dofVector.basis(), treePath)) {} - /// Constructor forwards to the treePath constructor, with empty TreePath - DiscreteFunction(DOFVector const& dofVector) - : DiscreteFunction(dofVector, Dune::TypeTree::hybridTreePath()) - {} - /// \brief Evaluate DiscreteFunction in global coordinates. NOTE: expensive Range operator()(Domain const& x) const; @@ -98,7 +130,6 @@ namespace AMDiS return entitySet_; } - public: /// \brief Return global basis bound to the DOFVector std::shared_ptr basis() const { @@ -125,29 +156,86 @@ namespace AMDiS }; -#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION - // Deduction guide for DiscreteFunction class - template - DiscreteFunction(DOFVector const& dofVector) - -> DiscreteFunction>; -#endif - /// A Generator for a \ref DiscreteFunction - template - auto makeDiscreteFunction(DOFVector const& dofVector, PreTreePath const& preTreePath) + /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction + template + class DiscreteFunction + : public DiscreteFunction { - auto treePath = makeTreePath(preTreePath); - return DiscreteFunction{dofVector, treePath}; - } + using Self = DiscreteFunction; + using Super = DiscreteFunction; - /// A Generator for a \ref DiscreteFunction - template - auto makeDiscreteFunction(DOFVector const& dofVector) - { - auto treePath = Dune::TypeTree::hybridTreePath(); - return DiscreteFunction>{dofVector, treePath}; - } + using GlobalBasis = GB; + using TreePath = TP; + + public: + /// Constructor. Stores a pointer to the mutable `dofvector`. + DiscreteFunction(DOFVector& dofVector, TP const& treePath = {}) + : Super(dofVector, treePath) + , mutableDofVector_(&dofVector) + {} + + public: + /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no + /// reference to this DOFVector in the expression. + /** + * **Example:** + * ``` + * auto v = makeDiscreteFunction(prob.solutionVector(),0); + * v.interpolate_noalias([](auto const& x) { return x[0]; }); + * ``` + **/ + template + void interpolate_noalias(Expr&& expr, Tag strategy = {}); + + /// \brief Interpolation of GridFunction to DOFVector + /** + * **Example:** + * ``` + * auto v = makeDiscreteFunction(prob.solutionVector(),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 + void interpolate(Expr&& expr, Tag strategy = {}); + + /// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate() + template + Self& operator<<(Expr&& expr) + { + interpolate(FWD(expr)); + return *this; + } + + /// \brief interpolate `(*this) + expr` to DOFVector + template + Self& operator+=(Expr&& expr) + { + interpolate((*this) + expr); + return *this; + } + + /// \brief interpolate `(*this) - expr` to DOFVector + template + Self& operator-=(Expr&& expr) + { + interpolate((*this) - expr); + return *this; + } + + /// Return the mutable DOFVector + DOFVector& coefficients() { return *mutableDofVector_; } + + /// Return the const DOFVector + using Super::coefficients; + + protected: + DOFVector* mutableDofVector_; + }; } // end namespace AMDiS +#include "DiscreteLocalFunction.inc.hpp" #include "DiscreteFunction.inc.hpp" diff --git a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp index 1a18c464..3759909c 100644 --- a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp +++ b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp @@ -1,406 +1,72 @@ #pragma once -#include -#include -#include -#include +#include -#include #include -namespace AMDiS { - -template -class DiscreteFunction::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: - /// Constructor. Stores a copy of the DiscreteFunction. - LocalFunction(DiscreteFunction const& globalFunction) - : globalFunction_(globalFunction) - , localView_(globalFunction_.basis()->localView()) - , subTree_(&child(localView_.tree(), globalFunction_.treePath())) - {} - - /// Copy constructor. - LocalFunction(LocalFunction const& other) - : globalFunction_(other.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); - - globalFunction_.coefficients().gather(localView_, localCoefficients_); - 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 - GradientLocalFunction makeDerivative(tag::gradient type) const - { - return GradientLocalFunction{globalFunction_, type}; - } - - DivergenceLocalFunction makeDerivative(tag::divergence type) const - { - return DivergenceLocalFunction{globalFunction_, type}; - } - - PartialLocalFunction makeDerivative(tag::partial type) const - { - return PartialLocalFunction{globalFunction_, type}; - } - - /// \brief The \ref polynomialDegree() of the LocalFunctions - int order() const - { - assert( bound_ ); - return polynomialDegree(*subTree_); - } - - /// \brief Return the bound element - Element const& localContext() const - { - assert( bound_ ); - return localView_.element(); - } - -private: - DiscreteFunction globalFunction_; - LocalView localView_; - SubTree const* subTree_; - - std::vector localCoefficients_; - bool bound_ = false; -}; +#include +#include +namespace AMDiS { +// Evaluate DiscreteFunction in global coordinates template -typename DiscreteFunction::Range DiscreteFunction:: -LocalFunction::operator()(Domain const& x) const +typename DiscreteFunction::Range DiscreteFunction:: + operator()(Domain const& x) const { - assert( bound_ ); - Range y(0); - - auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_; - for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp) - { - auto localBasisCache = makeNodeCache(node); - auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(), x); - std::size_t size = node.finiteElement().size(); - - // Get range entry associated to this node - auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y)); - - for (std::size_t i = 0; i < size; ++i) { - // Get coefficient associated to i-th shape function - auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]); - - // Get value of i-th shape function - auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]); + using Grid = typename GlobalBasis::GridView::Grid; + using IS = typename GlobalBasis::GridView::IndexSet; - std::size_t dimC = c.size(); - std::size_t dimV = v.size(); - assert(dimC*dimV == std::size_t(re.size())); - for(std::size_t j = 0; j < dimC; ++j) { - auto&& c_j = c[j]; - for(std::size_t k = 0; k < dimV; ++k) - re[j*dimV + k] += c_j*v[k]; - } - } - }); + auto const& gv = this->basis()->gridView(); + Dune::HierarchicSearch hsearch{gv.grid(), gv.indexSet()}; - return y; + auto element = hsearch.findEntity(x); + auto geometry = element.geometry(); + auto localFct = localFunction(*this); + localFct.bind(element); + return localFct(geometry.local(x)); } +// Interpolation of GridFunction to DOFVector template - template -class DiscreteFunction::DerivativeLocalFunctionBase -{ - using R = typename DiscreteFunction::Range; - using D = typename DiscreteFunction::Domain; - using RawSignature = typename Dune::Functions::SignatureTraits::RawSignature; - using Traits = DerivativeTraits; - -public: - using Domain = typename EntitySet::LocalCoordinate; - using Range = typename Traits::Range; - - enum { hasDerivative = false }; - -private: - using LocalView = typename GlobalBasis::LocalView; - using Element = typename EntitySet::Element; - using Geometry = typename Element::Geometry; - -public: - /// Constructor. Stores a copy of the DiscreteFunction. - DerivativeLocalFunctionBase(DiscreteFunction const& globalFunction, Type const& type) - : globalFunction_(globalFunction) - , type_(type) - , localView_(globalFunction_.basis()->localView()) - , subTree_(&child(localView_.tree(), globalFunction_.treePath())) - {} - - /// Copy constructor. - DerivativeLocalFunctionBase(DerivativeLocalFunctionBase const& other) - : globalFunction_(other.globalFunction_) - , type_(other.type_) - , localView_(globalFunction_.basis()->localView()) - , subTree_(&child(localView_.tree(), globalFunction_.treePath())) - {} - - void bind(Element const& element) - { - localView_.bind(element); - geometry_.emplace(element.geometry()); - - globalFunction_.coefficients().gather(localView_, localCoefficients_); - bound_ = true; - } - - void unbind() - { - localView_.unbind(); - geometry_.reset(); - bound_ = false; - } - - int order() const - { - assert( bound_ ); - return std::max(0, polynomialDegree(*subTree_)-1); - } - - /// Return the bound element - Element const& localContext() const - { - assert( bound_ ); - return localView_.element(); - } - - Geometry const& geometry() const - { - assert( bound_ ); - return geometry_.value(); - } - - LocalView const& localView() const - { - return localView_; - } - -protected: - DiscreteFunction globalFunction_; - Type type_; - LocalView localView_; - SubTree const* subTree_; - Dune::Std::optional geometry_; - std::vector localCoefficients_; - bool bound_ = false; -}; - - -template -class DiscreteFunction::GradientLocalFunction - : public DerivativeLocalFunctionBase -{ - using Super = DerivativeLocalFunctionBase; -public: - using Range = typename Super::Range; - using Domain = typename Super::Domain; - - // adopt constructor from base class - using DerivativeLocalFunctionBase::DerivativeLocalFunctionBase; - - /// Evaluate Gradient at bound element in local coordinates - Range operator()(Domain const& x) const - { - assert( this->bound_ ); - Range dy(0); - - auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_; - for_each_leaf_node(*this->subTree_, [&](auto const& node, auto const& tp) - { - auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry()); - auto const& gradients = localBasis.gradientsAt(x); - - // Get range entry associated to this node - auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy)); - - for (std::size_t i = 0; i < localBasis.size(); ++i) { - // Get coefficient associated to i-th shape function - auto c = Dune::Functions::flatVectorView(this->localCoefficients_[node.localIndex(i)]); - - // Get value of i-th transformed reference gradient - auto grad = Dune::Functions::flatVectorView(gradients[i]); - - std::size_t dimC = c.size(); - std::size_t dimV = grad.size(); - assert(dimC*dimV == std::size_t(re.size())); - for(std::size_t j = 0; j < dimC; ++j) { - auto&& c_j = c[j]; - for(std::size_t k = 0; k < dimV; ++k) - re[j*dimV + k] += c_j*grad[k]; - } - } - }); - - return dy; - } - - using Super::localCoefficients_; -}; - - -template -class DiscreteFunction::DivergenceLocalFunction - : public DerivativeLocalFunctionBase -{ - using Super = DerivativeLocalFunctionBase; - -public: - using Range = typename Super::Range; - using Domain = typename Super::Domain; - - // adopt constructor from base class - using DerivativeLocalFunctionBase::DerivativeLocalFunctionBase; - - /// Evaluate divergence at bound element in local coordinates - Range operator()(Domain const& x) const - { - return evaluate_node(x, bool_t{}); - } - -private: - Range evaluate_node(Domain const& x, std::false_type) const - { - error_exit("Cannot calculate divergence(node) where node is not a power node."); - return Range(0); - } - - Range evaluate_node(Domain const& x, std::true_type) const - { - static_assert(static_size_v == 1, "Implemented for scalar coefficients only."); - - assert( this->bound_ ); - Range dy(0); - - auto&& node = *this->subTree_; - - auto localBasis = makeLocalToGlobalBasisAdapter(node.child(0), this->geometry()); - auto const& gradients = localBasis.gradientsAt(x); - - auto re = Dune::Functions::flatVectorView(dy); - assert(int(re.size()) == 1); - for (std::size_t i = 0; i < localBasis.size(); ++i) { - auto grad = Dune::Functions::flatVectorView(gradients[i]); - - assert(int(grad.size()) == GridView::dimensionworld); - for (std::size_t j = 0; j < GridView::dimensionworld; ++j) - re[0] += localCoefficients_[node.child(j).localIndex(i)] * grad[j]; - } - - return dy; - } - - using Super::localCoefficients_; -}; - - -template -class DiscreteFunction::PartialLocalFunction - : public DerivativeLocalFunctionBase + template +void DiscreteFunction:: + interpolate_noalias(Expr&& expr, Tag strategy) { - using Super = DerivativeLocalFunctionBase; - -public: - using Range = typename Super::Range; - using Domain = typename Super::Domain; + auto const& basis = *this->basis(); + auto const& treePath = this->treePath(); - using DerivativeLocalFunctionBase::DerivativeLocalFunctionBase; + auto&& gf = makeGridFunction(FWD(expr), basis.gridView()); - /// Evaluate partial derivative at bound element in local coordinates - Range operator()(Domain const& x) const - { - assert( this->bound_ ); - Range dy(0); + if (std::is_same::value) { + auto counter = coefficients(); + AMDiS::interpolate(basis, coefficients(), gf, treePath, counter); - std::size_t comp = this->type_.comp; - - auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_; - for_each_leaf_node(*this->subTree_, [&](auto const& node, auto const& tp) + coefficients().forEach([&counter](std::size_t dof, auto& coeff) { - auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry()); - auto const& partial = localBasis.partialsAt(x, comp); - - // Get range entry associated to this node - auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy)); - - for (std::size_t i = 0; i < localBasis.size(); ++i) { - // Get coefficient associated to i-th shape function - auto c = Dune::Functions::flatVectorView(this->localCoefficients_[node.localIndex(i)]); - - // Get value of i-th transformed reference partial_derivative - auto d_comp = Dune::Functions::flatVectorView(partial[i]); - - std::size_t dimC = c.size(); - std::size_t dimV = d_comp.size(); - assert(dimC*dimV == std::size_t(re.size())); - for(std::size_t j = 0; j < dimC; ++j) { - auto&& c_j = c[j]; - for(std::size_t k = 0; k < dimV; ++k) - re[j*dimV + k] += c_j*d_comp[k]; - } - } + coeff /= std::max(double(counter.at(dof)), 1.0); }); - - return dy; + } else { + AMDiS::interpolate(basis, coefficients(), gf, treePath); } - - using Super::localCoefficients_; -}; +} +// Interpolation of GridFunction to DOFVector template -typename DiscreteFunction::Range DiscreteFunction:: -operator()(Domain const& x) const + template +void DiscreteFunction:: + interpolate(Expr&& expr, Tag strategy) { - using Grid = typename GlobalBasis::GridView::Grid; - using IS = typename GlobalBasis::GridView::IndexSet; + // create temporary copy of data + DOFVector tmp(coefficients()); - auto const& gv = this->basis()->gridView(); - Dune::HierarchicSearch hsearch{gv.grid(), gv.indexSet()}; + Self tmpView{tmp, this->treePath()}; + tmpView.interpolate_noalias(FWD(expr), strategy); - auto element = hsearch.findEntity(x); - auto geometry = element.geometry(); - auto localFct = localFunction(*this); - localFct.bind(element); - return localFct(geometry.local(x)); + // move data from temporary vector into stored DOFVector + coefficients().backend() = std::move(tmp.backend()); } } // end namespace AMDiS diff --git a/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp b/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp new file mode 100644 index 00000000..1439118e --- /dev/null +++ b/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp @@ -0,0 +1,383 @@ +#pragma once + +#include +#include +#include +#include + +#include + +namespace AMDiS { + +template +class DiscreteFunction::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: + /// Constructor. Stores a copy of the DiscreteFunction. + LocalFunction(DiscreteFunction const& globalFunction) + : globalFunction_(globalFunction) + , localView_(globalFunction_.basis()->localView()) + , subTree_(&child(localView_.tree(), globalFunction_.treePath())) + {} + + /// Copy constructor. + LocalFunction(LocalFunction const& other) + : globalFunction_(other.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); + + globalFunction_.coefficients().gather(localView_, localCoefficients_); + 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 + { + assert( bound_ ); + Range y(0); + + auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_; + for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp) + { + auto localBasisCache = makeNodeCache(node); + auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(), x); + std::size_t size = node.finiteElement().size(); + + // Get range entry associated to this node + auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y)); + + for (std::size_t i = 0; i < size; ++i) { + // Get coefficient associated to i-th shape function + auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]); + + // Get value of i-th shape function + auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]); + + std::size_t dimC = c.size(); + std::size_t dimV = v.size(); + assert(dimC*dimV == std::size_t(re.size())); + for(std::size_t j = 0; j < dimC; ++j) { + auto&& c_j = c[j]; + for(std::size_t k = 0; k < dimV; ++k) + re[j*dimV + k] += c_j*v[k]; + } + } + }); + + return y; + } + + /// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction + GradientLocalFunction makeDerivative(tag::gradient type) const + { + return GradientLocalFunction{globalFunction_, type}; + } + + DivergenceLocalFunction makeDerivative(tag::divergence type) const + { + return DivergenceLocalFunction{globalFunction_, type}; + } + + PartialLocalFunction makeDerivative(tag::partial type) const + { + return PartialLocalFunction{globalFunction_, type}; + } + + /// \brief The \ref polynomialDegree() of the LocalFunctions + int order() const + { + assert( bound_ ); + return polynomialDegree(*subTree_); + } + + /// \brief Return the bound element + Element const& localContext() const + { + assert( bound_ ); + return localView_.element(); + } + +private: + DiscreteFunction globalFunction_; + LocalView localView_; + SubTree const* subTree_; + + std::vector localCoefficients_; + bool bound_ = false; +}; + + +template + template +class DiscreteFunction::DerivativeLocalFunctionBase +{ + using R = typename DiscreteFunction::Range; + using D = typename DiscreteFunction::Domain; + using RawSignature = typename Dune::Functions::SignatureTraits::RawSignature; + using Traits = DerivativeTraits; + +public: + using Domain = typename EntitySet::LocalCoordinate; + using Range = typename Traits::Range; + + enum { hasDerivative = false }; + +private: + using LocalView = typename GlobalBasis::LocalView; + using Element = typename EntitySet::Element; + using Geometry = typename Element::Geometry; + +public: + /// Constructor. Stores a copy of the DiscreteFunction. + DerivativeLocalFunctionBase(DiscreteFunction const& globalFunction, Type const& type) + : globalFunction_(globalFunction) + , type_(type) + , localView_(globalFunction_.basis()->localView()) + , subTree_(&child(localView_.tree(), globalFunction_.treePath())) + {} + + /// Copy constructor. + DerivativeLocalFunctionBase(DerivativeLocalFunctionBase const& other) + : globalFunction_(other.globalFunction_) + , type_(other.type_) + , localView_(globalFunction_.basis()->localView()) + , subTree_(&child(localView_.tree(), globalFunction_.treePath())) + {} + + void bind(Element const& element) + { + localView_.bind(element); + geometry_.emplace(element.geometry()); + + globalFunction_.coefficients().gather(localView_, localCoefficients_); + bound_ = true; + } + + void unbind() + { + localView_.unbind(); + geometry_.reset(); + bound_ = false; + } + + int order() const + { + assert( bound_ ); + return std::max(0, polynomialDegree(*subTree_)-1); + } + + /// Return the bound element + Element const& localContext() const + { + assert( bound_ ); + return localView_.element(); + } + + Geometry const& geometry() const + { + assert( bound_ ); + return geometry_.value(); + } + + LocalView const& localView() const + { + return localView_; + } + +protected: + DiscreteFunction globalFunction_; + Type type_; + LocalView localView_; + SubTree const* subTree_; + Dune::Std::optional geometry_; + std::vector localCoefficients_; + bool bound_ = false; +}; + + +template +class DiscreteFunction::GradientLocalFunction + : public DerivativeLocalFunctionBase +{ + using Super = DerivativeLocalFunctionBase; +public: + using Range = typename Super::Range; + using Domain = typename Super::Domain; + + // adopt constructor from base class + using DerivativeLocalFunctionBase::DerivativeLocalFunctionBase; + + /// Evaluate Gradient at bound element in local coordinates + Range operator()(Domain const& x) const + { + assert( this->bound_ ); + Range dy(0); + + auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_; + for_each_leaf_node(*this->subTree_, [&](auto const& node, auto const& tp) + { + auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry()); + auto const& gradients = localBasis.gradientsAt(x); + + // Get range entry associated to this node + auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy)); + + for (std::size_t i = 0; i < localBasis.size(); ++i) { + // Get coefficient associated to i-th shape function + auto c = Dune::Functions::flatVectorView(this->localCoefficients_[node.localIndex(i)]); + + // Get value of i-th transformed reference gradient + auto grad = Dune::Functions::flatVectorView(gradients[i]); + + std::size_t dimC = c.size(); + std::size_t dimV = grad.size(); + assert(dimC*dimV == std::size_t(re.size())); + for(std::size_t j = 0; j < dimC; ++j) { + auto&& c_j = c[j]; + for(std::size_t k = 0; k < dimV; ++k) + re[j*dimV + k] += c_j*grad[k]; + } + } + }); + + return dy; + } + + using Super::localCoefficients_; +}; + + +template +class DiscreteFunction::DivergenceLocalFunction + : public DerivativeLocalFunctionBase +{ + using Super = DerivativeLocalFunctionBase; + +public: + using Range = typename Super::Range; + using Domain = typename Super::Domain; + + // adopt constructor from base class + using DerivativeLocalFunctionBase::DerivativeLocalFunctionBase; + + /// Evaluate divergence at bound element in local coordinates + Range operator()(Domain const& x) const + { + return evaluate_node(x, bool_t{}); + } + +private: + Range evaluate_node(Domain const& x, std::false_type) const + { + error_exit("Cannot calculate divergence(node) where node is not a power node."); + return Range(0); + } + + Range evaluate_node(Domain const& x, std::true_type) const + { + static_assert(static_size_v == 1, "Implemented for scalar coefficients only."); + + assert( this->bound_ ); + Range dy(0); + + auto&& node = *this->subTree_; + + auto localBasis = makeLocalToGlobalBasisAdapter(node.child(0), this->geometry()); + auto const& gradients = localBasis.gradientsAt(x); + + auto re = Dune::Functions::flatVectorView(dy); + assert(int(re.size()) == 1); + for (std::size_t i = 0; i < localBasis.size(); ++i) { + auto grad = Dune::Functions::flatVectorView(gradients[i]); + + assert(int(grad.size()) == GridView::dimensionworld); + for (std::size_t j = 0; j < GridView::dimensionworld; ++j) + re[0] += localCoefficients_[node.child(j).localIndex(i)] * grad[j]; + } + + return dy; + } + + using Super::localCoefficients_; +}; + + +template +class DiscreteFunction::PartialLocalFunction + : public DerivativeLocalFunctionBase +{ + using Super = DerivativeLocalFunctionBase; + +public: + using Range = typename Super::Range; + using Domain = typename Super::Domain; + + using DerivativeLocalFunctionBase::DerivativeLocalFunctionBase; + + /// Evaluate partial derivative at bound element in local coordinates + Range operator()(Domain const& x) const + { + assert( this->bound_ ); + Range dy(0); + + std::size_t comp = this->type_.comp; + + auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_; + for_each_leaf_node(*this->subTree_, [&](auto const& node, auto const& tp) + { + auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry()); + auto const& partial = localBasis.partialsAt(x, comp); + + // Get range entry associated to this node + auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy)); + + for (std::size_t i = 0; i < localBasis.size(); ++i) { + // Get coefficient associated to i-th shape function + auto c = Dune::Functions::flatVectorView(this->localCoefficients_[node.localIndex(i)]); + + // Get value of i-th transformed reference partial_derivative + auto d_comp = Dune::Functions::flatVectorView(partial[i]); + + std::size_t dimC = c.size(); + std::size_t dimV = d_comp.size(); + assert(dimC*dimV == std::size_t(re.size())); + for(std::size_t j = 0; j < dimC; ++j) { + auto&& c_j = c[j]; + for(std::size_t k = 0; k < dimV; ++k) + re[j*dimV + k] += c_j*d_comp[k]; + } + } + }); + + return dy; + } + + using Super::localCoefficients_; +}; + + +} // end namespace AMDiS diff --git a/src/amdis/io/FileWriterCreator.hpp b/src/amdis/io/FileWriterCreator.hpp index 1193f19b..fe46c730 100644 --- a/src/amdis/io/FileWriterCreator.hpp +++ b/src/amdis/io/FileWriterCreator.hpp @@ -53,7 +53,7 @@ namespace AMDiS private: template std::unique_ptr - create_impl(std::string type, std::string prefix, DiscreteFunction const& data, ValCat) const + create_impl(std::string type, std::string prefix, DiscreteFunction const& data, ValCat) const { // ParaView VTK format, writer from dune-grid if (type == "vtk") @@ -92,7 +92,7 @@ namespace AMDiS // The value-category is unknown, like a composite/hierarchic vector or any unknown type. template std::unique_ptr - create_impl(std::string type, std::string prefix, DiscreteFunction const& /*data*/, tag::unknown) const + create_impl(std::string type, std::string prefix, DiscreteFunction const& /*data*/, tag::unknown) const { // Backup writer, writing the grid and the solution vector if (type == "backup") diff --git a/src/amdis/io/VTKWriter.hpp b/src/amdis/io/VTKWriter.hpp index c9b6bbab..37074806 100644 --- a/src/amdis/io/VTKWriter.hpp +++ b/src/amdis/io/VTKWriter.hpp @@ -42,7 +42,9 @@ namespace AMDiS using GridView = typename GB::GridView; using Writer = Dune::VTKWriter; using SeqWriter = VTKSequenceWriter; - using Range = typename DiscreteFunction::Range; + + using Function = DiscreteFunction; + using Range = typename Function::Range; template static constexpr Dune::VTK::FieldInfo::Type VTKFieldType() { @@ -56,7 +58,7 @@ namespace AMDiS public: /// Constructor. - VTKWriter(std::string const& name, DiscreteFunction const& discreteFct) + VTKWriter(std::string const& name, Function const& discreteFct) : FileWriterBase(name) , discreteFct_(discreteFct) { @@ -103,7 +105,7 @@ namespace AMDiS } private: - DiscreteFunction discreteFct_; + Function discreteFct_; std::shared_ptr vtkWriter_; std::shared_ptr vtkSeqWriter_; -- GitLab From 13a9fc32c22a04a1a746e4e8e740ba4efa9269f0 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sun, 8 Dec 2019 11:22:35 +0100 Subject: [PATCH 02/10] renamed makeDiscreteFunction to discreteFunction --- src/amdis/DOFVector.hpp | 4 +- src/amdis/gridfunctions/DiscreteFunction.hpp | 12 +++--- src/amdis/io/FileWriterCreator.hpp | 2 +- test/DiscreteFunctionTest.cpp | 42 ++++++++++++-------- 4 files changed, 35 insertions(+), 25 deletions(-) diff --git a/src/amdis/DOFVector.hpp b/src/amdis/DOFVector.hpp index 4549edf9..99c257a8 100644 --- a/src/amdis/DOFVector.hpp +++ b/src/amdis/DOFVector.hpp @@ -88,14 +88,14 @@ namespace AMDiS auto child(TreePath const& path = {}) { auto&& tp = makeTreePath(path); - return makeDiscreteFunction(*this, tp); + return discreteFunction(*this, tp); } template auto child(TreePath const& path = {}) const { auto&& tp = makeTreePath(path); - return makeDiscreteFunction(*this, tp); + return discreteFunction(*this, tp); } diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp index f9de7bd9..92da26c8 100644 --- a/src/amdis/gridfunctions/DiscreteFunction.hpp +++ b/src/amdis/gridfunctions/DiscreteFunction.hpp @@ -47,8 +47,8 @@ namespace AMDiS /// A Generator for a const \ref DiscreteFunction template > - auto makeDiscreteFunction(DOFVector const& dofVector, - PreTreePath const& preTreePath = {}) + auto discreteFunction(DOFVector const& dofVector, + PreTreePath const& preTreePath = {}) { auto treePath = makeTreePath(preTreePath); return DiscreteFunction{dofVector, treePath}; @@ -57,8 +57,8 @@ namespace AMDiS /// A Generator for a mutable \ref DiscreteFunction template > - auto makeDiscreteFunction(DOFVector& dofVector, - PreTreePath const& preTreePath = {}) + auto discreteFunction(DOFVector& dofVector, + PreTreePath const& preTreePath = {}) { auto treePath = makeTreePath(preTreePath); return DiscreteFunction{dofVector, treePath}; @@ -181,7 +181,7 @@ namespace AMDiS /** * **Example:** * ``` - * auto v = makeDiscreteFunction(prob.solutionVector(),0); + * auto v = discreteFunction(prob.solutionVector(),0); * v.interpolate_noalias([](auto const& x) { return x[0]; }); * ``` **/ @@ -192,7 +192,7 @@ namespace AMDiS /** * **Example:** * ``` - * auto v = makeDiscreteFunction(prob.solutionVector(),0); + * auto v = discreteFunction(prob.solutionVector(),0); * v.interpolate(v + [](auto const& x) { return x[0]; }); * ``` * Allows to have a reference to the DOFVector in the expression, e.g. as diff --git a/src/amdis/io/FileWriterCreator.hpp b/src/amdis/io/FileWriterCreator.hpp index fe46c730..8bb50ea9 100644 --- a/src/amdis/io/FileWriterCreator.hpp +++ b/src/amdis/io/FileWriterCreator.hpp @@ -44,7 +44,7 @@ namespace AMDiS std::unique_ptr create(std::string type, std::string prefix, TreePath treePath = {}) const { - auto data = makeDiscreteFunction(*systemVector_, treePath); + auto data = discreteFunction(*systemVector_, treePath); using Range = typename TYPEOF(data)::Range; return create_impl(std::move(type), std::move(prefix), data, ValueCategory_t{}); diff --git a/test/DiscreteFunctionTest.cpp b/test/DiscreteFunctionTest.cpp index f3ad5286..5665630a 100644 --- a/test/DiscreteFunctionTest.cpp +++ b/test/DiscreteFunctionTest.cpp @@ -7,7 +7,6 @@ #include #include #include -#include #include #include "Tests.hpp" @@ -55,9 +54,20 @@ int main(int argc, char** argv) auto U1 = *prob.solutionVector(); auto U2 = *prob.solutionVector(); - auto u0 = makeDOFVectorView(U0); - auto u1 = makeDOFVectorView(U1); - auto u2 = makeDOFVectorView(U2); + auto u0 = discreteFunction(U0); + auto u1 = discreteFunction(U1); + auto u2 = discreteFunction(U2); + + // Test const and mutable construction + auto const& U_c = *prob.solutionVector(); + auto& U_m = *prob.solutionVector(); + +#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION + DiscreteFunction u0_c(U_c); + DiscreteFunction u0_m(U_m); +#endif + auto u1_c = discreteFunction(U_c); + auto u1_m = discreteFunction(U_m); using Range = typename decltype(u0)::Range; @@ -123,30 +133,30 @@ int main(int argc, char** argv) } auto V0 = makeDOFVector(prob.globalBasis()); - auto v0 = makeDOFVectorView(V0); + auto v0 = discreteFunction(V0); v0 << expr1; - // test makeDiscreteFunction + // test discreteFunction int preTP1 = 0; std::integral_constant preTP2; auto tp = treepath(preTP1); auto W0 = *prob.solutionVector(); auto W1 = *prob.solutionVector(); auto W2 = *prob.solutionVector(); - auto w0 = makeDiscreteFunction(W0, preTP1); - auto w1 = makeDiscreteFunction(W1, preTP2); - auto w2 = makeDiscreteFunction(W2, tp); + auto w0 = discreteFunction(W0, preTP1); + auto w1 = discreteFunction(W1, preTP2); + auto w2 = discreteFunction(W2, tp); - // test makeDOFVectorView with (pre)treepath argument + // test discreteFunction with (pre)treepath argument auto expr = [](auto const& x) { return 1 + x[0] + x[1]; }; auto W3 = *prob.solutionVector(); auto W4 = *prob.solutionVector(); auto W5 = *prob.solutionVector(); auto W7 = *prob.solutionVector(); auto W8 = *prob.solutionVector(); - auto w3 = makeDOFVectorView(W3, preTP1); - auto w4 = makeDOFVectorView(W4, preTP2); - auto w5 = makeDOFVectorView(W5, tp); + auto w3 = discreteFunction(W3, preTP1); + auto w4 = discreteFunction(W4, preTP2); + auto w5 = discreteFunction(W5, tp); auto w6 = prob.solution(tp); auto& W6 = *prob.solutionVector(); w3 << expr; @@ -158,9 +168,9 @@ int main(int argc, char** argv) AMDIS_TEST( comp(W3, W6) ); // test interpolation on subbasis - auto w7 = makeDOFVectorView(W7); - auto w8_0 = makeDOFVectorView(W8, 0); - auto w8_1 = makeDOFVectorView(W8, 1); + auto w7 = discreteFunction(W7); + auto w8_0 = discreteFunction(W8, 0); + auto w8_1 = discreteFunction(W8, 1); w7 << expr; w8_0 << expr; w8_1 << expr; -- GitLab From fe4b5fefae11b7566b08f1d349dd5ae7e773a481 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sun, 8 Dec 2019 12:03:35 +0100 Subject: [PATCH 03/10] fixed renamed function makeDOFVectorView to discreteFunction --- test/GradientTest.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/GradientTest.cpp b/test/GradientTest.cpp index 85023b41..015dc42d 100644 --- a/test/GradientTest.cpp +++ b/test/GradientTest.cpp @@ -28,9 +28,9 @@ void test() auto uVector = makeDOFVector(basis); - auto u = makeDOFVectorView(uVector); - auto u_0 = makeDOFVectorView(uVector, 0); - auto u_1 = makeDOFVectorView(uVector, 1); + auto u = discreteFunction(uVector); + auto u_0 = discreteFunction(uVector, 0); + auto u_1 = discreteFunction(uVector, 1); // eval a functor at global coordinates (at quadrature points) u_0 << evalAtQP([](auto const& x) { return x[0] + x[1]; }); @@ -65,9 +65,9 @@ void test() (integrate(partialAtQP(u_0,0) + partialAtQP(u_1,1), gridView))); auto vVector(uVector); - auto v = makeDOFVectorView(vVector); - auto v_0 = makeDOFVectorView(vVector, 0); - auto v_1 = makeDOFVectorView(vVector, 1); + auto v = discreteFunction(vVector); + auto v_0 = discreteFunction(vVector, 0); + auto v_1 = discreteFunction(vVector, 1); // test gradient v << evalAtQP([](auto const& x) { -- GitLab From 3259559cea8779bc86414efa983544e7d7d5d0d1 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Fri, 20 Dec 2019 22:36:52 +0100 Subject: [PATCH 04/10] add missing includes --- src/amdis/DOFVector.hpp | 2 +- src/amdis/io/VTKWriter.hpp | 1 + test/GradientTest.cpp | 1 + 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/amdis/DOFVector.hpp b/src/amdis/DOFVector.hpp index 99c257a8..d3b9c5b2 100644 --- a/src/amdis/DOFVector.hpp +++ b/src/amdis/DOFVector.hpp @@ -10,7 +10,7 @@ #include #include #include -#include +#include #include namespace Dune diff --git a/src/amdis/io/VTKWriter.hpp b/src/amdis/io/VTKWriter.hpp index 37074806..b85ff48c 100644 --- a/src/amdis/io/VTKWriter.hpp +++ b/src/amdis/io/VTKWriter.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include diff --git a/test/GradientTest.cpp b/test/GradientTest.cpp index 015dc42d..f5e566d5 100644 --- a/test/GradientTest.cpp +++ b/test/GradientTest.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include "Tests.hpp" -- GitLab From c2a6cdf0b3fbb0a3e09ba5bbc09ad2cb584f4d88 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Fri, 20 Dec 2019 22:55:19 +0100 Subject: [PATCH 05/10] added forward declaration of DOFVector --- src/amdis/gridfunctions/DiscreteFunction.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp index 92da26c8..5abcaee3 100644 --- a/src/amdis/gridfunctions/DiscreteFunction.hpp +++ b/src/amdis/gridfunctions/DiscreteFunction.hpp @@ -16,6 +16,9 @@ namespace AMDiS { + template + class DOFVector; + /// \class DiscreteFunction /// \brief A view on a subspace of a \ref DOFVector /** -- GitLab From c81efe926e1fa34088ebaaf888e06c99ee2f0b89 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Thu, 26 Dec 2019 00:21:10 +0100 Subject: [PATCH 06/10] renamed maker function discreteFunction() back to makeDiscreteFunction() for consistency --- src/amdis/DOFVector.hpp | 4 +-- src/amdis/gridfunctions/DiscreteFunction.hpp | 12 ++++---- src/amdis/io/FileWriterCreator.hpp | 2 +- test/DiscreteFunctionTest.cpp | 30 ++++++++++---------- test/GradientTest.cpp | 12 ++++---- 5 files changed, 30 insertions(+), 30 deletions(-) diff --git a/src/amdis/DOFVector.hpp b/src/amdis/DOFVector.hpp index d3b9c5b2..f70b4d91 100644 --- a/src/amdis/DOFVector.hpp +++ b/src/amdis/DOFVector.hpp @@ -88,14 +88,14 @@ namespace AMDiS auto child(TreePath const& path = {}) { auto&& tp = makeTreePath(path); - return discreteFunction(*this, tp); + return makeDiscreteFunction(*this, tp); } template auto child(TreePath const& path = {}) const { auto&& tp = makeTreePath(path); - return discreteFunction(*this, tp); + return makeDiscreteFunction(*this, tp); } diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp index 5abcaee3..2c9d610c 100644 --- a/src/amdis/gridfunctions/DiscreteFunction.hpp +++ b/src/amdis/gridfunctions/DiscreteFunction.hpp @@ -50,8 +50,8 @@ namespace AMDiS /// A Generator for a const \ref DiscreteFunction template > - auto discreteFunction(DOFVector const& dofVector, - PreTreePath const& preTreePath = {}) + auto makeDiscreteFunction(DOFVector const& dofVector, + PreTreePath const& preTreePath = {}) { auto treePath = makeTreePath(preTreePath); return DiscreteFunction{dofVector, treePath}; @@ -60,8 +60,8 @@ namespace AMDiS /// A Generator for a mutable \ref DiscreteFunction template > - auto discreteFunction(DOFVector& dofVector, - PreTreePath const& preTreePath = {}) + auto makeDiscreteFunction(DOFVector& dofVector, + PreTreePath const& preTreePath = {}) { auto treePath = makeTreePath(preTreePath); return DiscreteFunction{dofVector, treePath}; @@ -184,7 +184,7 @@ namespace AMDiS /** * **Example:** * ``` - * auto v = discreteFunction(prob.solutionVector(),0); + * auto v = makeDiscreteFunction(prob.solutionVector(),0); * v.interpolate_noalias([](auto const& x) { return x[0]; }); * ``` **/ @@ -195,7 +195,7 @@ namespace AMDiS /** * **Example:** * ``` - * auto v = discreteFunction(prob.solutionVector(),0); + * auto v = makeDiscreteFunction(prob.solutionVector(),0); * v.interpolate(v + [](auto const& x) { return x[0]; }); * ``` * Allows to have a reference to the DOFVector in the expression, e.g. as diff --git a/src/amdis/io/FileWriterCreator.hpp b/src/amdis/io/FileWriterCreator.hpp index 8bb50ea9..fe46c730 100644 --- a/src/amdis/io/FileWriterCreator.hpp +++ b/src/amdis/io/FileWriterCreator.hpp @@ -44,7 +44,7 @@ namespace AMDiS std::unique_ptr create(std::string type, std::string prefix, TreePath treePath = {}) const { - auto data = discreteFunction(*systemVector_, treePath); + auto data = makeDiscreteFunction(*systemVector_, treePath); using Range = typename TYPEOF(data)::Range; return create_impl(std::move(type), std::move(prefix), data, ValueCategory_t{}); diff --git a/test/DiscreteFunctionTest.cpp b/test/DiscreteFunctionTest.cpp index 5665630a..96a947b4 100644 --- a/test/DiscreteFunctionTest.cpp +++ b/test/DiscreteFunctionTest.cpp @@ -54,9 +54,9 @@ int main(int argc, char** argv) auto U1 = *prob.solutionVector(); auto U2 = *prob.solutionVector(); - auto u0 = discreteFunction(U0); - auto u1 = discreteFunction(U1); - auto u2 = discreteFunction(U2); + auto u0 = makeDiscreteFunction(U0); + auto u1 = makeDiscreteFunction(U1); + auto u2 = makeDiscreteFunction(U2); // Test const and mutable construction auto const& U_c = *prob.solutionVector(); @@ -66,8 +66,8 @@ int main(int argc, char** argv) DiscreteFunction u0_c(U_c); DiscreteFunction u0_m(U_m); #endif - auto u1_c = discreteFunction(U_c); - auto u1_m = discreteFunction(U_m); + auto u1_c = makeDiscreteFunction(U_c); + auto u1_m = makeDiscreteFunction(U_m); using Range = typename decltype(u0)::Range; @@ -133,7 +133,7 @@ int main(int argc, char** argv) } auto V0 = makeDOFVector(prob.globalBasis()); - auto v0 = discreteFunction(V0); + auto v0 = makeDiscreteFunction(V0); v0 << expr1; // test discreteFunction @@ -143,9 +143,9 @@ int main(int argc, char** argv) auto W0 = *prob.solutionVector(); auto W1 = *prob.solutionVector(); auto W2 = *prob.solutionVector(); - auto w0 = discreteFunction(W0, preTP1); - auto w1 = discreteFunction(W1, preTP2); - auto w2 = discreteFunction(W2, tp); + auto w0 = makeDiscreteFunction(W0, preTP1); + auto w1 = makeDiscreteFunction(W1, preTP2); + auto w2 = makeDiscreteFunction(W2, tp); // test discreteFunction with (pre)treepath argument auto expr = [](auto const& x) { return 1 + x[0] + x[1]; }; @@ -154,9 +154,9 @@ int main(int argc, char** argv) auto W5 = *prob.solutionVector(); auto W7 = *prob.solutionVector(); auto W8 = *prob.solutionVector(); - auto w3 = discreteFunction(W3, preTP1); - auto w4 = discreteFunction(W4, preTP2); - auto w5 = discreteFunction(W5, tp); + auto w3 = makeDiscreteFunction(W3, preTP1); + auto w4 = makeDiscreteFunction(W4, preTP2); + auto w5 = makeDiscreteFunction(W5, tp); auto w6 = prob.solution(tp); auto& W6 = *prob.solutionVector(); w3 << expr; @@ -168,9 +168,9 @@ int main(int argc, char** argv) AMDIS_TEST( comp(W3, W6) ); // test interpolation on subbasis - auto w7 = discreteFunction(W7); - auto w8_0 = discreteFunction(W8, 0); - auto w8_1 = discreteFunction(W8, 1); + auto w7 = makeDiscreteFunction(W7); + auto w8_0 = makeDiscreteFunction(W8, 0); + auto w8_1 = makeDiscreteFunction(W8, 1); w7 << expr; w8_0 << expr; w8_1 << expr; diff --git a/test/GradientTest.cpp b/test/GradientTest.cpp index f5e566d5..23a7005f 100644 --- a/test/GradientTest.cpp +++ b/test/GradientTest.cpp @@ -29,9 +29,9 @@ void test() auto uVector = makeDOFVector(basis); - auto u = discreteFunction(uVector); - auto u_0 = discreteFunction(uVector, 0); - auto u_1 = discreteFunction(uVector, 1); + auto u = makeDiscreteFunction(uVector); + auto u_0 = makeDiscreteFunction(uVector, 0); + auto u_1 = makeDiscreteFunction(uVector, 1); // eval a functor at global coordinates (at quadrature points) u_0 << evalAtQP([](auto const& x) { return x[0] + x[1]; }); @@ -66,9 +66,9 @@ void test() (integrate(partialAtQP(u_0,0) + partialAtQP(u_1,1), gridView))); auto vVector(uVector); - auto v = discreteFunction(vVector); - auto v_0 = discreteFunction(vVector, 0); - auto v_1 = discreteFunction(vVector, 1); + auto v = makeDiscreteFunction(vVector); + auto v_0 = makeDiscreteFunction(vVector, 0); + auto v_1 = makeDiscreteFunction(vVector, 1); // test gradient v << evalAtQP([](auto const& x) { -- GitLab From 376d791ddf80052d21c27a5bd6334bc38c638741 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Thu, 26 Dec 2019 10:20:56 +0100 Subject: [PATCH 07/10] Add child() function to DiscreteFunction --- src/amdis/gridfunctions/DiscreteFunction.hpp | 16 ++++++++++++++++ .../gridfunctions/DiscreteLocalFunction.inc.hpp | 8 ++++---- src/amdis/io/DuneVtkWriter.hpp | 5 +++-- test/DiscreteFunctionTest.cpp | 7 +++++++ 4 files changed, 30 insertions(+), 6 deletions(-) diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp index 2c9d610c..ea97c780 100644 --- a/src/amdis/gridfunctions/DiscreteFunction.hpp +++ b/src/amdis/gridfunctions/DiscreteFunction.hpp @@ -151,6 +151,13 @@ namespace AMDiS return *dofVector_; } + template + auto child(TreePath const& path = {}) const + { + auto&& tp = makeTreePath(path); + return makeDiscreteFunction(*dofVector_, cat(this->treePath_,tp)); + } + protected: DOFVector const* dofVector_; TreePath treePath_; @@ -234,6 +241,15 @@ namespace AMDiS /// Return the const DOFVector using Super::coefficients; + template + auto child(TreePath const& path = {}) + { + auto&& tp = makeTreePath(path); + return makeDiscreteFunction(*mutableDofVector_, cat(this->treePath_,tp)); + } + + using Super::child; + protected: DOFVector* mutableDofVector_; }; diff --git a/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp b/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp index 1439118e..899818c6 100644 --- a/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp +++ b/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp @@ -28,14 +28,14 @@ public: LocalFunction(DiscreteFunction const& globalFunction) : globalFunction_(globalFunction) , localView_(globalFunction_.basis()->localView()) - , subTree_(&child(localView_.tree(), globalFunction_.treePath())) + , subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_.treePath())) {} /// Copy constructor. LocalFunction(LocalFunction const& other) : globalFunction_(other.globalFunction_) , localView_(globalFunction_.basis()->localView()) - , subTree_(&child(localView_.tree(), globalFunction_.treePath())) + , subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_.treePath())) {} /// \brief Bind the LocalView to the element @@ -157,7 +157,7 @@ public: : globalFunction_(globalFunction) , type_(type) , localView_(globalFunction_.basis()->localView()) - , subTree_(&child(localView_.tree(), globalFunction_.treePath())) + , subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_.treePath())) {} /// Copy constructor. @@ -165,7 +165,7 @@ public: : globalFunction_(other.globalFunction_) , type_(other.type_) , localView_(globalFunction_.basis()->localView()) - , subTree_(&child(localView_.tree(), globalFunction_.treePath())) + , subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_.treePath())) {} void bind(Element const& element) diff --git a/src/amdis/io/DuneVtkWriter.hpp b/src/amdis/io/DuneVtkWriter.hpp index 60b9fd8a..e6164c35 100644 --- a/src/amdis/io/DuneVtkWriter.hpp +++ b/src/amdis/io/DuneVtkWriter.hpp @@ -22,10 +22,11 @@ namespace AMDiS using GridView = typename GB::GridView; using Writer = Dune::VtkWriter; using SeqWriter = Dune::PvdWriter; + using Function = DiscreteFunction; public: /// Constructor. - DuneVtkWriter(std::string const& name, DiscreteFunction const& discreteFct) + DuneVtkWriter(std::string const& name, Function const& discreteFct) : FileWriterBase(name) , discreteFct_(discreteFct) { @@ -75,7 +76,7 @@ namespace AMDiS } private: - DiscreteFunction discreteFct_; + Function discreteFct_; std::shared_ptr vtkWriter_; std::shared_ptr vtkSeqWriter_; diff --git a/test/DiscreteFunctionTest.cpp b/test/DiscreteFunctionTest.cpp index 96a947b4..58d51d00 100644 --- a/test/DiscreteFunctionTest.cpp +++ b/test/DiscreteFunctionTest.cpp @@ -69,6 +69,13 @@ int main(int argc, char** argv) auto u1_c = makeDiscreteFunction(U_c); auto u1_m = makeDiscreteFunction(U_m); + // sub-range view on DiscreteFunction + auto su1_c = u1_c.child(); + auto su1_c0 = u1_c.child(0); + + auto su1_m = u1_m.child(); + auto su1_m0 = u1_m.child(0); + using Range = typename decltype(u0)::Range; auto expr1 = [](auto const& x) { -- GitLab From 0d138ebe18940198f7d86058669c92ebef518b92 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sat, 28 Dec 2019 11:00:22 +0100 Subject: [PATCH 08/10] renamed HAVE_XXX into AMDIS_HAS_XXX and updated Hypre cmake configuration --- src/amdis/LinearAlgebra.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/amdis/LinearAlgebra.hpp b/src/amdis/LinearAlgebra.hpp index f8c5d4de..d1dda9f1 100644 --- a/src/amdis/LinearAlgebra.hpp +++ b/src/amdis/LinearAlgebra.hpp @@ -37,7 +37,7 @@ #endif -#include -#include +#include +#include #include #include -- GitLab From aba164a7803477ed85e18c2d37215d2287a6d18b Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sat, 28 Dec 2019 11:42:54 +0100 Subject: [PATCH 09/10] update facade name --- src/amdis/LinearAlgebra.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/amdis/LinearAlgebra.hpp b/src/amdis/LinearAlgebra.hpp index d1dda9f1..f8c5d4de 100644 --- a/src/amdis/LinearAlgebra.hpp +++ b/src/amdis/LinearAlgebra.hpp @@ -37,7 +37,7 @@ #endif -#include -#include +#include +#include #include #include -- GitLab From ad922b1ba9e90c0fad149cfd21565c11bcf69c0d Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sat, 28 Dec 2019 14:59:43 +0100 Subject: [PATCH 10/10] changed config parameter in Environment --- src/amdis/Environment.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/amdis/Environment.cpp b/src/amdis/Environment.cpp index c7625130..c024e3a4 100644 --- a/src/amdis/Environment.cpp +++ b/src/amdis/Environment.cpp @@ -4,7 +4,7 @@ #include "Environment.hpp" -#if HAVE_PETSC +#if AMDIS_HAS_PETSC #include #endif @@ -25,7 +25,7 @@ namespace AMDiS Environment::Environment(int& argc, char**& argv, std::string const& initFileName) : Environment(initFileName) { -#if HAVE_PETSC +#if AMDIS_HAS_PETSC PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL); petscInitialized_ = true; #endif @@ -45,7 +45,7 @@ namespace AMDiS Environment::~Environment() { -#if HAVE_PETSC +#if AMDIS_HAS_PETSC if (petscInitialized_) PetscFinalize(); #endif -- GitLab