diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5537eb4b20e63bcfa30b0c60eadc6fb2b7219ec6..a42deff0475c4d050155499389d209f37db16180 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,13 @@
# Master
+- Replaced the class `GlobalGeodesicFEFunction` by a new one called
+ `GlobalGFEFunction`. There are two important changes: First of all,
+ the new class implements the `dune-functions` interface rather than
+ the deprecated one based on inheritance from `VirtualGridViewFunction`.
+ Secondly, the new class does not hard-wire geodesic interpolation
+ anymore. Rather, you can give it interpolation classes which then
+ govern how interpolation is done.
+
- Added dune-gmsh4 as a dependency
- Build cosserat-continuum for different combinations of LFE-orders and
GFE-orders, the respective program is called
diff --git a/dune/gfe/globalgeodesicfefunction.hh b/dune/gfe/globalgeodesicfefunction.hh
deleted file mode 100644
index 1ed94f1aa4f12fbbed8e64d69c2a960b5e24de30..0000000000000000000000000000000000000000
--- a/dune/gfe/globalgeodesicfefunction.hh
+++ /dev/null
@@ -1,121 +0,0 @@
-#ifndef GLOBAL_GEODESIC_FINITE_ELEMENT_FUNCTION_HH
-#define GLOBAL_GEODESIC_FINITE_ELEMENT_FUNCTION_HH
-
-#include <vector>
-
-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-
-#include <dune/gfe/localgeodesicfefunction.hh>
-
-template <class dctype, int DimDomain, class rctype, int DimRange>
-struct DerivativeTypefier<Dune::FieldVector<dctype, DimDomain>, UnitVector<rctype,DimRange> >
-{
- typedef Dune::FieldMatrix<rctype, DimRange, DimDomain> DerivativeType;
-};
-
-template <class dctype, int DimDomain, class rctype, int DimRange>
-struct DerivativeTypefier<Dune::FieldVector<dctype, DimDomain>, RealTuple<rctype,DimRange> >
-{
- typedef Dune::FieldMatrix<rctype, DimRange, DimDomain> DerivativeType;
-};
-
-
-
-/** \brief Global geodesic finite element function.
- *
- * \tparam B - The global basis type.
- * \tparam TargetSpace - The manifold that this functions takes its values in.
- */
-template<class B, class TargetSpace>
-class GlobalGeodesicFEFunction
-: public VirtualGridViewFunction<typename B::GridView, TargetSpace>
-{
-
-public:
- typedef B Basis;
-
- typedef typename Basis::LocalFiniteElement LocalFiniteElement;
- typedef typename Basis::GridView GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::Grid::ctype ctype;
-
- typedef LocalGeodesicFEFunction<GridView::dimension, ctype, LocalFiniteElement, TargetSpace> LocalGFEFunction;
- typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
-
- //! Dimension of the grid.
- constexpr static int gridDim = GridView::dimension;
-
- //! Dimension of the embedded tangent space
- constexpr static int embeddedDim = EmbeddedTangentVector::dimension;
-
-
- //! Create global function by a global basis and the corresponding coefficient vector
- GlobalGeodesicFEFunction(const Basis& basis, const std::vector<TargetSpace>& coefficients) :
- VirtualGridViewFunction<typename B::GridView, TargetSpace>(basis.getGridView()),
- basis_(basis),
- coefficients_(coefficients)
- {}
-
-
- /** \brief Evaluate the function at local coordinates. */
- void evaluateLocal(const Element& element, const Dune::FieldVector<ctype,gridDim>& local, TargetSpace& out) const
- {
- int numOfBaseFct = basis_.getLocalFiniteElement(element).localBasis().size();
-
- // Extract local coefficients
- std::vector<TargetSpace> localCoeff(numOfBaseFct);
-
- for (int i=0; i<numOfBaseFct; i++)
- localCoeff[i] = coefficients_[basis_.index(element,i)];
-
- // create local gfe function
- LocalGFEFunction localGFE(basis_.getLocalFiniteElement(element),localCoeff);
- out = localGFE.evaluate(local);
- }
-
- /** \brief Evaluate the derivative of the function at local coordinates. */
- void evaluateDerivativeLocal(const Element& element, const Dune::FieldVector<ctype,gridDim>& local,
- Dune::FieldMatrix<ctype, embeddedDim, gridDim>& out) const
- {
- int numOfBaseFct = basis_.getLocalFiniteElement(element).localBasis().size();
-
- // Extract local coefficients
- std::vector<TargetSpace> localCoeff(numOfBaseFct);
-
- for (int i=0; i<numOfBaseFct; i++)
- localCoeff[i] = coefficients_[basis_.index(element,i)];
-
- // create local gfe function
- LocalGFEFunction localGFE(basis_.getLocalFiniteElement(element),localCoeff);
-
- // use it to evaluate the derivative
- Dune::FieldMatrix<ctype, embeddedDim, gridDim> refJac = localGFE.evaluateDerivative(local);
-
- out =0.0;
- //transform the gradient
- const auto jacInvTrans = element.geometry().jacobianInverseTransposed(local);
- for (size_t k=0; k< refJac.N(); k++)
- jacInvTrans.umv(refJac[k],out[k]);
-
- }
-
- /** \brief Export basis */
- const Basis& basis() const
- {
- return basis_;
- }
-
- /** \brief Export coefficients. */
- const std::vector<TargetSpace>& coefficients() const
- {
- return coefficients_;
- }
-
-private:
- //! The global basis
- const Basis& basis_;
- //! The coefficient vector
- const std::vector<TargetSpace>& coefficients_;
-};
-#endif
diff --git a/dune/gfe/globalgfefunction.hh b/dune/gfe/globalgfefunction.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bc18312af5be91d3e7310bd9e1e2293ae19d9c5f
--- /dev/null
+++ b/dune/gfe/globalgfefunction.hh
@@ -0,0 +1,459 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_GFE_GLOBALGFEFUNCTION_HH
+#define DUNE_GFE_GLOBALGFEFUNCTION_HH
+
+#include <memory>
+#include <optional>
+#include <vector>
+
+#include <dune/common/typetraits.hh>
+
+#include <dune/grid/utility/hierarchicsearch.hh>
+
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+#include <dune/functions/gridfunctions/gridfunction.hh>
+#include <dune/functions/backends/concepts.hh>
+
+namespace Dune::GFE {
+
+namespace Impl {
+
+/** \brief Common base class for GlobalGFEFunction and its derivative
+ *
+ * \tparam B Scalar(!) function-space basis
+ * \tparam V Container of coefficients
+ * \tparam LocalInterpolationRule How to interpolate manifold-valued data
+ */
+template<typename B, typename V, typename LocalInterpolationRule>
+class GlobalGFEFunctionBase
+{
+public:
+ using Basis = B;
+ using Vector = V;
+
+ // In order to make the cache work for proxy-references
+ // we have to use AutonomousValue<T> instead of std::decay_t<T>
+ using Coefficient = Dune::AutonomousValue<decltype(std::declval<Vector>()[std::declval<typename Basis::MultiIndex>()])>;
+
+ using GridView = typename Basis::GridView;
+ using EntitySet = Functions::GridViewEntitySet<GridView, 0>;
+ using Tree = typename Basis::LocalView::Tree;
+
+ using Domain = typename EntitySet::GlobalCoordinate;
+
+ using LocalDomain = typename EntitySet::LocalCoordinate;
+ using Element = typename EntitySet::Element;
+
+protected:
+
+ // This collects all data that is shared by all related
+ // global and local functions. This way we don't need to
+ // keep track of it individually.
+ struct Data
+ {
+ EntitySet entitySet;
+ std::shared_ptr<const Basis> basis;
+ std::shared_ptr<const Vector> coefficients;
+ };
+
+public:
+ class LocalFunctionBase
+ {
+ using LocalView = typename Basis::LocalView;
+ using size_type = typename Tree::size_type;
+
+ public:
+ using Domain = LocalDomain;
+ using Element = typename EntitySet::Element;
+
+ protected:
+ LocalFunctionBase(const std::shared_ptr<const Data>& data)
+ : data_(data)
+ , localView_(data_->basis->localView())
+ {
+ localDoFs_.reserve(localView_.maxSize());
+ }
+
+ /**
+ * \brief Copy-construct the local-function.
+ *
+ * This copy-constructor copies the cached local DOFs only
+ * if the `other` local-function is bound to an element.
+ **/
+ LocalFunctionBase(const LocalFunctionBase& other)
+ : data_(other.data_)
+ , localView_(other.localView_)
+ {
+ localDoFs_.reserve(localView_.maxSize());
+ if (bound())
+ localDoFs_ = other.localDoFs_;
+ }
+
+ /**
+ * \brief Copy-assignment of the local-function.
+ *
+ * Assign all members from `other` to `this`, except the
+ * local DOFs. Those are copied only if the `other`
+ * local-function is bound to an element.
+ **/
+ LocalFunctionBase& operator=(const LocalFunctionBase& other)
+ {
+ data_ = other.data_;
+ localView_ = other.localView_;
+ if (bound())
+ localDoFs_ = other.localDoFs_;
+ return *this;
+ }
+
+ public:
+ /**
+ * \brief Bind LocalFunction to grid element.
+ *
+ * You must call this method before `operator()`
+ * and after changes to the coefficient vector.
+ */
+ void bind(const Element& element)
+ {
+ localView_.bind(element);
+
+ localDoFs_.resize(localView_.size());
+ const auto& dofs = *data_->coefficients;
+ for (size_type i = 0; i < localView_.tree().size(); ++i)
+ {
+ // For a subspace basis the index-within-tree i
+ // is not the same as the localIndex within the
+ // full local view.
+ size_t localIndex = localView_.tree().localIndex(i);
+ localDoFs_[localIndex] = dofs[localView_.index(localIndex)];
+ }
+
+ // create local GFE function
+ // TODO Store this object by value
+ localInterpolationRule_ = std::make_unique<LocalInterpolationRule>(this->localView_.tree().finiteElement(),localDoFs_);
+ }
+
+ //! Unbind the local-function.
+ void unbind()
+ {
+ localView_.unbind();
+ }
+
+ //! Check if LocalFunction is already bound to an element.
+ bool bound() const
+ {
+ return localView_.bound();
+ }
+
+ //! Return the element the local-function is bound to.
+ const Element& localContext() const
+ {
+ return localView_.element();
+ }
+
+ protected:
+
+ std::shared_ptr<const Data> data_;
+ LocalView localView_;
+ std::vector<Coefficient> localDoFs_;
+ std::unique_ptr<LocalInterpolationRule> localInterpolationRule_;
+ };
+
+protected:
+ GlobalGFEFunctionBase(const std::shared_ptr<const Data>& data)
+ : data_(data)
+ {
+ /* Nothing. */
+ }
+
+public:
+
+ //! Return a const reference to the stored basis.
+ const Basis& basis() const
+ {
+ return *data_->basis;
+ }
+
+ //! Return the coefficients of this discrete function by reference.
+ const Vector& dofs() const
+ {
+ return *data_->coefficients;
+ }
+
+ //! Get associated set of entities the local-function can be bound to.
+ const EntitySet& entitySet() const
+ {
+ return data_->entitySet;
+ }
+
+protected:
+ std::shared_ptr<const Data> data_;
+};
+
+} // namespace Impl
+
+
+
+template<typename GGF>
+class GlobalGFEFunctionDerivative;
+
+/**
+ * \brief A global geometric finite element function
+ *
+ * \tparam B Type of global basis
+ * \tparam LIR Local interpolation rule for manifold-valued data
+ * \tparam TargetSpace Range type of this function
+ */
+template<typename B, typename LIR, typename TargetSpace>
+class GlobalGFEFunction
+ : public Impl::GlobalGFEFunctionBase<B, std::vector<TargetSpace>, LIR>
+{
+ using Base = Impl::GlobalGFEFunctionBase<B, std::vector<TargetSpace>, LIR>;
+ using Data = typename Base::Data;
+
+public:
+ using Basis = typename Base::Basis;
+ using Vector = typename Base::Vector;
+ using LocalInterpolationRule = LIR;
+
+ using Domain = typename Base::Domain;
+ using Range = typename TargetSpace::CoordinateType;
+
+ using Traits = Functions::Imp::GridFunctionTraits<Range(Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
+
+ class LocalFunction
+ : public Base::LocalFunctionBase
+ {
+ using LocalBase = typename Base::LocalFunctionBase;
+ using size_type = typename Base::Tree::size_type;
+
+ public:
+
+ using GlobalFunction = GlobalGFEFunction;
+ using Domain = typename LocalBase::Domain;
+ using Range = GlobalFunction::Range;
+ using Element = typename LocalBase::Element;
+
+ //! Create a local-function from the associated grid-function
+ LocalFunction(const GlobalGFEFunction& globalFunction)
+ : LocalBase(globalFunction.data_)
+ {
+ /* Nothing. */
+ }
+
+ /**
+ * \brief Evaluate this local-function in coordinates `x` in the bound element.
+ *
+ * The result of this method is undefined if you did
+ * not call bind() beforehand or changed the coefficient
+ * vector after the last call to bind(). In the latter case
+ * you have to call bind() again in order to make operator()
+ * usable.
+ */
+ Range operator()(const Domain& x) const
+ {
+ return this->localInterpolationRule_->evaluate(x).globalCoordinates();
+ }
+
+ //! Local function of the derivative
+ friend typename GlobalGFEFunctionDerivative<GlobalGFEFunction>::LocalFunction derivative(const LocalFunction& lf)
+ {
+ auto dlf = localFunction(GlobalGFEFunctionDerivative<GlobalGFEFunction>(lf.data_));
+ if (lf.bound())
+ dlf.bind(lf.localContext());
+ return dlf;
+ }
+ };
+
+ //! Create a grid-function, by wrapping the arguments in `std::shared_ptr`.
+ template<class B_T, class V_T>
+ GlobalGFEFunction(B_T && basis, V_T && coefficients)
+ : Base(std::make_shared<Data>(Data{{basis.gridView()}, wrap_or_move(std::forward<B_T>(basis)), wrap_or_move(std::forward<V_T>(coefficients))}))
+ {}
+
+ //! Create a grid-function, by moving the arguments in `std::shared_ptr`.
+ GlobalGFEFunction(std::shared_ptr<const Basis> basis, std::shared_ptr<const Vector> coefficients)
+ : Base(std::make_shared<Data>(Data{{basis->gridView()}, basis, coefficients}))
+ {}
+
+ /** \brief Evaluate at a point given in world coordinates
+ *
+ * \warning This has to find the element that the evaluation point is in.
+ * It is therefore very slow.
+ */
+ Range operator() (const Domain& x) const
+ {
+ HierarchicSearch search(this->data_->basis->gridView().grid(), this->data_->basis->gridView().indexSet());
+
+ const auto e = search.findEntity(x);
+ auto localThis = localFunction(*this);
+ localThis.bind(e);
+ return localThis(e.geometry().local(x));
+ }
+
+ //! Derivative of the `GlobalGFEFunction`
+ friend GlobalGFEFunctionDerivative<GlobalGFEFunction> derivative(const GlobalGFEFunction& f)
+ {
+ return GlobalGFEFunctionDerivative<GlobalGFEFunction>(f.data_);
+ }
+
+ /**
+ * \brief Construct local function from a GlobalGFEFunction.
+ *
+ * The obtained local function satisfies the concept
+ * `Dune::Functions::Concept::LocalFunction`. It must be bound
+ * to an entity from the entity set of the GlobalGFEFunction
+ * before it can be used.
+ */
+ friend LocalFunction localFunction(const GlobalGFEFunction& t)
+ {
+ return LocalFunction(t);
+ }
+};
+
+
+/**
+ * \brief Derivative of a `GlobalGFEFunction`
+ *
+ * Function returning the derivative of the given `GlobalGFEFunction`
+ * with respect to global coordinates.
+ *
+ * \tparam GGF instance of the `GlobalGFEFunction` this is a derivative of
+ */
+template<typename GGF>
+class GlobalGFEFunctionDerivative
+ : public Impl::GlobalGFEFunctionBase<typename GGF::Basis, typename GGF::Vector, typename GGF::LocalInterpolationRule>
+{
+ using Base = Impl::GlobalGFEFunctionBase<typename GGF::Basis, typename GGF::Vector, typename GGF::LocalInterpolationRule>;
+ using Data = typename Base::Data;
+
+public:
+ using GlobalGFEFunction = GGF;
+
+ using Basis = typename Base::Basis;
+ using Vector = typename Base::Vector;
+
+ using Domain = typename Base::Domain;
+ using Range = typename Functions::SignatureTraits<typename GlobalGFEFunction::Traits::DerivativeInterface>::Range;
+
+ using Traits = Functions::Imp::GridFunctionTraits<Range(Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
+
+ /**
+ * \brief local function evaluating the derivative in reference coordinates
+ *
+ * Note that the function returns the derivative with respect to global
+ * coordinates even when the point is given in reference coordinates on
+ * an element.
+ */
+ class LocalFunction
+ : public Base::LocalFunctionBase
+ {
+ using LocalBase = typename Base::LocalFunctionBase;
+ using size_type = typename Base::Tree::size_type;
+
+ public:
+ using GlobalFunction = GlobalGFEFunctionDerivative;
+ using Domain = typename LocalBase::Domain;
+ using Range = GlobalFunction::Range;
+ using Element = typename LocalBase::Element;
+
+ //! Create a local function from the associated grid function
+ LocalFunction(const GlobalFunction& globalFunction)
+ : LocalBase(globalFunction.data_)
+ {
+ /* Nothing. */
+ }
+
+ /**
+ * \brief Bind LocalFunction to grid element.
+ *
+ * You must call this method before `operator()`
+ * and after changes to the coefficient vector.
+ */
+ void bind(const Element& element)
+ {
+ LocalBase::bind(element);
+ geometry_.emplace(element.geometry());
+ }
+
+ //! Unbind the local-function.
+ void unbind()
+ {
+ geometry_.reset();
+ LocalBase::unbind();
+ }
+
+ /**
+ * \brief Evaluate this local-function in coordinates `x` in the bound element.
+ *
+ * The result of this method is undefined if you did
+ * not call bind() beforehand or changed the coefficient
+ * vector after the last call to bind(). In the latter case
+ * you have to call bind() again in order to make operator()
+ * usable.
+ *
+ * Note that the function returns the derivative with respect to global
+ * coordinates even though the evaluation point is given in reference coordinates
+ * on the current element.
+ */
+ Range operator()(const Domain& x) const
+ {
+ // Jacobian with respect to local coordinates
+ auto refJac = this->localInterpolationRule_->evaluateDerivative(x);
+
+ // Transform to world coordinates
+ return refJac * geometry_->jacobianInverse(x);
+ }
+
+ //! Not implemented
+ friend typename Traits::LocalFunctionTraits::DerivativeInterface derivative(const LocalFunction&)
+ {
+ DUNE_THROW(NotImplemented, "derivative of derivative is not implemented");
+ }
+
+ private:
+ std::optional<typename Element::Geometry> geometry_;
+ };
+
+ /**
+ * \brief create object from `GlobalGFEFunction` data
+ *
+ * Please call `derivative(globalGFEFunction)` to create an instance
+ * of this class.
+ */
+ GlobalGFEFunctionDerivative(const std::shared_ptr<const Data>& data)
+ : Base(data)
+ {
+ /* Nothing. */
+ }
+
+ /** \brief Evaluate the discrete grid-function derivative in global coordinates
+ *
+ * \warning This has to find the element that the evaluation point is in.
+ * It is therefore very slow.
+ */
+ Range operator()(const Domain& x) const
+ {
+ HierarchicSearch search(this->data_->basis->gridView().grid(), this->data_->basis->gridView().indexSet());
+
+ const auto e = search.findEntity(x);
+ auto localThis = localFunction(*this);
+ localThis.bind(e);
+ return localThis(e.geometry().local(x));
+ }
+
+ friend typename Traits::DerivativeInterface derivative(const GlobalGFEFunctionDerivative& f)
+ {
+ DUNE_THROW(NotImplemented, "derivative of derivative is not implemented");
+ }
+
+ //! Construct local function from a `GlobalGFEFunctionDerivative`
+ friend LocalFunction localFunction(const GlobalGFEFunctionDerivative& f)
+ {
+ return LocalFunction(f);
+ }
+};
+
+} // namespace Dune::GFE
+
+#endif // DUNE_GFE_GLOBALGFEFUNCTION_HH
diff --git a/dune/gfe/l2distancesquaredenergy.hh b/dune/gfe/l2distancesquaredenergy.hh
index 2546e70d5c5867f1f33bde7a8dc917e6dcec1c3b..05e460c9f6dad9546912aaa86dff055130352ffa 100644
--- a/dune/gfe/l2distancesquaredenergy.hh
+++ b/dune/gfe/l2distancesquaredenergy.hh
@@ -5,6 +5,7 @@
#include <dune/geometry/quadraturerules.hh>
+#include <dune/gfe/globalgfefunction.hh>
#include <dune/gfe/localenergy.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
@@ -20,10 +21,12 @@ class L2DistanceSquaredEnergy
// some other sizes
constexpr static int gridDim = GridView::dimension;
+ using LocalInterpolationRule = LocalGeodesicFEFunction<gridDim, DT, typename Basis::LocalView::Tree::FiniteElement, typename TargetSpace::template rebind<double>::other>;
+
public:
// This is the function that we are computing the L2-distance to
- std::shared_ptr<VirtualGridViewFunction<GridView,typename TargetSpace::template rebind<double>::other > > origin_;
+ std::shared_ptr<Dune::GFE::GlobalGFEFunction<Basis, LocalInterpolationRule, typename TargetSpace::template rebind<double>::other > > origin_;
/** \brief Assemble the energy for a single element */
RT energy (const typename Basis::LocalView& localView,
@@ -36,11 +39,13 @@ public:
typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
+ const auto element = localView.element();
+ auto localOrigin = localFunction(*origin_);
+ localOrigin.bind(element);
+
// Just guessing an appropriate quadrature order
auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim;
- const auto element = localView.element();
-
const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++)
@@ -54,21 +59,14 @@ public:
// The function value
auto value = localGeodesicFEFunction.evaluate(quadPos);
- typename TargetSpace::template rebind<double>::other originValue;
- origin_->evaluateLocal(element,quadPos, originValue);
+ auto originValue = localOrigin(quadPos);
// The derivative of the 'origin' function
// First: as function defined on the reference element
- typename VirtualGridViewFunction<GridView,typename TargetSpace::template rebind<double>::other>::DerivativeType originReferenceDerivative;
- origin_->evaluateDerivativeLocal(element,quadPos,originReferenceDerivative);
+ auto originReferenceDerivative = derivative(localOrigin)(quadPos);
// The derivative of the function defined on the actual element
- typename VirtualGridViewFunction<GridView,typename TargetSpace::template rebind<double>::other >::DerivativeType originDerivative(0);
-
- auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
- for (size_t comp=0; comp<originReferenceDerivative.N(); comp++)
- jacobianInverseTransposed.umv(originReferenceDerivative[comp], originDerivative[comp]);
+ auto originDerivative = originReferenceDerivative * element.geometry().jacobianInverse(quadPos);
double weightFactor = originDerivative.frobenius_norm();
// un-comment the following line to switch off the weight factor
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index ea1be930de232f2fb8b5f7dd6ccaa777f195886f..7ca443c6aeee032a7acd084179e888b7a2cf8753 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -25,11 +25,10 @@
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
-#include <dune/fufem/discretizationerror.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
@@ -39,7 +38,7 @@
#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/geodesicfeassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/globalgeodesicfefunction.hh>
+#include <dune/gfe/globalgfefunction.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
#include <dune/gfe/harmonicenergy.hh>
#include <dune/gfe/l2distancesquaredenergy.hh>
@@ -137,22 +136,19 @@ int main (int argc, char *argv[]) try
}
grid->globalRefine(numLevels-1);
+ auto gridView = grid->leafGridView();
//////////////////////////////////////////////////////////////////////////////////
// Construct the scalar function space basis corresponding to the GFE space
//////////////////////////////////////////////////////////////////////////////////
typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, order> FEBasis;
- //typedef Dune::Functions::Periodic1DPQ1NodalBasis<typename GridType::LeafGridView> FEBasis;
FEBasis feBasis(grid->leafGridView());
///////////////////////////////////////////
// Read Dirichlet values
///////////////////////////////////////////
- typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
- FufemFEBasis fufemFeBasis(feBasis);
-
BitSetVector<1> dirichletVertices(grid->leafGridView().size(dim), false);
const auto& indexSet = grid->leafGridView().indexSet();
@@ -171,7 +167,7 @@ int main (int argc, char *argv[]) try
BoundaryPatch<GridType::LeafGridView> dirichletBoundary(grid->leafGridView(), dirichletVertices);
BitSetVector<blocksize> dirichletNodes(feBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,fufemFeBasis,dirichletNodes);
+ constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
////////////////////////////
// Initial iterate
@@ -182,7 +178,16 @@ int main (int argc, char *argv[]) try
auto pythonInitialIterate = Python::make_function<TargetSpace::CoordinateType>(module.get("f"));
std::vector<TargetSpace::CoordinateType> v;
- Functions::interpolate(feBasis, v, pythonInitialIterate);
+ using namespace Functions::BasisFactory;
+
+ auto powerBasis = makeBasis(
+ gridView,
+ power<TargetSpace::CoordinateType::dimension>(
+ lagrange<order>(),
+ blockedInterleaved()
+ ));
+
+ Functions::interpolate(powerBasis, v, pythonInitialIterate);
SolutionType x(feBasis.size());
@@ -259,7 +264,7 @@ int main (int argc, char *argv[]) try
for (int i=0; i<numTimeSteps; i++)
{
- auto previousTimeStepFct = std::make_shared<GlobalGeodesicFEFunction<FufemFEBasis,TargetSpace> >(fufemFeBasis,previousTimeStep);
+ auto previousTimeStepFct = std::make_shared<GFE::GlobalGFEFunction<FEBasis,GeodesicInterpolationRule::rebind<TargetSpace>::other,TargetSpace> >(feBasis,previousTimeStep);
l2DistanceSquaredEnergy->origin_ = previousTimeStepFct;
solver.setInitialIterate(x);