Merge branch 'feature/vtkfunction' into 'master'

Feature/vtkfunction

See merge request spraetor/dune-vtk!3

dune/vtk/defaultvtkfunction.hh

+#pragma once
+
+#include "vtklocalfunctioninterface.hh"
+
+namespace Dune { namespace experimental
+{
+  /// Type erasure for dune-functions LocalFunction interface
+  template <class GridView, class LocalFunction>
+  class LocalFunctionWrapper final
+      : public VtkLocalFunctionInterface<GridView>
+  {
+    using Self = LocalFunctionWrapper;
+    using Interface = VtkLocalFunctionInterface<GridView>;
+    using Entity = typename Interface::Entity;
+    using LocalCoordinate = typename Interface::LocalCoordinate;
+
+    template <class F, class D>
+    using Range = std::decay_t<decltype(std::declval<F>()(std::declval<D>()))>;
+
+    template <class F, class D>
+    using VectorValued = decltype(std::declval<Range<F,D>>()[0u]);
+
+  public:
+    template <class LocalFct, disableCopyMove<Self, LocalFct> = 0>
+    LocalFunctionWrapper (LocalFct&& localFct)
+      : localFct_(std::forward<LocalFct>(localFct))
+    {}
+
+    virtual void bind (Entity const& entity) override
+    {
+      localFct_.bind(entity);
+    }
+
+    virtual void unbind () override
+    {
+      localFct_.unbind();
+    }
+
+    virtual double evaluate (int comp, LocalCoordinate const& xi) const override
+    {
+      return evaluateImpl(comp, xi, Std::is_detected<VectorValued,LocalFunction,LocalCoordinate>{});
+    }
+
+  private:
+    // Evaluate a component of a vector valued data
+    double evaluateImpl (int comp, LocalCoordinate const& xi, std::true_type) const
+    {
+      auto y = localFct_(xi);
+      return comp < y.size() ? y[comp] : 0.0;
+    }
+
+    // Return the scalar values
+    double evaluateImpl (int comp, LocalCoordinate const& xi, std::false_type) const
+    {
+      assert(comp == 0);
+      return localFct_(xi);
+    }
+
+  private:
+    LocalFunction localFct_;
+  };
+
+}} // end namespace Dune::experimental

dune/vtk/legacyvtkfunction.hh

+#pragma once
+
+#include <memory>
+
+#include <dune/grid/io/file/vtk/function.hh>
+
+#include "vtklocalfunctioninterface.hh"
+
+namespace Dune { namespace experimental
+{
+  /// Type erasure for Legacy VTKFunction
+  template <class GridView>
+  class VTKLocalFunctionWrapper final
+      : public VtkLocalFunctionInterface<GridView>
+  {
+    using Interface = VtkLocalFunctionInterface<GridView>;
+    using Entity = typename Interface::Entity;
+    using LocalCoordinate = typename Interface::LocalCoordinate;
+
+  public:
+    VTKLocalFunctionWrapper (std::shared_ptr<VTKFunction<GridView> const> const& fct)
+      : fct_(fct)
+    {}
+
+    virtual void bind (Entity const& entity) override
+    {
+      entity_ = &entity;
+    }
+
+    virtual void unbind () override
+    {
+      entity_ = nullptr;
+    }
+
+    virtual double evaluate (int comp, LocalCoordinate const& xi) const override
+    {
+      return fct_->evaluate(comp, *entity_, xi);
+    }
+
+  private:
+    std::shared_ptr<VTKFunction<GridView> const> fct_;
+    Entity const* entity_;
+  };
+
+}} // end namespace Dune::experimental

dune/vtk/vtkfunction.hh

@@ -4,201 +4,60 @@
 
 #include <dune/common/std/type_traits.hh>
 #include <dune/functions/common/signature.hh>
