diff --git a/README.md b/README.md
index 2dcaa97c2d8e62b772e3d6adc6b3c24416a33b87..c1fa70cb49bcefa200ea3f5da5db1133b6f89dc9 100644
--- a/README.md
+++ b/README.md
@@ -2,14 +2,14 @@
File reader and writer for the VTK Format
## Summary
-Provides structured and unstructured file writers for the VTK XML File Formats
+Provides structured and unstructured file writers for the VTK XML File Formats
that can be opened in the popular ParaView visualization application. Additionally
a file reader is provided to import VTK files into Dune grid and data objects.
## Usage
-The VTK writer works similar to the dune-grid `VTKWriter`. It needs to be bound
+The VTK writer works similar to the dune-grid `VTKWriter`. It needs to be bound
to a GridView and then data can be added to the points or cells in the grid.
-Points are not necessarily grid vertices, but any coordinates placed inside the
+Points are not necessarily grid vertices, but any coordinates placed inside the
grid cells, so the data must be provided as GridViewFunction to allow the local
evaluation in arbitrary local coordinates.
@@ -21,11 +21,11 @@ class Vtk[Type]Writer
public:
// Constructor
Vtk[Type]Writer(GridView, Vtk::FormatTypes = Vtk::BINARY, Vtk::DataTypes = Vtk::FLOAT32);
-
+
// Bind data to the writer
Vtk[Type]Writer& addPointData(Function [, std::string name, int numComponents, Vtk::FormatTypes]);
Vtk[Type]Writer& addCellData(Function [, std::string name, int numComponents, Vtk::FormatTypes]);
-
+
// Write file with filename
void write(std::string filename);
};
@@ -64,74 +64,114 @@ proposed one. A comparions:
| Timeseries | - | **x** |
## Writers and Readers
-Dune-Vtk provides nearly all file formats specified in VTK + 2 time series formats:
+Dune-Vtk provides nearly all file formats specified in VTK + 2 time series formats:
PVD and VTK-Timeseries.
### VtkUnstructuredGridWriter
-Implements a VTK file format for unstructured grids with arbitrary element types
-in 1d, 2d, and 3d. Coordinates are specified explicitly and a connectivity table +
-element types are specified for all grid elements (of codim 0). Can be used with
+Implements a VTK file format for unstructured grids with arbitrary element types
+in 1d, 2d, and 3d. Coordinates are specified explicitly and a connectivity table +
+element types are specified for all grid elements (of codim 0). Can be used with
all Dune grid types.
### VtkStructuredGridWriter
-Implements a writer for grid composed of cube elements (lines, pixels, voxels) with
-local numbering similar to Dunes `cube(d)` numbering. The coordinates of the vertices
-can be arbitrary but the connectivity is implicitly given and equals that of
-`Dune::YaspGrid` or `Dune::SPGrid`. Might be chosen as writer for a transformed
-structured grid, using, e.g., a `GeometryGrid` meta-grid. See `src/geometrygrid.cc`
+Implements a writer for grid composed of cube elements (lines, pixels, voxels) with
+local numbering similar to Dunes `cube(d)` numbering. The coordinates of the vertices
+can be arbitrary but the connectivity is implicitly given and equals that of
+`Dune::YaspGrid` or `Dune::SPGrid`. Might be chosen as writer for a transformed
+structured grid, using, e.g., a `GeometryGrid` meta-grid. See `src/geometrygrid.cc`
for an example.
### VtkRectilinearGridWriter
-Rectilinear grids are tensor-product grids with given coordinates along the x, y,
-and z axes. Therefore, the grid must allow to extract these 1d coordinates and a
-specialization for a `StructuredDataCollector` must be provided, that implements
-the `ordinates()` function. By default, it assumes constant grid spacing starting
-from a lower left corner. For `YaspGrid` a specialization is implemented if the
-coordinates type is `TensorProductCoordinates`. See `src/structuredgridwriter.cc`
+Rectilinear grids are tensor-product grids with given coordinates along the x, y,
+and z axes. Therefore, the grid must allow to extract these 1d coordinates and a
+specialization for a `StructuredDataCollector` must be provided, that implements
+the `ordinates()` function. By default, it assumes constant grid spacing starting
+from a lower left corner. For `YaspGrid` a specialization is implemented if the
+coordinates type is `TensorProductCoordinates`. See `src/structuredgridwriter.cc`
for an example.
### VtkImageDataWriter
-The *most structured* grid is composed of axis-parallel cube elements with constant
-size along each axis. The is implemented in the VtkImageDataWriter. A specialization
-of the `StructuredDataCollector` must implement `origin()` for the lower left corner,
-`wholeExtent()` for the range of cell numbers along each axis in the global grid,
-`extent()` for the range in the local grid, and `spacing()` for the constant grid
+The *most structured* grid is composed of axis-parallel cube elements with constant
+size along each axis. The is implemented in the VtkImageDataWriter. A specialization
+of the `StructuredDataCollector` must implement `origin()` for the lower left corner,
+`wholeExtent()` for the range of cell numbers along each axis in the global grid,
+`extent()` for the range in the local grid, and `spacing()` for the constant grid
spacing in each direction.
### PvdWriter
-A sequence writer, i.e. a collection of timestep files, in the ParaView Data (PVD)
-format. Supports all VtkWriters for the timestep output. In each timestep a collection
+A sequence writer, i.e. a collection of timestep files, in the ParaView Data (PVD)
+format. Supports all VtkWriters for the timestep output. In each timestep a collection
(.pvd) file is created.
### VtkTimseriesWriter
-A timeseries is a collection of timesteps stored in one file, instead of separate
-files for each timestep value. Since in the `Vtk::APPENDED` mode, the data is written
-as binary blocks in the appended section of the file and references by an offset
-in the XML DataArray attributes, it allows to reuse written data. An example of
-usage is when the grid points and cells do not change over time, but just the
+A timeseries is a collection of timesteps stored in one file, instead of separate
+files for each timestep value. Since in the `Vtk::APPENDED` mode, the data is written
+as binary blocks in the appended section of the file and references by an offset
+in the XML DataArray attributes, it allows to reuse written data. An example of
+usage is when the grid points and cells do not change over time, but just the
point-/cell-data. Then, the grid is written only once and the data is just appended.
-Timeseries file are create a bit differently from other Vtk file. There, in the
-first write the grid points and cells are stored in a separate file, and in each
-timestep just the data is written also to temporary files. When you need the timeseries
-file, these stored temporaries are collected and combined to one VTK file. Thus,
-only the minimum amount of data is written in each timestep. The intermediate files
+Timeseries file are create a bit differently from other Vtk file. There, in the
+first write the grid points and cells are stored in a separate file, and in each
+timestep just the data is written also to temporary files. When you need the timeseries
+file, these stored temporaries are collected and combined to one VTK file. Thus,
+only the minimum amount of data is written in each timestep. The intermediate files
are stored, by default, in a `/tmp` folder, with (hopefully) fast write access.
-### VtkReader
-Read in unstructured grid files (.vtu files) and create a new grid, using a GridFactory.
-The reader allows to create the grid in multiple ways, by providing a `GridCreator`
-template parameter. The `ContinuousGridCreator` reads the connectivity of the grid
-as it is and assumes that the elements are already connected correctly. On the other
-hand, a `DiscontinuousGridCreator` reconnects separated elements, by identifying
-matching coordinates of the cell vertices.
+## VtkReader
+Reading unstructured grid files (.vtu files) and creating a new grid, using a GridFactory,
+can be performed using the `VtkReader` class. The reader allows to create the grid in
+multiple ways, by providing a `GridCreator` template parameter. The `ContinuousGridCreator`
+reads the connectivity of the grid as it is and assumes that the elements are already
+connected correctly. On the other hand, a `DiscontinuousGridCreator` reconnects separated
+elements, by identifying matching coordinates of the cell vertices. See more possible
+grid-creators in the directory `dune/vtk/gridcreators`.
-The VtkReader supports grid creation in parallel. If a partition file .pvtu is
-provided, all partitions can be read by 1. one processor and distributed later on
-(`SerialGridCreator`) or read directly in parallel (`ParallelGridCreator`). The later
-is currently only available in dune-alugrid 2.6.
+General interface of a VtkReader
+```c++
+template <class Grid, class GridCreator = ContinuousGridCreator<Grid>, class FieldType = double>
+class VtkReader
+{
+public:
+ // Constructors
+ VtkReader(); // Construct a GridCreator with internal stored GridFactory
+ VtkReader(GridFactory<Grid>&); // Construct a GridCreator referencing the passed GridFactory
+ VtkReader(GridCreator&); // Reference the passed GridCreator
-**TODO:**
+ // Read the data from a file with filename
+ void read(std::string filename);
-- Provide an interface to read the points-/cell-data from the file
-- Extent the implementation to other file formats
+ // Construct the Grid from the data read before.
+ std::unique_ptr<Grid> createGrid() const;
+
+ // Static method to construct the Grid directly
+ static std::unique_ptr<Grid> createGridFromFile(std::string file);
+
+ // Extract data from the reader
+ _PointDataGridFunction_ getPointData(std::string name) const;
+ _CellDataGridFunction_ getCellData(std::string name) const;
+};
+```
+where `Grid` is a dune grid type providing a `GridFactory` specialization, `GridCreator is
+the policy type implementing how the raw data from the file is transformed in data that can
+be passed to the GridFactory and `FieldType` is the data for for internal storage of value
+data from the file, i.e. point-data or cell-data.
+
+The grid can either be created using the static method `createGridFromFile()` or by first
+constructing the `VtkReader` with its `GridCreator`, calling `read()` and finally `createGrid()`.
+The latter allows to access data stored in the reader, like point-data or cell-data.
+
+Value from point-data or cell-data cannot be accessed directly, but through the interface
+of a grid-function. These grid-functions are provided by the `getPointData()` or `getCellData()`
+member functions of the reader. The interface of a dune-functions grid-function concept is
+implemented by these two types. The reason why the reader does not provide the data directly
+is, that it is quiet complicated to associate the specific value to a DOF in the grid, since
+the GridFactory is allows to change the global indexing and even to change to local indexing
+in the elements such that even the local element coordinates might need a transformation
+compared to that of the element stored in the file. All these renumbering and coordinate
+transformations are performed by the grid-functions internally.
+
+The VtkReader supports grid creation in parallel. If a partition file .pvtu is
+provided, all partitions can be read by either one processor and distributed later on
+(`SerialGridCreator`) or read directly in parallel (`ParallelGridCreator`). The later
+is currently only available in dune-alugrid 2.6.
diff --git a/dune/vtk/CMakeLists.txt b/dune/vtk/CMakeLists.txt
index 4715cce1f4e9c2b036dd026a18bdf7ae5134f911..3b2bf43315bcf2242ce2f3ed1f3ba0521d61ab4b 100644
--- a/dune/vtk/CMakeLists.txt
+++ b/dune/vtk/CMakeLists.txt
@@ -27,6 +27,7 @@ install(FILES
add_subdirectory(datacollectors)
add_subdirectory(gridcreators)
+add_subdirectory(gridfunctions)
add_subdirectory(test)
add_subdirectory(utility)
add_subdirectory(writers)
diff --git a/dune/vtk/defaultvtkfunction.hh b/dune/vtk/defaultvtkfunction.hh
index 2dd750bf9167748c21be6a25ad0d7c28e0c31403..38706d3e6fa3891e8a478f8c9cbf5559ecdd608b 100644
--- a/dune/vtk/defaultvtkfunction.hh
+++ b/dune/vtk/defaultvtkfunction.hh
@@ -67,13 +67,22 @@ namespace Dune
return comp < N ? vec[comp] : 0.0;
}
- // Return the scalar values
- template <class T>
- double evaluateImpl (int comp, T const& value) const
- {
- assert(comp == 0);
- return value;
- }
+ // Evaluate a component of a vector valued data
+ template <class T,
+ std::enable_if_t<IsIndexable<T,int>::value, int> = 0>
+ double evaluateImpl (int comp, T const& value) const
+ {
+ return value[comp];
+ }
+
+ // Return the scalar values
+ template <class T,
+ std::enable_if_t<not IsIndexable<T,int>::value, int> = 0>
+ double evaluateImpl (int comp, T const& value) const
+ {
+ assert(comp == 0);
+ return value;
+ }
private:
LocalFunction localFct_;
diff --git a/dune/vtk/gridcreatorinterface.hh b/dune/vtk/gridcreatorinterface.hh
index a92e88bdcf637ac426ee7920953dce56f16d6e78..68bd18fe1934635a5cdf671893cc40dd4ea21637 100644
--- a/dune/vtk/gridcreatorinterface.hh
+++ b/dune/vtk/gridcreatorinterface.hh
@@ -29,9 +29,21 @@ namespace Dune
using Derived = DerivedType;
public:
- /// Constructor. Stores a reference to the passed GridFactory
+ /// Constructor. Stores a reference to the passed GridFactory.
GridCreatorInterface (GridFactory<Grid>& factory)
- : factory_(&factory)
+ : factory_(stackobject_to_shared_ptr(factory))
+ {}
+
+ /// Constructor. Store the shared_ptr to the GridFactory.
+ GridCreatorInterface (std::shared_ptr<GridFactory<Grid>> factory)
+ : factory_(std::move(factory))
+ {}
+
+ /// Constructor. Construct a new GridFactory from the passed arguments.
+ template <class... Args,
+ std::enable_if_t<std::is_constructible<GridFactory<Grid>, Args...>::value,int> = 0>
+ GridCreatorInterface (Args&&... args)
+ : factory_(std::make_shared<GridFactory<Grid>>(std::forward<Args>(args)...))
{}
/// Insert all points as vertices into the factory
@@ -112,7 +124,7 @@ namespace Dune
}
protected:
- GridFactory<Grid>* factory_;
+ std::shared_ptr<GridFactory<Grid>> factory_;
};
} // end namespace Vtk
diff --git a/dune/vtk/gridcreators/continuousgridcreator.hh b/dune/vtk/gridcreators/continuousgridcreator.hh
index b259132aac0ece102538ab11376233fafbc082dc..555161bef18450f2a626e7e2c17640516c3876e1 100644
--- a/dune/vtk/gridcreators/continuousgridcreator.hh
+++ b/dune/vtk/gridcreators/continuousgridcreator.hh
@@ -11,6 +11,7 @@
#include <dune/vtk/types.hh>
#include <dune/vtk/gridcreatorinterface.hh>
+#include <dune/vtk/gridfunctions/continuousgridfunction.hh>
namespace Dune
{
@@ -26,10 +27,7 @@ namespace Dune
using Super = GridCreatorInterface<Grid, ContinuousGridCreator>;
using GlobalCoordinate = typename Super::GlobalCoordinate;
- ContinuousGridCreator (GridFactory<Grid>& factory)
- : Super(factory)
- {}
-
+ using Super::Super;
using Super::factory;
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
@@ -69,5 +67,16 @@ namespace Dune
}
};
+ // deduction guides
+ template <class Grid>
+ ContinuousGridCreator(GridFactory<Grid>&)
+ -> ContinuousGridCreator<Grid>;
+
+ template <class GridType, class FieldType, class Context>
+ struct AssociatedGridFunction<ContinuousGridCreator<GridType>, FieldType, Context>
+ {
+ using type = ContinuousGridFunction<GridType, FieldType, Context>;
+ };
+
} // end namespace Vtk
} // end namespace Dune
diff --git a/dune/vtk/gridcreators/derivedgridcreator.hh b/dune/vtk/gridcreators/derivedgridcreator.hh
index 946aee644d206c991ad2c542085b4e0915fccd4a..ca0480d8e344ba04b8120544963765c92d9924a3 100644
--- a/dune/vtk/gridcreators/derivedgridcreator.hh
+++ b/dune/vtk/gridcreators/derivedgridcreator.hh
@@ -22,9 +22,11 @@ namespace Dune
using Grid = typename GridCreator::Grid;
using GlobalCoordinate = typename Super::GlobalCoordinate;
- DerivedGridCreator (GridFactory<Grid>& factory)
- : Super(factory)
- , gridCreator_(factory)
+ template <class... Args,
+ disableCopyMove<DerivedGridCreator, Args...> = 0>
+ DerivedGridCreator (Args&&... args)
+ : Super(std::forward<Args>(args)...)
+ , gridCreator_(Super::factory())
{}
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
diff --git a/dune/vtk/gridcreators/discontinuousgridcreator.hh b/dune/vtk/gridcreators/discontinuousgridcreator.hh
index ca4d29c70b95d208d5ce06da8ab63c1716bcdc1e..4ea2e89f4998130732b742c7200df0125f574fc2 100644
--- a/dune/vtk/gridcreators/discontinuousgridcreator.hh
+++ b/dune/vtk/gridcreators/discontinuousgridcreator.hh
@@ -40,11 +40,9 @@ namespace Dune
}
};
- DiscontinuousGridCreator (GridFactory<Grid>& factory)
- : Super(factory)
- {}
-
+ using Super::Super;
using Super::factory;
+
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
std::vector<std::uint64_t> const& /*point_ids*/)
{
@@ -97,5 +95,10 @@ namespace Dune
std::map<GlobalCoordinate, std::size_t, CoordLess> uniquePoints_;
};
+ // deduction guides
+ template <class Grid>
+ DiscontinuousGridCreator(GridFactory<Grid>&)
+ -> DiscontinuousGridCreator<Grid>;
+
} // end namespace Vtk
} // end namespace Dune
diff --git a/dune/vtk/gridcreators/lagrangegridcreator.hh b/dune/vtk/gridcreators/lagrangegridcreator.hh
index bbe6010b3ad6c321eed8ded565680c67b69bb318..b5daeb0abe182c4f97593e24427094daf440c275 100644
--- a/dune/vtk/gridcreators/lagrangegridcreator.hh
+++ b/dune/vtk/gridcreators/lagrangegridcreator.hh
@@ -15,6 +15,7 @@
#include <dune/vtk/types.hh>
#include <dune/vtk/gridcreatorinterface.hh>
+#include <dune/vtk/gridfunctions/lagrangegridfunction.hh>
#include <dune/vtk/utility/lagrangepoints.hh>
namespace Dune
@@ -58,11 +59,11 @@ namespace Dune
class LocalFunction;
public:
- using Super::factory;
+ using LocalGeometry = MultiLinearGeometry<typename Element::Geometry::ctype,Element::dimension,Element::dimension>;
- LagrangeGridCreator (GridFactory<GridType>& factory)
- : Super(factory)
- {}
+ public:
+ using Super::Super;
+ using Super::factory;
/// Implementation of the interface function `insertVertices()`
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
@@ -163,22 +164,41 @@ namespace Dune
* The returned LocalParametrization is a mapping `GlobalCoordinate(LocalCoordinate)`
* where `LocalCoordinate is w.r.t. the local coordinate system in the passed element
* and `GlobalCoordinate` a world coordinate in the parametrized grid.
- *
- * Note, when an element is passed, it might have a different local coordinate system
- * than the coordinate system used to defined the element parametrization. Thus
- * coordinate transform is internally chained to the evaluation of the local
- * parametrization. This local geometry transform is obtained by figuring out the
- * permutation of corners in the element corresponding to the inserted corner
- * vertices.
**/
LocalParametrization localParametrization (Element const& element) const
{
- assert(!nodes_.empty() && !parametrization_.empty());
+ VTK_ASSERT(!nodes_.empty() && !parametrization_.empty());
unsigned int insertionIndex = factory().insertionIndex(element);
auto const& localParam = parametrization_.at(insertionIndex);
- assert(element.type() == localParam.type);
+ VTK_ASSERT(element.type() == localParam.type);
+
+ return {nodes_, localParam, order(localParam), localGeometry(element, localParam)};
+ }
+
+ /// \brief Construct a transformation of local element coordinates
+ /**
+ * An element might have a different local coordinate system than the coordinate
+ * system used to defined the element parametrization. Thus coordinate transform
+ * of the local parametrization is needed for element-local evaluations. This
+ * local geometry transform is obtained by figuring out the permutation of corners
+ * in the element corresponding to the inserted corner vertices.
+ **/
+ LocalGeometry localGeometry (Element const& element) const
+ {
+ VTK_ASSERT(!nodes_.empty() && !parametrization_.empty());
+
+ unsigned int insertionIndex = factory().insertionIndex(element);
+ auto const& localParam = parametrization_.at(insertionIndex);
+ VTK_ASSERT(element.type() == localParam.type);
+
+ return localGeometry(element, localParam);
+ }
+ private:
+ // implementation details of localGeometry()
+ LocalGeometry localGeometry (Element const& element, ElementParametrization const& localParam) const
+ {
// collect indices of vertices
std::vector<unsigned int> indices(element.subEntities(GridType::dimension));
for (unsigned int i = 0; i < element.subEntities(GridType::dimension); ++i)
@@ -188,31 +208,41 @@ namespace Dune
std::vector<unsigned int> permutation(indices.size());
for (std::size_t i = 0; i < indices.size(); ++i) {
auto it = std::find(localParam.corners.begin(), localParam.corners.end(), indices[i]);
- assert(it != localParam.corners.end());
+ VTK_ASSERT(it != localParam.corners.end());
permutation[i] = std::distance(localParam.corners.begin(), it);
}
- return LocalParametrization{nodes_, localParam, order(localParam), permutation};
+ auto refElem = referenceElement<typename Element::Geometry::ctype,Element::dimension>(localParam.type);
+ std::vector<LocalCoordinate> corners(permutation.size());
+ for (std::size_t i = 0; i < permutation.size(); ++i)
+ corners[i] = refElem.position(permutation[i], Element::dimension);
+
+ return {localParam.type, corners};
}
+ public:
+
/// Determine lagrange order from number of points
- template <class LocalParam>
- int order (LocalParam const& localParam) const
+ int order (GeometryType type, std::size_t nNodes) const
{
- GeometryType type = localParam.type;
- int nNodes = localParam.nodes.size();
- for (int o = 1; o <= nNodes; ++o)
+ for (int o = 1; o <= int(nNodes); ++o)
if (numLagrangePoints(type.id(), type.dim(), o) == std::size_t(nNodes))
return o;
return 1;
}
+ int order (ElementParametrization const& localParam) const
+ {
+ return order(localParam.type, localParam.nodes.size());
+ }
+
/// Determine lagrange order from number of points from the first element parametrization
int order () const
{
assert(!parametrization_.empty());
- return order(parametrization_.front());
+ auto const& localParam = parametrization_.front();
+ return order(localParam);
}
public:
@@ -274,6 +304,16 @@ namespace Dune
Parametrization parametrization_;
};
+ // deduction guides
+ template <class Grid>
+ LagrangeGridCreator(GridFactory<Grid>&)
+ -> LagrangeGridCreator<Grid>;
+
+ template <class GridType, class FieldType, class Context>
+ struct AssociatedGridFunction<LagrangeGridCreator<GridType>, FieldType, Context>
+ {
+ using type = LagrangeGridFunction<GridType, FieldType, Context>;
+ };
template <class Grid>
class LagrangeGridCreator<Grid>::LocalParametrization
@@ -300,16 +340,11 @@ namespace Dune
}
/// Construct a local element parametrization for elements with permuted corners
- template <class Nodes, class LocalParam, class Permutation>
- LocalParametrization (Nodes const& nodes, LocalParam const& param, int order, Permutation const& permutation)
+ template <class Nodes, class LocalParam, class LG>
+ LocalParametrization (Nodes const& nodes, LocalParam const& param, int order, LG&& localGeometry)
: LocalParametrization(nodes, param, order)
{
- auto refElem = referenceElement<ctype,Grid::dimension>(param.type);
- std::vector<LocalCoordinate> corners(permutation.size());
- for (std::size_t i = 0; i < permutation.size(); ++i)
- corners[i] = refElem.position(permutation[i], Grid::dimension);
-
- localGeometry_.emplace(param.type, corners);
+ localGeometry_.emplace(std::forward<LG>(localGeometry));
}
/// Evaluate the local parametrization in local coordinates
diff --git a/dune/vtk/gridcreators/parallelgridcreator.hh b/dune/vtk/gridcreators/parallelgridcreator.hh
index 164a09e8e6a99f58a1d6e2fb10fae08b73934cb4..5f10058091685add56750d974d8363a79c8fe5a5 100644
--- a/dune/vtk/gridcreators/parallelgridcreator.hh
+++ b/dune/vtk/gridcreators/parallelgridcreator.hh
@@ -27,9 +27,9 @@ namespace Dune
// The GridFactory must support insertion of global vertex IDs
static_assert(Std::is_detected<HasInsertVertex, GridFactory<Grid>, GlobalCoordinate, VertexId>{}, "");
- ParallelGridCreator (GridFactory<Grid>& factory)
- : Super(factory)
- {}
+ public:
+
+ using Super::Super;
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
std::vector<std::uint64_t> const& point_ids)
@@ -48,5 +48,10 @@ namespace Dune
}
};
+ // deduction guides
+ template <class Grid>
+ ParallelGridCreator(GridFactory<Grid>&)
+ -> ParallelGridCreator<Grid>;
+
} // end namespace Vtk
} // end namespace Dune
diff --git a/dune/vtk/gridcreators/serialgridcreator.hh b/dune/vtk/gridcreators/serialgridcreator.hh
index de9814954357e4380d12d283cc94992f4972bb46..369f4b9b4e1a306183aeb013c46d24e3c302d2c2 100644
--- a/dune/vtk/gridcreators/serialgridcreator.hh
+++ b/dune/vtk/gridcreators/serialgridcreator.hh
@@ -21,9 +21,9 @@ namespace Dune
using Super = GridCreatorInterface<Grid, Self>;
using GlobalCoordinate = typename Super::GlobalCoordinate;
- SerialGridCreator (GridFactory<Grid>& factory)
- : Super(factory)
- {}
+ public:
+
+ using Super::Super;
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
std::vector<std::uint64_t> const& /*point_ids*/)
@@ -72,5 +72,10 @@ namespace Dune
std::vector<std::int64_t> shift_;
};
+ // deduction guides
+ template <class Grid>
+ SerialGridCreator(GridFactory<Grid>&)
+ -> SerialGridCreator<Grid>;
+
} // end namespace Vtk
} // end namespace Dune
diff --git a/dune/vtk/gridfunctions/CMakeLists.txt b/dune/vtk/gridfunctions/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f4d7d41e5595aeaaf0317e183a0d517949fd1196
--- /dev/null
+++ b/dune/vtk/gridfunctions/CMakeLists.txt
@@ -0,0 +1,6 @@
+#install headers
+install(FILES
+ common.hh
+ continuousgridfunction.hh
+ lagrangegridfunction.hh
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vtk/gridfunctions)
diff --git a/dune/vtk/gridfunctions/common.hh b/dune/vtk/gridfunctions/common.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2a57b5f63a0a373abc0ed258123de804fbabdb9a
--- /dev/null
+++ b/dune/vtk/gridfunctions/common.hh
@@ -0,0 +1,31 @@
+#pragma once
+
+namespace Dune
+{
+ namespace Vtk
+ {
+ /// Context indicating that a GridFunction generates a local-function from point data
+ struct PointContext {};
+
+ /// Context indicating that a GridFunction generates a local-function from cell data
+ struct CellContext {};
+
+ /// \brief Type-Traits to associate a GridFunction to a GridCreator.
+ /**
+ * Each GridCreator type should specialize this template and set `type` to the
+ * corresponding GridFunction type, e.g. Vtk::ContinuousGridFunction or
+ * Vtk::LagrangeGridFunction.
+ *
+ * \tparam GridCreator A Type implementing the GridCreatorInterface
+ * \tparam FieldType Coefficient type of the data extracted from the file.
+ * \tparam Context A context-type for specialization of the local-function, e.g.,
+ * Vtk::PointContext or Vtk::CellContext.
+ **/
+ template <class GridCreator, class FieldType, class Context>
+ struct AssociatedGridFunction
+ {
+ using type = void;
+ };
+
+ } // end namespace Vtk
+} // end namespace Dune
diff --git a/dune/vtk/gridfunctions/continuousgridfunction.hh b/dune/vtk/gridfunctions/continuousgridfunction.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8d9c6fa4a7d2859d75bb868782eca47a0ba2293d
--- /dev/null
+++ b/dune/vtk/gridfunctions/continuousgridfunction.hh
@@ -0,0 +1,193 @@
+#pragma once
+
+#include <vector>
+
+#include <dune/common/dynvector.hh>
+#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
+#include <dune/vtk/gridfunctions/common.hh>
+
+namespace Dune
+{
+ namespace Vtk
+ {
+ /// \brief A GridFunction representing data stored on the grid vertices in a continuous manner.
+ template <class GridType, class FieldType, class Context>
+ class ContinuousGridFunction
+ {
+ using Grid = GridType;
+ using Field = FieldType;
+
+ using Factory = GridFactory<Grid>;
+
+ public:
+ struct EntitySet
+ {
+ using Grid = GridType;
+ using Element = typename GridType::template Codim<0>::Entity;
+ using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+ using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
+ };
+
+ using Domain = typename EntitySet::GlobalCoordinate;
+ using Range = DynamicVector<Field>;
+ using Signature = Range(Domain);
+
+ private:
+ template <class LC>
+ class PointDataLocalFunction
+ {
+ // NOTE: local finite-element fixed to Lagrange P1/Q1
+ using LFECache = LagrangeLocalFiniteElementCache<Field,Field,LC::mydimension,1>;
+ using LFE = typename LFECache::FiniteElementType;
+ using LB = typename LFE::Traits::LocalBasisType;
+
+ public:
+ using LocalContext = LC;
+ using Domain = typename LC::Geometry::LocalCoordinate;
+ using Range = DynamicVector<Field>;
+ using Signature = Range(Domain);
+
+ public:
+ PointDataLocalFunction (Factory const& factory, std::vector<Field> const& values, unsigned int comp)
+ : factory_(factory)
+ , values_(values)
+ , comp_(comp)
+ , cache_()
+ {}
+
+ void bind (LocalContext const& element)
+ {
+ lfe_ = &cache_.get(element.type());
+
+ // collect values on vertices
+ // NOTE: assumes, that Lagrange nodes are ordered like element vertices
+ localValues_.resize(element.subEntities(Grid::dimension));
+ for (unsigned int i = 0; i < element.subEntities(Grid::dimension); ++i) {
+ unsigned int idx = factory_.insertionIndex(element.template subEntity<Grid::dimension>(i));
+ DynamicVector<Field>& v = localValues_[i];
+ v.resize(comp_);
+ for (unsigned int j = 0; j < comp_; ++j)
+ v[j] = values_[comp_*idx + j];
+ }
+ }
+
+ void unbind ()
+ {
+ lfe_ = nullptr;
+ }
+
+ Range operator() (Domain const& local) const
+ {
+ assert(!!lfe_);
+ auto const& lb = lfe_->localBasis();
+ lb.evaluateFunction(local, shapeValues_);
+ assert(shapeValues_.size() == localValues_.size());
+
+ Range y(comp_, Field(0));
+ for (std::size_t i = 0; i < shapeValues_.size(); ++i)
+ y.axpy(shapeValues_[i], localValues_[i]);
+
+ return y;
+ }
+
+ private:
+ Factory const& factory_;
+ std::vector<Field> const& values_;
+ unsigned int comp_;
+ LFECache cache_;
+
+ // Local Finite-Element
+ LFE const* lfe_ = nullptr;
+
+ // cache of local values
+ std::vector<DynamicVector<Field>> localValues_;
+ mutable std::vector<typename LB::Traits::RangeType> shapeValues_;
+ };
+
+ template <class LC>
+ class CellDataLocalFunction
+ {
+ public:
+ using LocalContext = LC;
+ using Domain = typename LC::Geometry::LocalCoordinate;
+ using Range = DynamicVector<Field>;
+ using Signature = Range(Domain);
+
+ public:
+ CellDataLocalFunction (Factory const& factory, std::vector<Field> const& values, unsigned int comp)
+ : factory_(factory)
+ , values_(values)
+ , comp_(comp)
+ {}
+
+ void bind (LocalContext const& element)
+ {
+ unsigned int idx = factory_.insertionIndex(element);
+
+ // collect values on cells
+ DynamicVector<Field>& v = localValue_;
+ v.resize(comp_);
+
+ for (unsigned int j = 0; j < comp_; ++j)
+ v[j] = values_[comp_*idx + j];
+ }
+
+ void unbind ()
+ {}
+
+ Range operator() (Domain const& local) const
+ {
+ return localValue_;
+ }
+
+ private:
+ Factory const& factory_;
+ std::vector<Field> const& values_;
+ unsigned int comp_;
+
+ // cache of local values
+ DynamicVector<Field> localValue_;
+ };
+
+ template <class LC>
+ using LocalFunction = std::conditional_t< std::is_same_v<Context,PointContext>,
+ PointDataLocalFunction<LC>,
+ CellDataLocalFunction<LC>>;
+
+ public:
+ template <class GridCreator>
+ ContinuousGridFunction (GridCreator const& creator, std::vector<Field> const& values, unsigned int comp,
+ std::vector<std::uint8_t> const& /*types*/,
+ std::vector<std::int64_t> const& /*offsets*/,
+ std::vector<std::int64_t> const& /*connectivity*/)
+ : factory_(creator.factory())
+ , values_(values)
+ , comp_(comp)
+ {}
+
+ Range operator() (Domain const& global) const
+ {
+ DUNE_THROW(Dune::NotImplemented, "Evaluation in global coordinates not implemented.");
+ return Range(comp_, 0);
+ }
+
+ EntitySet const& entitySet () const
+ {
+ return entitySet_;
+ }
+
+ friend LocalFunction<typename EntitySet::Element> localFunction (ContinuousGridFunction const& gf)
+ {
+ return {gf.factory_, gf.values_, gf.comp_};
+ }
+
+ private:
+ Factory const& factory_;
+ std::vector<Field> const& values_;
+ unsigned int comp_;
+
+ EntitySet entitySet_;
+ };
+
+ } // end namespace Vtk
+} // end namespace Dune
diff --git a/dune/vtk/gridfunctions/lagrangegridfunction.hh b/dune/vtk/gridfunctions/lagrangegridfunction.hh
new file mode 100644
index 0000000000000000000000000000000000000000..15d4f69b8e04228305e6ec04741950a63dcc0af5
--- /dev/null
+++ b/dune/vtk/gridfunctions/lagrangegridfunction.hh
@@ -0,0 +1,255 @@
+#pragma once
+
+#include <optional>
+#include <vector>
+
+#include <dune/common/dynvector.hh>
+#include <dune/localfunctions/lagrange.hh>
+#include <dune/vtk/gridfunctions/common.hh>
+#include <dune/vtk/gridfunctions/continuousgridfunction.hh>
+#include <dune/vtk/utility/errors.hh>
+#include <dune/vtk/utility/lagrangepoints.hh>
+
+namespace Dune
+{
+ namespace Vtk
+ {
+ template <class GridType>
+ class LagrangeGridCreator;
+
+ /// \brief Grid-function representing values from a VTK file with local Lagrange
+ /// interpolation of the values stored on the Lagrange nodes.
+ template <class GridType, class FieldType, class Context>
+ class LagrangeGridFunction
+ {
+ using Grid = GridType;
+ using Field = FieldType;
+
+ using Factory = GridFactory<Grid>;
+ using GridCreator = LagrangeGridCreator<Grid>;
+
+ public:
+ struct EntitySet
+ {
+ using Grid = GridType;
+ using Element = typename GridType::template Codim<0>::Entity;
+ using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+ using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
+ };
+
+ using Domain = typename EntitySet::GlobalCoordinate;
+ using Range = DynamicVector<Field>;
+ using Signature = Range(Domain);
+
+ private:
+ template <class LC>
+ class PointDataLocalFunction
+ {
+ using LFE = LagrangeLocalFiniteElement<LagrangePointSet, LC::dimension, FieldType, FieldType>;
+ using LB = typename LFE::Traits::LocalBasisType;
+
+ public:
+ using LocalContext = LC;
+ using Domain = typename LC::Geometry::LocalCoordinate;
+ using Range = DynamicVector<Field>;
+ using Signature = Range(Domain);
+
+ public:
+ /// Constructor. Stores references to the passed data.
+ PointDataLocalFunction (GridCreator const& creator, std::vector<Field> const& values, unsigned int comp,
+ std::vector<std::uint8_t> const& types,
+ std::vector<std::int64_t> const& offsets,
+ std::vector<std::int64_t> const& connectivity)
+ : creator_(creator)
+ , values_(values)
+ , comp_(comp)
+ , types_(types)
+ , offsets_(offsets)
+ , connectivity_(connectivity)
+ {}
+
+ /// Binding the local-function to an element.
+ /**
+ * Constructs a new local finite-element with a polynomial order given
+ * by the number of Lagrange nodes on the element. Extracts values on all
+ * Lagrange nodes stored in file and stores these local coefficients in
+ * a class member variable.
+ **/
+ void bind (LocalContext const& element)
+ {
+ unsigned int insertionIndex = creator_.factory().insertionIndex(element);
+
+ std::int64_t shift = (insertionIndex == 0 ? 0 : offsets_[insertionIndex-1]);
+ std::int64_t numNodes = offsets_[insertionIndex] - shift;
+ [[maybe_unused]] std::int64_t maxNumNodes = numLagrangePoints(element.type().id(), element.type().dim(), 20);
+ VTK_ASSERT(numNodes > 0 && numNodes < maxNumNodes);
+
+ int order = creator_.order(element.type(), numNodes);
+ VTK_ASSERT(order > 0 && order < 20);
+
+ // construct a local finite-element with the corresponding order and Lagrange points
+ // as stored in the file
+ lfe_.emplace(LFE{element.type(), (unsigned int)(order)});
+ lgeo_.emplace(creator_.localGeometry(element));
+
+ // collect values on lagrange nodes
+ localValues_.resize(numNodes);
+ for (std::int64_t i = shift, i0 = 0; i < offsets_[insertionIndex]; ++i, ++i0) {
+ std::int64_t idx = connectivity_[i];
+ DynamicVector<Field>& v = localValues_[i0];
+ v.resize(comp_);
+ for (unsigned int j = 0; j < comp_; ++j)
+ v[j] = values_[comp_*idx + j];
+ }
+ }
+
+ /// Unbind the local-function and the local finite-element from the element.
+ void unbind ()
+ {
+ lfe_.reset();
+ lgeo_.reset();
+ }
+
+ /// Evaluation in element local coordinates. Essentially, a local Lagrange
+ /// interpolation with coefficients extracted in \ref bind().
+ // NOTE: do we need to transform the local coordinates?
+ Range operator() (Domain const& local) const
+ {
+ assert(!!lfe_);
+ auto const& lb = lfe_->localBasis();
+ lb.evaluateFunction(lgeo_->global(local), shapeValues_);
+ assert(shapeValues_.size() == localValues_.size());
+
+ Range y(comp_, Field(0));
+ for (std::size_t i = 0; i < shapeValues_.size(); ++i)
+ y.axpy(shapeValues_[i], localValues_[i]);
+
+ return y;
+ }
+
+ private:
+ GridCreator const& creator_;
+ std::vector<Field> const& values_;
+ unsigned int comp_;
+ std::vector<std::uint8_t> const& types_;
+ std::vector<std::int64_t> const& offsets_;
+ std::vector<std::int64_t> const& connectivity_;
+
+ // Local Finite-Element
+ std::optional<LFE> lfe_ = std::nullopt;
+ std::optional<typename GridCreator::LocalGeometry> lgeo_ = std::nullopt;
+
+ // cache of local values
+ std::vector<DynamicVector<Field>> localValues_;
+ mutable std::vector<typename LB::Traits::RangeType> shapeValues_;
+ };
+
+ /// Evaluation of data on the cells of the grid
+ template <class LC>
+ class CellDataLocalFunction
+ {
+ public:
+ using LocalContext = LC;
+ using Domain = typename LC::Geometry::LocalCoordinate;
+ using Range = DynamicVector<Field>;
+ using Signature = Range(Domain);
+
+ public:
+ /// Constructor. Stores references to the passed data.
+ CellDataLocalFunction (GridCreator const& creator, std::vector<Field> const& values, unsigned int comp,
+ std::vector<std::uint8_t> const& /*types*/,
+ std::vector<std::int64_t> const& /*offsets*/,
+ std::vector<std::int64_t> const& /*connectivity*/)
+ : factory_(creator.factory())
+ , values_(values)
+ , comp_(comp)
+ {}
+
+ /// Binding the local-function to an element extract the cell-value from the vector
+ /// of data.
+ void bind (LocalContext const& element)
+ {
+ unsigned int idx = factory_.insertionIndex(element);
+
+ // collect values on cells
+ DynamicVector<Field>& v = localValue_;
+ v.resize(comp_);
+
+ for (unsigned int j = 0; j < comp_; ++j)
+ v[j] = values_[comp_*idx + j];
+ }
+
+ /// Unbinds from the bound element. Does nothing
+ void unbind ()
+ {}
+
+ /// Evaluation in local element coordinates. Returns the constant value
+ /// extracted in \ref bind().
+ Range operator() (Domain const& local) const
+ {
+ return localValue_;
+ }
+
+ private:
+ Factory const& factory_;
+ std::vector<Field> const& values_;
+ unsigned int comp_;
+
+ // cache of local values
+ DynamicVector<Field> localValue_;
+ };
+
+ /// Type switch for local-functions depending on the Context
+ template <class LC>
+ using LocalFunction = std::conditional_t< std::is_same_v<Context,PointContext>,
+ PointDataLocalFunction<LC>,
+ CellDataLocalFunction<LC>>;
+
+ public:
+ /// Construct a grid-function. Passed in data is stroed by reference, thus must have
+ /// a life-time greater than that of the grid-function and corresponding local-function.
+ LagrangeGridFunction (GridCreator const& creator, std::vector<Field> const& values, unsigned int comp,
+ std::vector<std::uint8_t> const& types,
+ std::vector<std::int64_t> const& offsets,
+ std::vector<std::int64_t> const& connectivity)
+ : creator_(creator)
+ , values_(values)
+ , comp_(comp)
+ , types_(types)
+ , offsets_(offsets)
+ , connectivity_(connectivity)
+ {}
+
+ /// Global evaluation. Not supported!
+ Range operator() (Domain const& global) const
+ {
+ DUNE_THROW(Dune::NotImplemented, "Evaluation in global coordinates not implemented.");
+ return Range(comp_, 0);
+ }
+
+ /// Return a type that defines the element that can be iterated.
+ EntitySet const& entitySet () const
+ {
+ return entitySet_;
+ }
+
+ /// Construct a local-function depending on the Context type either PointDataLocalFunction
+ /// or CellDataLocalFunction
+ friend LocalFunction<typename EntitySet::Element> localFunction (LagrangeGridFunction const& gf)
+ {
+ return {gf.creator_, gf.values_, gf.comp_, gf.types_, gf.offsets_, gf.connectivity_};
+ }
+
+ private:
+ GridCreator const& creator_;
+ std::vector<Field> const& values_;
+ unsigned int comp_;
+ std::vector<std::uint8_t> const& types_;
+ std::vector<std::int64_t> const& offsets_;
+ std::vector<std::int64_t> const& connectivity_;
+
+ EntitySet entitySet_;
+ };
+
+ } // end namespace Vtk
+} // end namespace Dune
diff --git a/dune/vtk/utility/lagrangepoints.impl.hh b/dune/vtk/utility/lagrangepoints.impl.hh
index e49934513c16e2161bdebad4141c597e9e7b341b..6d3f5c4210ee5e314887a1e21025b5a2f3241d49 100644
--- a/dune/vtk/utility/lagrangepoints.impl.hh
+++ b/dune/vtk/utility/lagrangepoints.impl.hh
@@ -136,7 +136,7 @@ void LagrangePointSetBuilder<K,2>::buildTriangle (std::size_t nPoints, int order
std::array<int,3> bindex;
- double order_d = double(order);
+ K order_d = K(order);
for (std::size_t p = 0; p < nPoints; ++p) {
barycentricIndex(p, bindex, order);
Vec point{bindex[0] / order_d, bindex[1] / order_d};
@@ -198,8 +198,8 @@ void LagrangePointSetBuilder<K,2>::buildQuad(std::size_t nPoints, int order, Poi
for (int m = 0; m <= orders[0]; ++m) {
// int idx = pointIndexFromIJ(m,n,orders);
- const double r = double(m) / orders[0];
- const double s = double(n) / orders[1];
+ const K r = K(m) / orders[0];
+ const K s = K(n) / orders[1];
Vec p = (1.0 - r) * (nodes[3] * s + nodes[0] * (1.0 - s)) +
r * (nodes[2] * s + nodes[1] * (1.0 - s));
@@ -317,7 +317,7 @@ void LagrangePointSetBuilder<K,3>::buildTetra (std::size_t nPoints, int order, P
std::array<int,4> bindex;
- double order_d = double(order);
+ K order_d = K(order);
for (std::size_t p = 0; p < nPoints; ++p) {
barycentricIndex(p, bindex, order);
Vec point{bindex[0] / order_d, bindex[1] / order_d, bindex[2] / order_d};
@@ -407,9 +407,9 @@ void LagrangePointSetBuilder<K,3>::buildHex (std::size_t nPoints, int order, Poi
for (int p = 0; p <= orders[2]; ++p) {
for (int n = 0; n <= orders[1]; ++n) {
for (int m = 0; m <= orders[0]; ++m) {
- const double r = double(m) / orders[0];
- const double s = double(n) / orders[1];
- const double t = double(p) / orders[2];
+ const K r = K(m) / orders[0];
+ const K s = K(n) / orders[1];
+ const K t = K(p) / orders[2];
Vec point = (1.0-r) * ((nodes[3] * (1.0-t) + nodes[7] * t) * s + (nodes[0] * (1.0-t) + nodes[4] * t) * (1.0-s)) +
r * ((nodes[2] * (1.0-t) + nodes[6] * t) * s + (nodes[1] * (1.0-t) + nodes[5] * t) * (1.0-s));
diff --git a/dune/vtk/vtkreader.hh b/dune/vtk/vtkreader.hh
index e5f683a92a8e09960e66d2ad070a76dddab664b4..2eb6f942a94abd0cfdcd2fc8d74494be8cb991ac 100644
--- a/dune/vtk/vtkreader.hh
+++ b/dune/vtk/vtkreader.hh
@@ -10,9 +10,11 @@
#include <dune/vtk/filereader.hh>
#include <dune/vtk/types.hh>
+#include <dune/vtk/utility/errors.hh>
// default GridCreator
#include <dune/vtk/gridcreators/continuousgridcreator.hh>
+#include <dune/vtk/gridfunctions/common.hh>
namespace Dune
{
@@ -21,13 +23,14 @@ namespace Dune
* Reads .vtu files and constructs a grid from the cells stored in the file
* Additionally, stored data can be read.
*
- * \tparam Grid Type of the grid to construct.
+ * NOTE: Assumption on the file structure: Each XML tag must be on a separate line.
+ *
+ * \tparam Grid The type of the grid to construct.
* \tparam GridCreator Policy type to control what to pass to a grid factory with
* data given from the file. [ContinuousGridCreator]
- *
- * Assumption on the file structure: Each XML tag must be on a separate line.
+ * \tparam FieldType Type of the components of the data to extract from the file [default: double]
**/
- template <class Grid, class GridCreator = Vtk::ContinuousGridCreator<Grid>>
+ template <class Grid, class GridCreator = Vtk::ContinuousGridCreator<Grid>, class FieldType = double>
class VtkReader
: public Vtk::FileReader<Grid, VtkReader<Grid, GridCreator>>
{
@@ -37,15 +40,25 @@ namespace Dune
POINTS, POINTS_DATA_ARRAY, CELLS, CELLS_DATA_ARRAY, APPENDED_DATA, XML_NAME, XML_NAME_ASSIGN, XML_VALUE
};
+ // Type storing information about read data
struct DataArrayAttributes
{
Vtk::DataTypes type;
- std::size_t components = 1;
+ unsigned int components = 1;
std::uint64_t offset = 0;
+ Sections section = NO_SECTION;
};
+ // Type of global world coordinates
using GlobalCoordinate = typename GridCreator::GlobalCoordinate;
+ // Template representing a grid-function that is created in getPointData() and getCellData()
+ // with Context either Vtk::PointContext or Vek::CellContext, respectively.
+ // To each GridCreate a GridFunction is associated, see, e.g. Vtk::ContinuousGridFunction
+ // or Vtk::LagrangeGridFunction.
+ template <class Context>
+ using GridFunction = typename Vtk::AssociatedGridFunction<GridCreator, FieldType, Context>::type;
+
public:
/// Constructor. Creates a new GridCreator with the passed factory
template <class... Args,
@@ -66,30 +79,52 @@ namespace Dune
/// Read the grid from file with `filename` into the GridCreator
/**
- * \param filename The name of the input file
- * \param fillCreator If `false`, only fill internal data structures, if `true`, pass the internal data to the GridCreator. [true]
+ * This function fills internal data containers representing the information from the
+ * passed file.
+ *
+ * \param filename The name of the input file
+ * \param fillCreator If `false`, only fill internal data structures, if `true`, pass
+ * the internal data to the GridCreator. [true]
**/
void read (std::string const& filename, bool fillCreator = true);
- /// Implementation of \ref FileReader interface
- static void fillFactoryImpl (GridFactory<Grid>& factory, std::string const& filename)
- {
- VtkReader reader{factory};
- reader.read(filename);
- }
-
- /// Obtains the creator of the grid
+ /// Obtains the creator of the reader
GridCreator& gridCreator ()
{
return *creator_;
}
/// Construct the actual grid using the GridCreator
+ /// NOTE: requires an aforgoing call to \ref read()
std::unique_ptr<Grid> createGrid () const
{
return creator_->createGrid();
}
+ /// Construct a grid-function representing the point-data with the given name
+ /// NOTE: requires an aforgoing call to \ref read()
+ GridFunction<Vtk::PointContext> getPointData (std::string const& name) const
+ {
+ auto const& data = dataArray_.at(name);
+ VTK_ASSERT_MSG(data.section == POINT_DATA,
+ "The data to extract is not point-data. Use `getCellData()` instead!");
+
+ return {*creator_, pointData_.at(name), data.components,
+ vec_types, vec_offsets, vec_connectivity};
+ }
+
+ /// Construct a grid-function representing the cell-data with the given name
+ /// NOTE: requires an aforgoing call to \ref read()
+ GridFunction<Vtk::CellContext> getCellData (std::string const& name) const
+ {
+ auto const& data = dataArray_.at(name);
+ VTK_ASSERT_MSG(data.section == CELL_DATA,
+ "The data to extract is not cell-data. Use `getPointData()` instead!");
+
+ return {*creator_, cellData_.at(name), data.components,
+ vec_types, vec_offsets, vec_connectivity};
+ }
+
/// Advanced read methods
/// @{
@@ -108,7 +143,7 @@ namespace Dune
void readParallelFileFromStream (std::ifstream& input, int rank, int size, bool create = true);
/// Insert all internal data to the GridCreator
- // NOTE: requires the internal data structures to be filled by an aforgoing call to read()
+ /// NOTE: requires an aforgoing call to \ref read()
void fillGridCreator (bool insertPieces = true);
/// @}
@@ -119,21 +154,28 @@ namespace Dune
return pieces_;
}
+#ifndef DOXYGEN
+ // Implementation of the FileReader interface
+ static void fillFactoryImpl (GridFactory<Grid>& factory, std::string const& filename)
+ {
+ VtkReader reader{factory};
+ reader.read(filename);
+ }
+#endif
+
private:
// Read values stored on the cells with name `name`
+ Sections readCellData (std::ifstream& /*input*/, std::string /*name*/);
+
template <class T>
- Sections readCellData (std::ifstream& /*input*/, std::vector<T>& /*values*/, std::string /*name*/)
- {
- /* does not read anything */
- return CD_DATA_ARRAY;
- }
+ void readCellDataAppended (std::ifstream& /*input*/, std::string /*name*/);
+
+ // Read values stored on the points with name `name`
+ Sections readPointData (std::ifstream& /*input*/, std::string /*name*/);
template <class T>
- Sections readPointData (std::ifstream& /*input*/, std::vector<T>& /*values*/, std::string /*name*/)
- {
- /* does not read anything */
- return PD_DATA_ARRAY;
- }
+ void readPointDataAppended (std::ifstream& /*input*/, std::string /*name*/);
+
// Read vertex coordinates from `input` stream and store in into `factory`
Sections readPoints (std::ifstream& input, std::string name);
@@ -141,11 +183,13 @@ namespace Dune
template <class T>
void readPointsAppended (std::ifstream& input);
+
// Read cell type, cell offsets and connectivity from `input` stream
Sections readCells (std::ifstream& input, std::string name);
void readCellsAppended (std::ifstream& input);
+
// Read data from appended section in vtk file, starting from `offset`
template <class T>
void readAppended (std::ifstream& input, std::vector<T>& values, std::uint64_t offset);
@@ -196,6 +240,10 @@ namespace Dune
// map Name -> {DataType,NumberOfComponents,Offset}
std::map<std::string, DataArrayAttributes> dataArray_;
+ // storage for internal point and cell data
+ std::map<std::string, std::vector<FieldType>> pointData_;
+ std::map<std::string, std::vector<FieldType>> cellData_;
+
// vector of filenames of parallel pieces
std::vector<std::string> pieces_;
diff --git a/dune/vtk/vtkreader.impl.hh b/dune/vtk/vtkreader.impl.hh
index e5507a50e4552f0d5c413a39738f89eacd8ba67e..e459d5c853f26e5d4ae2f3ec57ab514cb0a407e4 100644
--- a/dune/vtk/vtkreader.impl.hh
+++ b/dune/vtk/vtkreader.impl.hh
@@ -16,8 +16,8 @@
namespace Dune {
-template <class Grid, class Creator>
-void VtkReader<Grid,Creator>::read (std::string const& filename, bool fillCreator)
+template <class Grid, class Creator, class Field>
+void VtkReader<Grid,Creator,Field>::read (std::string const& filename, bool fillCreator)
{
// check whether file exists!
if (!Vtk::exists(filename))
@@ -38,14 +38,14 @@ void VtkReader<Grid,Creator>::read (std::string const& filename, bool fillCreato
}
-template <class Grid, class Creator>
-void VtkReader<Grid,Creator>::readSerialFileFromStream (std::ifstream& input, bool fillCreator)
+template <class Grid, class Creator, class Field>
+void VtkReader<Grid,Creator,Field>::readSerialFileFromStream (std::ifstream& input, bool fillCreator)
{
clear();
std::string compressor = "";
std::string data_name = "", data_format = "";
Vtk::DataTypes data_type = Vtk::UNKNOWN;
- std::size_t data_components = 0;
+ unsigned int data_components = 0;
std::uint64_t data_offset = 0;
Sections section = NO_SECTION;
@@ -135,7 +135,7 @@ void VtkReader<Grid,Creator>::readSerialFileFromStream (std::ifstream& input, bo
}
// Store attributes of DataArray
- dataArray_[data_name] = {data_type, data_components, data_offset};
+ dataArray_[data_name] = {data_type, data_components, data_offset, section};
// Skip section in appended mode
if (data_format == "appended") {
@@ -185,6 +185,23 @@ void VtkReader<Grid,Creator>::readSerialFileFromStream (std::ifstream& input, bo
readPointsAppended<double>(input);
readCellsAppended(input);
+
+ // read point an cell data
+ for (auto const& [name,data] : dataArray_) {
+ if (data.section == POINT_DATA) {
+ if (data.type == Vtk::FLOAT32)
+ readPointDataAppended<float>(input, name);
+ else
+ readPointDataAppended<double>(input, name);
+ }
+ else if (data.section == CELL_DATA) {
+ if (data.type == Vtk::FLOAT32)
+ readCellDataAppended<float>(input, name);
+ else
+ readCellDataAppended<double>(input, name);
+ }
+ }
+
section = NO_SECTION; // finish reading after appended section
}
else if (isSection(line, "/AppendedData", section, APPENDED_DATA))
@@ -192,25 +209,13 @@ void VtkReader<Grid,Creator>::readSerialFileFromStream (std::ifstream& input, bo
switch (section) {
case PD_DATA_ARRAY:
- if (data_type == Vtk::FLOAT32) {
- std::vector<float> values;
- section = readPointData(input, values, data_name);
- } else if (data_type == Vtk::FLOAT64) {
- std::vector<double> values;
- section = readPointData(input, values, data_name);
- }
+ section = readPointData(input, data_name);
break;
case POINTS_DATA_ARRAY:
section = readPoints(input, data_name);
break;
case CD_DATA_ARRAY:
- if (data_type == Vtk::FLOAT32) {
- std::vector<float> values;
- section = readCellData(input, values, data_name);
- } else if (data_type == Vtk::FLOAT64) {
- std::vector<double> values;
- section = readCellData(input, values, data_name);
- }
+ section = readCellData(input, data_name);
break;
case CELLS_DATA_ARRAY:
section = readCells(input, data_name);
@@ -232,8 +237,8 @@ void VtkReader<Grid,Creator>::readSerialFileFromStream (std::ifstream& input, bo
}
-template <class Grid, class Creator>
-void VtkReader<Grid,Creator>::readParallelFileFromStream (std::ifstream& input, int commRank, int commSize, bool fillCreator)
+template <class Grid, class Creator, class Field>
+void VtkReader<Grid,Creator,Field>::readParallelFileFromStream (std::ifstream& input, int commRank, int commSize, bool fillCreator)
{
clear();
@@ -275,8 +280,7 @@ void VtkReader<Grid,Creator>::readParallelFileFromStream (std::ifstream& input,
break;
}
- if (section != NO_SECTION)
- DUNE_THROW(Dune::VtkError, "VTK-File is incomplete. It must end with </VTKFile>!");
+ VTK_ASSERT_MSG(section == NO_SECTION, "VTK-File is incomplete. It must end with </VTKFile>!");
if (fillCreator)
fillGridCreator();
@@ -330,9 +334,80 @@ Sections skipRestOfDataArray (IStream& input, Sections section, Sections parent_
// @}
-template <class Grid, class Creator>
-typename VtkReader<Grid,Creator>::Sections
-VtkReader<Grid,Creator>::readPoints (std::ifstream& input, std::string name)
+// Read values stored on the cells with name `name`
+template <class Grid, class Creator, class Field>
+typename VtkReader<Grid,Creator,Field>::Sections
+VtkReader<Grid,Creator,Field>::readCellData (std::ifstream& input, std::string name)
+{
+ VTK_ASSERT(numberOfCells_ > 0);
+ unsigned int components = dataArray_[name].components;
+
+ Sections sec;
+ std::vector<Field>& values = cellData_[name];
+ sec = readDataArray(input, values, components*numberOfCells_, CD_DATA_ARRAY, CELL_DATA);
+ if (sec != CELL_DATA)
+ sec = skipRestOfDataArray(input, CD_DATA_ARRAY, CELL_DATA);
+ VTK_ASSERT(sec == CELL_DATA);
+ VTK_ASSERT(values.size() == components*numberOfCells_);
+
+ return sec;
+}
+
+
+template <class Grid, class Creator, class Field>
+ template <class T>
+void VtkReader<Grid,Creator,Field>::readCellDataAppended (std::ifstream& input, std::string name)
+{
+ VTK_ASSERT(numberOfCells_ > 0);
+ unsigned int components = dataArray_[name].components;
+
+ std::vector<T> values;
+ readAppended(input, values, dataArray_[name].offset);
+ VTK_ASSERT(values.size() == components*numberOfCells_);
+
+ cellData_[name].resize(values.size());
+ std::copy(values.begin(), values.end(), cellData_[name].begin());
+}
+
+
+template <class Grid, class Creator, class Field>
+typename VtkReader<Grid,Creator,Field>::Sections
+VtkReader<Grid,Creator,Field>::readPointData (std::ifstream& input, std::string name)
+{
+ VTK_ASSERT(numberOfPoints_ > 0);
+ unsigned int components = dataArray_[name].components;
+
+ Sections sec;
+ std::vector<Field>& values = pointData_[name];
+ sec = readDataArray(input, values, components*numberOfPoints_, PD_DATA_ARRAY, POINT_DATA);
+ if (sec != POINT_DATA)
+ sec = skipRestOfDataArray(input, PD_DATA_ARRAY, POINT_DATA);
+ VTK_ASSERT(sec == POINT_DATA);
+ VTK_ASSERT(values.size() == components*numberOfPoints_);
+
+ return sec;
+}
+
+
+template <class Grid, class Creator, class Field>
+ template <class T>
+void VtkReader<Grid,Creator,Field>::readPointDataAppended (std::ifstream& input, std::string name)
+{
+ VTK_ASSERT(numberOfPoints_ > 0);
+ unsigned int components = dataArray_[name].components;
+
+ std::vector<T> values;
+ readAppended(input, values, dataArray_[name].offset);
+ VTK_ASSERT(values.size() == components*numberOfPoints_);
+
+ pointData_[name].resize(values.size());
+ std::copy(values.begin(), values.end(), pointData_[name].begin());
+}
+
+
+template <class Grid, class Creator, class Field>
+typename VtkReader<Grid,Creator,Field>::Sections
+VtkReader<Grid,Creator,Field>::readPoints (std::ifstream& input, std::string name)
{
using T = typename GlobalCoordinate::value_type;
VTK_ASSERT(numberOfPoints_ > 0);
@@ -362,9 +437,9 @@ VtkReader<Grid,Creator>::readPoints (std::ifstream& input, std::string name)
}
-template <class Grid, class Creator>
+template <class Grid, class Creator, class Field>
template <class T>
-void VtkReader<Grid,Creator>::readPointsAppended (std::ifstream& input)
+void VtkReader<Grid,Creator,Field>::readPointsAppended (std::ifstream& input)
{
VTK_ASSERT(numberOfPoints_ > 0);
VTK_ASSERT(dataArray_["points"].components == 3u);
@@ -385,9 +460,9 @@ void VtkReader<Grid,Creator>::readPointsAppended (std::ifstream& input)
}
-template <class Grid, class Creator>
-typename VtkReader<Grid,Creator>::Sections
-VtkReader<Grid,Creator>::readCells (std::ifstream& input, std::string name)
+template <class Grid, class Creator, class Field>
+typename VtkReader<Grid,Creator,Field>::Sections
+VtkReader<Grid,Creator,Field>::readCells (std::ifstream& input, std::string name)
{
Sections sec = CELLS_DATA_ARRAY;
@@ -409,8 +484,8 @@ VtkReader<Grid,Creator>::readCells (std::ifstream& input, std::string name)
}
-template <class Grid, class Creator>
-void VtkReader<Grid,Creator>::readCellsAppended (std::ifstream& input)
+template <class Grid, class Creator, class Field>
+void VtkReader<Grid,Creator,Field>::readCellsAppended (std::ifstream& input)
{
VTK_ASSERT(numberOfCells_ > 0);
auto types_data = dataArray_["types"];
@@ -487,9 +562,9 @@ void read_compressed (T* buffer, unsigned char* buffer_in,
// @}
-template <class Grid, class Creator>
+template <class Grid, class Creator, class Field>
template <class T>
-void VtkReader<Grid,Creator>::readAppended (std::ifstream& input, std::vector<T>& values, std::uint64_t offset)
+void VtkReader<Grid,Creator,Field>::readAppended (std::ifstream& input, std::vector<T>& values, std::uint64_t offset)
{
input.seekg(offset0_ + offset);
@@ -538,8 +613,8 @@ void VtkReader<Grid,Creator>::readAppended (std::ifstream& input, std::vector<T>
}
-template <class Grid, class Creator>
-void VtkReader<Grid,Creator>::fillGridCreator (bool insertPieces)
+template <class Grid, class Creator, class Field>
+void VtkReader<Grid,Creator,Field>::fillGridCreator (bool insertPieces)
{
VTK_ASSERT(vec_points.size() == numberOfPoints_);
VTK_ASSERT(vec_types.size() == numberOfCells_);
@@ -554,8 +629,8 @@ void VtkReader<Grid,Creator>::fillGridCreator (bool insertPieces)
}
// Assume input already read the line <AppendedData ...>
-template <class Grid, class Creator>
-std::uint64_t VtkReader<Grid,Creator>::findAppendedDataPosition (std::ifstream& input) const
+template <class Grid, class Creator, class Field>
+std::uint64_t VtkReader<Grid,Creator,Field>::findAppendedDataPosition (std::ifstream& input) const
{
char c;
while (input.get(c) && std::isblank(c)) { /*do nothing*/ }
@@ -568,8 +643,8 @@ std::uint64_t VtkReader<Grid,Creator>::findAppendedDataPosition (std::ifstream&
}
-template <class Grid, class Creator>
-std::map<std::string, std::string> VtkReader<Grid,Creator>::parseXml (std::string const& line, bool& closed)
+template <class Grid, class Creator, class Field>
+std::map<std::string, std::string> VtkReader<Grid,Creator,Field>::parseXml (std::string const& line, bool& closed)
{
closed = false;
std::map<std::string, std::string> attr;
@@ -622,8 +697,8 @@ std::map<std::string, std::string> VtkReader<Grid,Creator>::parseXml (std::strin
}
-template <class Grid, class Creator>
-void VtkReader<Grid,Creator>::clear ()
+template <class Grid, class Creator, class Field>
+void VtkReader<Grid,Creator,Field>::clear ()
{
vec_points.clear();
vec_point_ids.clear();
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 0c09ee229182b20448be80877d9cb5b765aaadda..9c490c66c232d5293249e412ae963324dc26905d 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -19,6 +19,10 @@ dune_add_test(SOURCES datacollector.cc
LINK_LIBRARIES dunevtk
CMAKE_GUARD dune-functions_FOUND)
+dune_add_test(SOURCES datareader.cc
+ LINK_LIBRARIES dunevtk
+ CMAKE_GUARD dune-functions_FOUND)
+
dune_add_test(SOURCES lagrangepoints.cc
LINK_LIBRARIES dunevtk)
diff --git a/src/datareader.cc b/src/datareader.cc
new file mode 100644
index 0000000000000000000000000000000000000000..08df1196b1ddb7ba593d1499a7ee36e78cd4a261
--- /dev/null
+++ b/src/datareader.cc
@@ -0,0 +1,162 @@
+// -*- 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 <cmath>
+#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/gridfunctions/analyticgridviewfunction.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/vtk/vtkreader.hh>
+#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
+#include <dune/vtk/gridcreators/continuousgridcreator.hh>
+#include <dune/vtk/gridcreators/lagrangegridcreator.hh>
+#include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
+
+#if HAVE_DUNE_UGGRID
+#include <dune/grid/uggrid.hh>
+#endif
+
+#if HAVE_DUNE_ALUGRID
+#include <dune/alugrid/grid.hh>
+#endif
+
+using namespace Dune;
+
+template <class GridType, class GridCreator>
+void run(std::string const& prefix)
+{
+ std::tuple f_tuple{0,
+ [](auto const& x) { return std::sin(x[0]*x[0]*x[0] + 20); },
+ [](auto const& x) { return std::sin(x[0]*x[0]) + std::cos(x[1] + x[0]); },
+ [](auto const& x) { return std::sin(x[0] + x[1] + x[2]) + std::cos(x[0]*x[1]*x[2]) + std::sin(x[0]*x[2])*std::cos(x[1] + x[2]); }
+ };
+
+ constexpr int dim = GridType::dimension;
+ auto f = std::get<dim>(f_tuple);
+
+ using GridView = typename GridType::LeafGridView;
+ {
+ 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;
+
+ GridView gridView = grid.leafGridView();
+ auto data = Functions::makeAnalyticGridViewFunction(f, gridView);
+
+ VtkUnstructuredGridWriter<GridView> vtkWriter1(gridView, Vtk::ASCII);
+ vtkWriter1.addPointData(data, "f");
+ vtkWriter1.addCellData(data, "g");
+ vtkWriter1.write(prefix + "_ascii.vtu");
+
+ VtkUnstructuredGridWriter<GridView> vtkWriter2(gridView, Vtk::BINARY);
+ vtkWriter2.addPointData(data, "f");
+ vtkWriter2.addCellData(data, "g");
+ vtkWriter2.write(prefix + "_binary.vtu");
+
+ VtkUnstructuredGridWriter<GridView> vtkWriter3(gridView, Vtk::COMPRESSED, Vtk::FLOAT64);
+ vtkWriter3.addPointData(data, "f");
+ vtkWriter3.addCellData(data, "g");
+ vtkWriter3.write(prefix + "_compressed.vtu");
+
+ if constexpr(std::is_same_v<GridCreator, Vtk::LagrangeGridCreator<GridType>>) {
+ VtkUnstructuredGridWriter<GridView, Vtk::LagrangeDataCollector<GridView,3>> vtkWriter4(gridView, Vtk::BINARY, Vtk::FLOAT64);
+ vtkWriter4.addPointData(data, "f");
+ vtkWriter4.addCellData(data, "g");
+ vtkWriter4.write(prefix + "_order3.vtu");
+ }
+ }
+
+ {
+ std::cout << "read '" << prefix << "_ascii.vtu'..." << std::endl;
+ VtkReader<GridType, GridCreator, float> reader;
+ reader.read(prefix + "_ascii.vtu");
+ auto gridPtr = reader.createGrid();
+ auto& grid = *gridPtr;
+
+ auto gf = reader.getPointData("f");
+ auto gg = reader.getCellData("g");
+
+ VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::ASCII);
+ vtkWriter.addPointData(gf, "f");
+ vtkWriter.addCellData(gg, "g");
+ vtkWriter.write(prefix + "_ascii_out.vtu");
+ }
+
+ {
+ std::cout << "read '" << prefix << "_binary.vtu'..." << std::endl;
+ VtkReader<GridType, GridCreator, float> reader;
+ reader.read(prefix + "_binary.vtu");
+ auto gridPtr = reader.createGrid();
+ auto& grid = *gridPtr;
+
+ auto gf = reader.getPointData("f");
+ auto gg = reader.getCellData("g");
+
+ VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::ASCII);
+ vtkWriter.addPointData(gf, "f");
+ vtkWriter.addCellData(gg, "g");
+ vtkWriter.write(prefix + "_binary_out.vtu");
+ }
+
+ {
+ std::cout << "read '" << prefix << "_compressed.vtu'..." << std::endl;
+ VtkReader<GridType, GridCreator> reader;
+ reader.read(prefix + "_compressed.vtu");
+ auto gridPtr = reader.createGrid();
+ auto& grid = *gridPtr;
+
+ auto gf = reader.getPointData("f");
+ auto gg = reader.getCellData("g");
+
+ VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::ASCII);
+ vtkWriter.addPointData(gf, "f");
+ vtkWriter.addCellData(gg, "g");
+ vtkWriter.write(prefix + "_compressed_out.vtu");
+ }
+
+ if constexpr(std::is_same_v<GridCreator, Vtk::LagrangeGridCreator<GridType>>) {
+ std::cout << "read '" << prefix << "_order3.vtu'..." << std::endl;
+ VtkReader<GridType, GridCreator> reader;
+ reader.read(prefix + "_order3.vtu");
+ auto gridPtr = reader.createGrid();
+ auto& grid = *gridPtr;
+
+ auto gf = reader.getPointData("f");
+ auto gg = reader.getCellData("g");
+
+ VtkUnstructuredGridWriter<GridView, Vtk::LagrangeDataCollector<GridView,3>> vtkWriter(grid.leafGridView(), Vtk::ASCII);
+ vtkWriter.addPointData(gf, "f");
+ vtkWriter.addCellData(gg, "g");
+ vtkWriter.write(prefix + "_order3_out.vtu");
+ }
+}
+
+template <class GridType>
+void run(std::string const& prefix)
+{
+ run<GridType, Vtk::ContinuousGridCreator<GridType>>(prefix + "_continuous");
+ run<GridType, Vtk::LagrangeGridCreator<GridType>>(prefix + "_lagrange");
+}
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ run<UGGrid<2>>("uggrid_2d");
+ run<UGGrid<3>>("uggrid_3d");
+
+#if HAVE_DUNE_ALUGRID
+ run<Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming>>("alugrid_2d");
+ run<Dune::ALUGrid<3,3,Dune::simplex,Dune::conforming>>("alugrid_3d");
+#endif
+}