Merge branch 'feature/vintage-gridfunction' into 'master'

Add VintageBasisGridFunction

See merge request !85

dune/fufem/functions/CMakeLists.txt

@@ -15,6 +15,7 @@ install(FILES
     portablepixelmap.hh
     sumfunction.hh
     sumgridfunction.hh
+    vintagebasisgridfunction.hh
     virtualdifferentiablefunction.hh
     virtualfunctiontoboundarysegmentadapter.hh
     virtualgridfunction.hh

dune/fufem/functions/vintagebasisgridfunction.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FUFEM_FUNCTIONS_VINTAGEBASISGRIDFUNCTION_HH
+#define DUNE_FUFEM_FUNCTIONS_VINTAGEBASISGRIDFUNCTION_HH
+
+#include <type_traits>
+#include <optional>
+#include <memory>
+
+#include <dune/common/shared_ptr.hh>
+#include <dune/common/exceptions.hh>
+
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+#include <dune/functions/gridfunctions/gridfunction.hh>
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+
+namespace Dune::Fufem {
+
+namespace Impl {
+
+  /**
+   * \brief A dune-functions grid function implemented using a dune-fufem basis
+   *
+   * This allows to use a dune-fufem basis as dune-functions
+   * grid function. Moreover you can customize the local function
+   * evaluation.
+   */
+  template<typename B, typename V, typename LEV, typename R = typename V::value_type>
+  class VintageBasisGridFunction
+  {
+  public:
+    using Basis = B;
+    using Vector = V;
+    using LocalEvaluator = LEV;
+
+    using Coefficient = std::decay_t<decltype(std::declval<Vector>()[0])>;
+
+    using GridView = typename Basis::GridView;
+    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
+
+    using Domain = typename EntitySet::GlobalCoordinate;
+    using Range = R;
+
+    using LocalDomain = typename EntitySet::LocalCoordinate;
+    using Element = typename EntitySet::Element;
+
+    using Traits = Dune::Functions::Imp::GridFunctionTraits<Range(Domain), EntitySet, Dune::Functions::DefaultDerivativeTraits, 16>;
+
+    class LocalFunction
+    {
+    public:
+      using GlobalFunction = VintageBasisGridFunction;
+      using Domain = GlobalFunction::LocalDomain;
+      using Range = GlobalFunction::Range;
+      using Element = GlobalFunction::Element;
+
+      LocalFunction(const VintageBasisGridFunction& globalFunction)
+        : globalFunction_(&globalFunction)
+      {}
+
+      LocalFunction(const LocalFunction& other) = default;
+      LocalFunction(LocalFunction&& other) = default;
+      LocalFunction& operator=(const LocalFunction& other) = default;
+      LocalFunction& operator=(LocalFunction&& other) = default;
+
+      /**
+       * \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) {
+        element_ = element;
+      }
+
+      void unbind() {
+        element_.reset();
+      }
+
+      bool bound() const {
+        return static_cast<bool>(element_);
+      }
+
+      Range operator()(const Domain& x) const {
+        auto y = Range(0);
+        globalFunction_->localEvaluator()(
+            globalFunction_->basis(),
+            globalFunction_->dofs(),
+            *element_,
+            x,
+            y);
+        return y;
+      }
+
+      const Element& localContext() const {
+        return *element_;
+      }
+
+      friend typename Traits::LocalFunctionTraits::DerivativeInterface derivative(const LocalFunction& t) {
+        DUNE_THROW(NotImplemented,"not implemented");
+      }
+
+    private:
+
+      const VintageBasisGridFunction* globalFunction_;
+      std::optional<Element> element_;
+    };
+
+    template<class B_T, class V_T, class LEV_T>
+    VintageBasisGridFunction(B_T&& basis, V_T&& v, LEV_T&& localEvaluator)
+      : entitySet_(basis.getGridView())
+      , basis_(wrap_or_move(std::forward<B_T>(basis)))
+      , coefficients_(wrap_or_move(std::forward<V_T>(v)))
+      , localEvaluator_(std::forward<LEV_T>(localEvaluator))
+    {}
+
+    template<class B_T, class V_T>
+    VintageBasisGridFunction(B_T&& basis, V_T&& v)
+      : entitySet_(basis.getGridView())
+      , basis_(wrap_or_move(std::forward<B_T>(basis)))
+      , coefficients_(wrap_or_move(std::forward<V_T>(v)))
+      , localEvaluator_()
+    {}
+
+    const Basis& basis() const {
+      return *basis_;
+    }
+
+    const V& dofs() const {
+      return *coefficients_;
+    }
+
+    const LEV& localEvaluator() const {
+      return localEvaluator_;
+    }
+
+    Range operator() (const Domain& x) const {
+      DUNE_THROW(NotImplemented,"not implemented");
+    }
+
+    friend typename Traits::DerivativeInterface derivative(const VintageBasisGridFunction& t) {
+      DUNE_THROW(NotImplemented,"not implemented");
+    }
+
+    friend LocalFunction localFunction(const VintageBasisGridFunction& t) {
+      return LocalFunction(t);
+    }
+
+    const EntitySet& entitySet() const {
+      return entitySet_;
+    }
+
+  private:
+
+    EntitySet entitySet_;
+    std::shared_ptr<const Basis> basis_;
+    std::shared_ptr<const V> coefficients_;
+    LocalEvaluator localEvaluator_;
+  };
+
+  /**
+   * \brief Construction of local functions from a temporary VintageBasisGridFunction (forbidden)
+   *
+   * Since a VintageBasisGridFunction::LocalFunction stores a reference
+   * to the global VintageBasisGridFunction its life time is bound to
+   * the latter. Hence construction from a temporary VintageBasisGridFunction
+   * would lead to a dangling reference and is thus forbidden/deleted.
+   */
+  template<typename... TT>
+  void localFunction(VintageBasisGridFunction<TT...>&& t) = delete;
+
+
+} // end namespace Impl
+
+
+/**
+ * \brief Create a VintageBasisGridFunction
+ *
+ * This allows to use a dune-fufem basis as dune-functions
+ * grid function. Moreover you can customize the local function
+ * evaluation.
+ *
+ * \tparam R Range type of function
+ * \tparam B Type of dune-fufem basis
+ * \tparam V Type of coefficient vector
+ * \tparam LEV Type of local evaluator
+ *
+ * \param basis Dune-fufem basis to interpret the coefficient vector
+ * \param v Coefficient vector of FE function wrt the given basis
+ * \param localEvaluator A function doing the local evaluation.
+ *
+ * The actual local function implementation is forwarded to the
+ * callback localEvaluator whose signature
+ * localEvaluator(basis, coefficients, element, localX, y)
+ * mimics the LocalEvaluator::evaluateLocal() function.
+ * By implementing a custom callback you can e.g.
+ * represent partial derivatives.
+ */
+template<typename R, typename B, typename V, typename LEV>
+auto makeVintageBasisGridFunction(B&& basis, V&& v, LEV&& localEvaluator)
+{
+  using Basis = std::decay_t<B>;
+  using Vector = std::decay_t<V>;
+  using LocalEvaluator = std::decay_t<LEV>;
+  return Impl::VintageBasisGridFunction<Basis, Vector, LocalEvaluator, R>(std::forward<B>(basis), std::forward<V>(v), std::forward<LEV>(localEvaluator));
+}
+
+/**
+ * \brief Create a VintageBasisGridFunction
+ *
+ * This allows to use a dune-fufem basis as dune-functions
+ * grid function.
+ *
+ * \tparam R Range type of function
+ * \tparam B Type of dune-fufem basis
+ * \tparam V Type of coefficient vector
+ *
+ * \param basis Dune-fufem basis to interpret the coefficient vector
+ * \param v Coefficient vector of FE function wrt the given basis
+ */
+template<typename R, typename B, typename V>
+auto makeVintageBasisGridFunction(B&& basis, V&& v)
+{
+  using Basis = std::decay_t<B>;
+  using Vector = std::decay_t<V>;
+  auto localEvaluator = [](const auto& b, const auto& v, const auto&e, const auto& x, auto& y) {
+    LocalEvaluator<Basis, Vector, R>::evaluateLocal(b, v, e, x, y);
+  };
+  using LocalEvaluator = decltype(localEvaluator);
+  return Impl::VintageBasisGridFunction<Basis, Vector, LocalEvaluator, R>(std::forward<B>(basis), std::forward<V>(v), std::move(localEvaluator));
+}
+
+} // end namespace Dune::Fufem
+
+
+
+
+
+
+
+#endif // DUNE_FUFEM_FUNCTIONS_VINTAGEBASISGRIDFUNCTION_HH
+
+

dune/fufem/test/CMakeLists.txt

@@ -28,6 +28,7 @@ set(GRID_BASED_TESTS
   sumfunctiontest
   tensortest
   transferoperatorassemblertest
+  vintagebasisgridfunctiontest
   )
 
 if (ADOLC_FOUND)

dune/fufem/test/vintagebasisgridfunctiontest.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/test/gridfunctiontest.hh>
+
+#include <dune/fufem/functions/vintagebasisgridfunction.hh>
+#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
+
+using namespace Dune;
+using namespace Dune::Functions;
+using namespace Dune::Functions::Test;
+
+
+
+
+int main (int argc, char* argv[]) try
+{
+  using namespace Dune;
+  using namespace Dune::Functions;
+  using namespace Dune::Functions::BasisFactory;
+
+  Dune::MPIHelper::instance(argc, argv);
+  bool passed = true;
+
+
+  // Generate grid for testing
+  const int dim = 2;
+  
+  using Grid = Dune::YaspGrid<dim>;
+
+  auto l = Dune::FieldVector<double,dim>(1);
+  auto elements = std::array<int,dim>{{10, 10}};
+  auto grid = Grid(l,elements);
+  grid.globalRefine(2);
+
+  using GridView = Grid::LeafGridView;
+  using Domain = GridView::template Codim<0>::Geometry::GlobalCoordinate;
+
+  auto gridView = grid.leafGridView();
+
+
+
+  {
+    auto duneFunctionsBasis = makeBasis(gridView, lagrange<3>());
+    auto duneFufemBasis = DuneFunctionsBasis<decltype(duneFunctionsBasis)>(duneFunctionsBasis);
+
+    using Range = double;
+    using Vector = std::vector<Range>;
+    auto f = [](const Domain& x){
+      return x[0];
+    };
+    double exactIntegral = 0.5;
+
+    // Construct some coefficient vector by interpolation
+    auto x = Vector{};
+    Dune::Functions::interpolate(duneFunctionsBasis, x, f);
+
+    // Create and use VintageBasisGridFunction
+    auto gf = Dune::Fufem::makeVintageBasisGridFunction<Range>(duneFufemBasis, x);
+    auto y = Vector{};
+    Dune::Functions::interpolate(duneFunctionsBasis, y, gf);
+
+    passed = passed and (x==y);
+
+    passed = passed and Dune::Functions::Test::checkGridViewFunction(gridView, gf, exactIntegral);
+  }
+
+  return passed ? 0: 1;
+} catch ( Dune::Exception &e )
+{
+  std::cerr << "Dune reported error: " << e << std::endl;
+  return 1;
+}
+catch(...)
+{
+  std::cerr << "Unknown exception thrown!" << std::endl;
+  return 1;
+}