-#include <dune/functions/common/typeerasure.hh>
+
+#include "vtklocalfunction.hh"
 
 namespace Dune { namespace experimental
 {
-  /// An abstract base class for LocalFunctions
-  template <class GridView>
-  class VTKLocalFunctionInterface
-  {
-  public:
-    using Entity = typename GridView::template Codim<0>::Entity;
-    using LocalCoordinate = typename Entity::Geometry::LocalCoordinate;
-
-    /// Bind the function to the grid entity
-    virtual void bind (Entity const& entity) = 0;
-
-    /// Unbind from the currently bound entity
-    virtual void unbind () = 0;
-
-    /// Evaluate single component comp in the entity at local coordinates xi
-    virtual double evaluate (int comp, LocalCoordinate const& xi) const = 0;
-
-    /// Virtual destructor
-    virtual ~VTKLocalFunctionInterface () = default;
-  };
-
-
-  template <class GridView>
-  struct VTKLocalFunctionImpl
-  {
-    template <class Wrapper>
-    class Model : public Wrapper
-    {
-    public:
-      using Wrapper::Wrapper;
-      using Function = typename Wrapper::Wrapped;
-      using Interface = VTKLocalFunctionInterface<GridView>;
-
-      using Entity = typename Interface::Entity;
-      using LocalCoordinate = typename Interface::LocalCoordinate;
-
-      template <class F, class D>
-      using Range = std::decay_t<decltype(std::declval<F>()(std::declval<D>()))>;
-
-      template <class F, class D>
-      using VectorValued = decltype(std::declval<Range<F,D>>()[0u]);
-
-      virtual void bind (Entity const& entity) override
-      {
-        this->get().bind(entity);
-      }
-
-      virtual void unbind () override
-      {
-        this->get().unbind();
-      }
-
-      virtual double evaluate (int comp, LocalCoordinate const& xi) const override
-      {
-        return evaluateImpl(comp, xi, Std::is_detected<VectorValued,Function,LocalCoordinate>{});
-      }
-
-    private:
-      // Evaluate a component of a vector valued data
-      double evaluateImpl (int comp, LocalCoordinate const& xi, std::true_type) const
-      {
-        auto y = this->get()(xi);
-        return comp < y.size() ? y[comp] : 0.0;
-      }
-
-      // Return the scalar values
-      double evaluateImpl (int comp, LocalCoordinate const& xi, std::false_type) const
-      {
-        assert(comp == 0);
-        return this->get()(xi);
-      }
-    };
-  };
-
-
-  template <class GridView>
-  class VTKLocalFunction
-      : public Functions::TypeErasureBase<VTKLocalFunctionInterface<GridView>,
-                                          VTKLocalFunctionImpl<GridView>::template Model>
-  {
-    using Super = Functions::TypeErasureBase<VTKLocalFunctionInterface<GridView>,
-                                             VTKLocalFunctionImpl<GridView>::template Model>;
-
-    using Entity = typename GridView::template Codim<0>::Entity;
-    using LocalCoordinate = typename Entity::Geometry::LocalCoordinate;
-
-  public:
-    template <class F, disableCopyMove<VTKLocalFunction, F> = 0>
-    VTKLocalFunction (F&& f)
-      : Super(std::forward<F>(f))
-    {}
-
-    VTKLocalFunction () = default;
-
-    /// Bind the function to the grid entity
-    void bind (Entity const& entity)
-    {
-      this->asInterface().bind(entity);
-    }
-
-    /// Unbind from the currently bound entity
-    void unbind ()
-    {
-      this->asInterface().unbind();
-    }
-
-    /// Evaluate the `comp` component of the Range value at local coordinate `xi`
-    double evaluate (int comp, LocalCoordinate const& xi) const
-    {
-      return this->asInterface().evaluate(comp, xi);
-    }
-  };
-
-  // ---------------------------------------------------------------------------
-
-  /// An abstract base class for GlobalFunctions
-  template <class GridView>
-  class VTKFunctionInterface
-  {
-  public:
-    /// Create a local function
-    virtual VTKLocalFunction<GridView> makeLocalFunction () const = 0;
-
-    /// Virtual destructor
-    virtual ~VTKFunctionInterface () = default;
-  };
-
-
-  template <class GridView>
-  struct VTKFunctionImpl
-  {
-    template <class Wrapper>
-    class Model : public Wrapper
-    {
-    public:
-      using Wrapper::Wrapper;
-      virtual VTKLocalFunction<GridView> makeLocalFunction () const override
-      {
-        return VTKLocalFunction<GridView>{localFunction(this->get())};
-      }
-    };
-  };
-
-
   template <class GridView>
-  class VTKFunction
-      : public Functions::TypeErasureBase<VTKFunctionInterface<GridView>,
-                                          VTKFunctionImpl<GridView>::template Model>
+  class VtkFunction
   {
-    using Super = Functions::TypeErasureBase<VTKFunctionInterface<GridView>,
-                                             VTKFunctionImpl<GridView>::template Model>;
-
     template <class F>
     using HasLocalFunction = decltype(localFunction(std::declval<F>()));
 
     template <class F>
-    using Signature = typename std::decay_t<F>::Signature;
+    using Domain = typename std::decay_t<F>::EntitySet::GlobalCoordinate;
+
+    template <class F>
+    using Range = std::decay_t<decltype(std::declval<F>()(std::declval<Domain<F>>()))>;
 
   public:
+    /// Constructor VtkFunction from legacy VTKFunction
+    VtkFunction (std::shared_ptr<VTKFunction<GridView> const> const& fct)
+      : localFct_(fct)
+      , name_(fct->name())
+      , ncomps_(fct->ncomps())
+      , type_(Vtk::FLOAT64)
+    {}
+
+    /// Construct VtkFunction from dune-functions GridFunction with Signature
     template <class F,
-      class Range = typename Functions::SignatureTraits<Signature<F>>::Range>
-    VTKFunction (F&& f, std::string name, int ncomps = 1,
-                 Vtk::DataTypes type = Vtk::Map::type<Range>)
-      : Super(std::forward<F>(f))
+      std::enable_if_t<Std::is_detected<HasLocalFunction,F>::value, int> = 0,
+      std::enable_if_t<Std::is_detected<Range,F>::value,int> = 0>
+    VtkFunction (F&& fct, std::string name, int ncomps = 1, Vtk::DataTypes type = Vtk::Map::type<Range<F>>)
+      : localFct_(localFunction(std::forward<F>(fct)))
       , name_(std::move(name))
-      , ncomps_(ncomps > 3 ? 9 : ncomps > 1 ? 3 : 1) // tensor, vector, or scalar
+      , ncomps_(ncomps)
       , type_(type)
-    {
-      static_assert(Std::is_detected<HasLocalFunction,F>::value,
-        "Requires A GridFunction to be passed to the VTKFunction.");
-    }
+    {}
 
+    /// Construct VtkFunction from dune-functions GridFunction without Signature
     template <class F,
-      std::enable_if_t<not Std::is_detected<Signature,F>::value,int> = 0>
-    VTKFunction (F&& f, std::string name, int ncomps = 1,
-                 Vtk::DataTypes type = Vtk::FLOAT32)
-      : Super(std::forward<F>(f))
+      std::enable_if_t<Std::is_detected<HasLocalFunction,F>::value, int> = 0,
+      std::enable_if_t<not Std::is_detected<Range,F>::value,int> = 0>
+    VtkFunction (F const& fct, std::string name, int ncomps = 1, Vtk::DataTypes type = Vtk::FLOAT32)
+      : localFct_(localFunction(std::forward<F>(fct)))
       , name_(std::move(name))
-      , ncomps_(ncomps > 3 ? 9 : ncomps > 1 ? 3 : 1) // tensor, vector, or scalar
+      , ncomps_(ncomps)
       , type_(type)
-    {
-      static_assert(Std::is_detected<HasLocalFunction,F>::value,
-        "Requires A GridFunction to be passed to the VTKFunction.");
-    }
+    {}
 
-    VTKFunction () = default;
+    VtkFunction () = default;
 
     /// Create a LocalFunction
-    friend VTKLocalFunction<GridView> localFunction (VTKFunction const& self)
+    friend VtkLocalFunction<GridView> localFunction (VtkFunction const& self)
     {
-      return self.asInterface().makeLocalFunction();
+      return self.localFct_;
     }
 
     /// Return a name associated with the function
@@ -220,9 +79,10 @@ namespace Dune { namespace experimental
     }
 
   private:
+    VtkLocalFunction<GridView> localFct_;
     std::string name_;
     int ncomps_ = 1;
-    Vtk::DataTypes type_;
+    Vtk::DataTypes type_ = Vtk::FLOAT32;
   };
 
 }} // end namespace Dune::experimental

