diff --git a/dune/vtk/defaultvtkfunction.hh b/dune/vtk/defaultvtkfunction.hh new file mode 100644 index 0000000000000000000000000000000000000000..078909bba5134bf80e0a129f54a41df957c53c2a --- /dev/null +++ b/dune/vtk/defaultvtkfunction.hh @@ -0,0 +1,63 @@ +#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 diff --git a/dune/vtk/legacyvtkfunction.hh b/dune/vtk/legacyvtkfunction.hh new file mode 100644 index 0000000000000000000000000000000000000000..2d3f7979f963b65013c98cbaa8179f166955e115 --- /dev/null +++ b/dune/vtk/legacyvtkfunction.hh @@ -0,0 +1,45 @@ +#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 diff --git a/dune/vtk/vtkfunction.hh b/dune/vtk/vtkfunction.hh index 11a61ece0afd77e03ed6ede88d087a3f157de1c9..08ea87c97a5be89a7f8c96cfff6359bd2ce891b8 100644 --- a/dune/vtk/vtkfunction.hh +++ b/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 diff --git a/dune/vtk/vtklocalfunction.hh b/dune/vtk/vtklocalfunction.hh new file mode 100644 index 0000000000000000000000000000000000000000..9b5dc85927eeee975252c7f9cc80710f2316626c --- /dev/null +++ b/dune/vtk/vtklocalfunction.hh @@ -0,0 +1,59 @@ +#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 diff --git a/dune/vtk/vtklocalfunctioninterface.hh b/dune/vtk/vtklocalfunctioninterface.hh new file mode 100644 index 0000000000000000000000000000000000000000..221dbc1b229f17fff22b0a08a6e188eb23c95bb1 --- /dev/null +++ b/dune/vtk/vtklocalfunctioninterface.hh @@ -0,0 +1,26 @@ +#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 diff --git a/dune/vtk/writers/vtkwriterinterface.hh b/dune/vtk/writers/vtkwriterinterface.hh index caf6e43ca84b70c92e483f720970122a878042d5..26e9202d75901cfa4d87649d3fa1f930c846de7a 100644 --- a/dune/vtk/writers/vtkwriterinterface.hh +++ b/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 diff --git a/dune/vtk/writers/vtkwriterinterface.impl.hh b/dune/vtk/writers/vtkwriterinterface.impl.hh index 5b39fb25fb068424f56e188213817420fecea52d..353fb08c7147d7c4054955e545d12223cdfd435e 100644 --- a/dune/vtk/writers/vtkwriterinterface.impl.hh +++ b/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"); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 23038f9f240071f74ca572c096a6cd58a2f76168..96f90d21ea5c9a0f365b93ed92ba4baa3e59230b 100644 --- a/src/CMakeLists.txt +++ b/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") diff --git a/src/legacyvtkwriter.cc b/src/legacyvtkwriter.cc new file mode 100644 index 0000000000000000000000000000000000000000..8631e0e74c21342ce9f2bf1576d94dacc8174b9f --- /dev/null +++ b/src/legacyvtkwriter.cc @@ -0,0 +1,77 @@ +// -*- 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); +}