dune/vtk/vtklocalfunction.hh

+#pragma once
+
+#include <memory>
+#include <type_traits>
+
+#include <dune/common/std/type_traits.hh>
+
+#include "vtklocalfunctioninterface.hh"
+#include "legacyvtkfunction.hh"
+#include "defaultvtkfunction.hh"
+
+namespace Dune { namespace experimental
+{
+  template <class GridView>
+  class VtkLocalFunction
+  {
+    using Self = VtkLocalFunction;
+    using Entity = typename GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Entity::Geometry::LocalCoordinate;
+
+    template <class LF, class E>
+    using HasBind = decltype(std::declval<LF>().bind(std::declval<E>()));
+
+  public:
+    template <class LF, disableCopyMove<Self, LF> = 0,
+      std::enable_if_t<Std::is_detected<HasBind,LF,Entity>::value, int> = 0>
+    VtkLocalFunction (LF&& lf)
+      : localFct_(std::make_shared<LocalFunctionWrapper<GridView,LF>>(std::forward<LF>(lf)))
+    {}
+
+    VtkLocalFunction (std::shared_ptr<VTKFunction<GridView> const> const& lf)
+      : localFct_(std::make_shared<VTKLocalFunctionWrapper<GridView>>(lf))
+    {}
+
+    VtkLocalFunction () = default;
+
+    /// Bind the function to the grid entity
+    void bind (Entity const& entity)
+    {
+      localFct_->bind(entity);
+    }
+
+    /// Unbind from the currently bound entity
+    void unbind ()
+    {
+      localFct_->unbind();
+    }
+
+    /// Evaluate the `comp` component of the Range value at local coordinate `xi`
+    double evaluate (int comp, LocalCoordinate const& xi) const
+    {
+      return localFct_->evaluate(comp, xi);
+    }
+
+  private:
+    std::shared_ptr<VtkLocalFunctionInterface<GridView>> localFct_;
+  };
+
+}} // end namespace Dune::experimental

dune/vtk/vtklocalfunctioninterface.hh

+#pragma once
+
+namespace Dune { namespace experimental
+{
+  /// An abstract base class for LocalFunctions
+  template <class GridView>
+  class VtkLocalFunctionInterface
+  {
+  public:
+    using Entity = typename GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Entity::Geometry::LocalCoordinate;
+
+    /// Bind the function to the grid entity
+    virtual void bind (Entity const& entity) = 0;
+
+    /// Unbind from the currently bound entity
+    virtual void unbind () = 0;
+
+    /// Evaluate single component comp in the entity at local coordinates xi
+    virtual double evaluate (int comp, LocalCoordinate const& xi) const = 0;
+
+    /// Virtual destructor
+    virtual ~VtkLocalFunctionInterface () = default;
+  };
+
+}} // end namespace Dune::experimental

dune/vtk/writers/vtkwriterinterface.hh

@@ -21,8 +21,7 @@ namespace Dune { namespace experimental
   protected:
     static constexpr int dimension = GridView::dimension;
 
-    using GlobalFunction = VTKFunction<GridView>;
-    using LocalFunction = VTKLocalFunction<GridView>;
+    using VtkFunction = Dune::experimental::VtkFunction<GridView>;
     using pos_type = typename std::ostream::pos_type;
 
     enum PositionTypes {
@@ -47,23 +46,19 @@ namespace Dune { namespace experimental
                 Vtk::FormatTypes format,
                 Vtk::DataTypes datatype = Vtk::FLOAT32);
 
-    /// Attach point data to the writer
-    template <class GridViewFunction>
-    VtkWriterInterface& addPointData (GridViewFunction const& gridViewFct,
-                                      std::string const& name = {},
-                                      int ncomps = 1)
+    /// Attach point data to the writer, \see VtkFunction for possible arguments
+    template <class Function, class... Args>
+    VtkWriterInterface& addPointData (Function const& fct, Args&&... args)
     {
-      pointData_.emplace_back(gridViewFct, name, ncomps);
+      pointData_.emplace_back(fct, std::forward<Args>(args)...);
       return *this;
     }
 
-    /// Attach cell data to the writer
-    template <class GridViewFunction>
-    VtkWriterInterface& addCellData (GridViewFunction const& gridViewFct,
-                                     std::string const& name = {},
-                                     int ncomps = 1)
+    /// Attach cell data to the writer, \see VtkFunction for possible arguments
+    template <class Function, class... Args>
+    VtkWriterInterface& addCellData (Function const& fct, Args&&... args)
     {
-      cellData_.emplace_back(gridViewFct, name, ncomps);
+      cellData_.emplace_back(fct, std::forward<Args>(args)...);
       return *this;
     }
 
@@ -84,14 +79,14 @@ namespace Dune { namespace experimental
     // attributes "offset" in the vector `offsets`.
     void writeData (std::ofstream& out,
                     std::vector<pos_type>& offsets,
-                    GlobalFunction const& fct,
+                    VtkFunction const& fct,
                     PositionTypes type) const;
 
     // Collect point or cell data (depending on \ref PositionTypes) and pass
     // the resulting vector to \ref writeAppended.
     template <class T>
     std::uint64_t writeDataAppended (std::ofstream& out,
-                                     GlobalFunction const& fct,
+                                     VtkFunction const& fct,
                                      PositionTypes type) const;
 
     // Write the coordinates of the vertices to the output stream `out`. In case
@@ -110,7 +105,7 @@ namespace Dune { namespace experimental
     std::uint64_t writeAppended (std::ofstream& out, std::vector<T> const& values) const;
 
     /// Return PointData/CellData attributes for the name of the first scalar/vector/tensor DataArray
-    std::string getNames (std::vector<GlobalFunction> const& data) const
+    std::string getNames (std::vector<VtkFunction> const& data) const
     {
       auto scalar = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 1; });
       auto vector = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 3; });
@@ -135,8 +130,8 @@ namespace Dune { namespace experimental
     Vtk::DataTypes datatype_;
 
     // attached data
-    std::vector<GlobalFunction> pointData_;
-    std::vector<GlobalFunction> cellData_;
+    std::vector<VtkFunction> pointData_;
+    std::vector<VtkFunction> cellData_;
 
     std::size_t const block_size = 1024*32;
     int compression_level = -1; // in [0,9], -1 ... use default value

dune/vtk/writers/vtkwriterinterface.impl.hh

@@ -65,7 +65,7 @@ void VtkWriterInterface<GV,DC>
 template <class GV, class DC>
 void VtkWriterInterface<GV,DC>
   ::writeData (std::ofstream& out, std::vector<pos_type>& offsets,
-               GlobalFunction const& fct, PositionTypes type) const
+               VtkFunction const& fct, PositionTypes type) const
 {
   out << "<DataArray Name=\"" << fct.name() << "\" type=\"" << to_string(fct.type()) << "\""
       << " NumberOfComponents=\"" << fct.ncomps() << "\" format=\"" << (format_ == Vtk::ASCII ? "ascii\">\n" : "appended\"");
@@ -116,7 +116,7 @@ void VtkWriterInterface<GV,DC>
 template <class GV, class DC>
   template <class T>
 std::uint64_t VtkWriterInterface<GV,DC>
-  ::writeDataAppended (std::ofstream& out, GlobalFunction const& fct, PositionTypes type) const
+  ::writeDataAppended (std::ofstream& out, VtkFunction const& fct, PositionTypes type) const
 {
   assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
 

src/CMakeLists.txt

@@ -22,6 +22,10 @@ add_executable("structuredgridwriter" structuredgridwriter.cc)
 target_link_dune_default_libraries("structuredgridwriter")
 target_link_libraries("structuredgridwriter" dunevtk)
 
+add_executable("legacyvtkwriter" legacyvtkwriter.cc)
+target_link_dune_default_libraries("legacyvtkwriter")
+target_link_libraries("legacyvtkwriter" dunevtk)
+
 if (dune-polygongrid_FOUND)
   add_executable("polygongrid" polygongrid.cc)
   target_link_dune_default_libraries("polygongrid")

src/legacyvtkwriter.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <iostream>
+#include <vector>
+
+#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
+#include <dune/common/exceptions.hh> // We use exceptions
+#include <dune/common/filledarray.hh>
+
+#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+
+#include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
+#include <dune/vtk/legacyvtkfunction.hh>
+
+using namespace Dune;
+using namespace Dune::experimental;
+using namespace Dune::Functions;
+
+#define GRID_TYPE 1
+
+int main(int argc, char** argv)
+{
+  Dune::MPIHelper::instance(argc, argv);
+
+  const int dim = 3;
+
+#if GRID_TYPE == 1
+  using GridType = YaspGrid<dim>;
+  FieldVector<double,dim> upperRight; upperRight = 1.0;
+  auto numElements = filledArray<dim,int>(8);
+  GridType grid(upperRight,numElements);
+#elif GRID_TYPE == 2
+  using GridType = UGGrid<dim>;
+  FieldVector<double,dim> lowerLeft; lowerLeft = 0.0;
+  FieldVector<double,dim> upperRight; upperRight = 1.0;
+  auto numElements = filledArray<dim,unsigned int>(4);
+  auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);
+  auto& grid = *gridPtr;
+#endif
+
+  using GridView = typename GridType::LeafGridView;
+  GridView gridView = grid.leafGridView();
+
+  using namespace BasisFactory;
+  auto basis = makeBasis(gridView, lagrange<1>());
+
+  std::vector<double> p1function(basis.dimension());
+
+  interpolate(basis, p1function, [](auto const& x) {
+    return 100*x[0] + 10*x[1] + 1*x[2];
+  });
+
+  using P1Function = P1VTKFunction<GridView,std::vector<double>>;
+  std::shared_ptr<VTKFunction<GridView> const> p1FctWrapped(new P1Function(gridView, p1function, "p1"));
+
+  using Writer = VtkUnstructuredGridWriter<GridView>;
+  Writer vtkWriter(gridView);
+  vtkWriter.addPointData(p1FctWrapped);
+
+  vtkWriter.write("test_ascii_float32.vtu", Vtk::ASCII);
+  vtkWriter.write("test_binary_float32.vtu", Vtk::BINARY);
+  vtkWriter.write("test_compressed_float32.vtu", Vtk::COMPRESSED);
+  vtkWriter.write("test_ascii_float64.vtu", Vtk::ASCII, Vtk::FLOAT64);
+  vtkWriter.write("test_binary_float64.vtu", Vtk::BINARY, Vtk::FLOAT64);
+  vtkWriter.write("test_compressed_float64.vtu", Vtk::COMPRESSED, Vtk::FLOAT64);
+}