Commit f10ba622 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Added data reader to VtkReader and a grid-function representation of the read data

parent 28a8ed70
...@@ -118,20 +118,60 @@ file, these stored temporaries are collected and combined to one VTK file. Thus, ...@@ -118,20 +118,60 @@ 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 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. are stored, by default, in a `/tmp` folder, with (hopefully) fast write access.
### VtkReader ## VtkReader
Read in unstructured grid files (.vtu files) and create a new grid, using a GridFactory. Reading unstructured grid files (.vtu files) and creating a new grid, using a GridFactory,
The reader allows to create the grid in multiple ways, by providing a `GridCreator` can be performed using the `VtkReader` class. The reader allows to create the grid in
template parameter. The `ContinuousGridCreator` reads the connectivity of the grid multiple ways, by providing a `GridCreator` template parameter. The `ContinuousGridCreator`
as it is and assumes that the elements are already connected correctly. On the other reads the connectivity of the grid as it is and assumes that the elements are already
hand, a `DiscontinuousGridCreator` reconnects separated elements, by identifying connected correctly. On the other hand, a `DiscontinuousGridCreator` reconnects separated
matching coordinates of the cell vertices. elements, by identifying matching coordinates of the cell vertices. See more possible
grid-creators in the directory `dune/vtk/gridcreators`.
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
// Read the data from a file with filename
void read(std::string filename);
// 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 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 provided, all partitions can be read by either one processor and distributed later on
(`SerialGridCreator`) or read directly in parallel (`ParallelGridCreator`). The later (`SerialGridCreator`) or read directly in parallel (`ParallelGridCreator`). The later
is currently only available in dune-alugrid 2.6. is currently only available in dune-alugrid 2.6.
**TODO:**
- Provide an interface to read the points-/cell-data from the file
- Extent the implementation to other file formats
...@@ -27,6 +27,7 @@ install(FILES ...@@ -27,6 +27,7 @@ install(FILES
add_subdirectory(datacollectors) add_subdirectory(datacollectors)
add_subdirectory(gridcreators) add_subdirectory(gridcreators)
add_subdirectory(gridfunctions)
add_subdirectory(test) add_subdirectory(test)
add_subdirectory(utility) add_subdirectory(utility)
add_subdirectory(writers) add_subdirectory(writers)
...@@ -67,8 +67,17 @@ namespace Dune ...@@ -67,8 +67,17 @@ namespace Dune
return comp < N ? vec[comp] : 0.0; return comp < N ? vec[comp] : 0.0;
} }
// 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 // Return the scalar values
template <class T> template <class T,
std::enable_if_t<not IsIndexable<T,int>::value, int> = 0>
double evaluateImpl (int comp, T const& value) const double evaluateImpl (int comp, T const& value) const
{ {
assert(comp == 0); assert(comp == 0);
......
...@@ -29,9 +29,21 @@ namespace Dune ...@@ -29,9 +29,21 @@ namespace Dune
using Derived = DerivedType; using Derived = DerivedType;
public: public:
/// Constructor. Stores a reference to the passed GridFactory /// Constructor. Stores a reference to the passed GridFactory.
GridCreatorInterface (GridFactory<Grid>& factory) 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 /// Insert all points as vertices into the factory
...@@ -112,7 +124,7 @@ namespace Dune ...@@ -112,7 +124,7 @@ namespace Dune
} }
protected: protected:
GridFactory<Grid>* factory_; std::shared_ptr<GridFactory<Grid>> factory_;
}; };
} // end namespace Vtk } // end namespace Vtk
......
...@@ -11,6 +11,7 @@ ...@@ -11,6 +11,7 @@
#include <dune/vtk/types.hh> #include <dune/vtk/types.hh>
#include <dune/vtk/gridcreatorinterface.hh> #include <dune/vtk/gridcreatorinterface.hh>
#include <dune/vtk/gridfunctions/continuousgridfunction.hh>
namespace Dune namespace Dune
{ {
...@@ -26,10 +27,7 @@ namespace Dune ...@@ -26,10 +27,7 @@ namespace Dune
using Super = GridCreatorInterface<Grid, ContinuousGridCreator>; using Super = GridCreatorInterface<Grid, ContinuousGridCreator>;
using GlobalCoordinate = typename Super::GlobalCoordinate; using GlobalCoordinate = typename Super::GlobalCoordinate;
ContinuousGridCreator (GridFactory<Grid>& factory) using Super::Super;
: Super(factory)
{}
using Super::factory; using Super::factory;
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points, void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
...@@ -69,5 +67,16 @@ namespace Dune ...@@ -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 Vtk
} // end namespace Dune } // end namespace Dune
...@@ -22,9 +22,11 @@ namespace Dune ...@@ -22,9 +22,11 @@ namespace Dune
using Grid = typename GridCreator::Grid; using Grid = typename GridCreator::Grid;
using GlobalCoordinate = typename Super::GlobalCoordinate; using GlobalCoordinate = typename Super::GlobalCoordinate;
DerivedGridCreator (GridFactory<Grid>& factory) template <class... Args,
: Super(factory) disableCopyMove<DerivedGridCreator, Args...> = 0>
, gridCreator_(factory) DerivedGridCreator (Args&&... args)
: Super(std::forward<Args>(args)...)
, gridCreator_(Super::factory())
{} {}
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points, void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
......
...@@ -40,11 +40,9 @@ namespace Dune ...@@ -40,11 +40,9 @@ namespace Dune
} }
}; };
DiscontinuousGridCreator (GridFactory<Grid>& factory) using Super::Super;
: Super(factory)
{}
using Super::factory; using Super::factory;
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points, void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
std::vector<std::uint64_t> const& /*point_ids*/) std::vector<std::uint64_t> const& /*point_ids*/)
{ {
...@@ -97,5 +95,10 @@ namespace Dune ...@@ -97,5 +95,10 @@ namespace Dune
std::map<GlobalCoordinate, std::size_t, CoordLess> uniquePoints_; std::map<GlobalCoordinate, std::size_t, CoordLess> uniquePoints_;
}; };
// deduction guides
template <class Grid>
DiscontinuousGridCreator(GridFactory<Grid>&)
-> DiscontinuousGridCreator<Grid>;
} // end namespace Vtk } // end namespace Vtk
} // end namespace Dune } // end namespace Dune
...@@ -15,6 +15,7 @@ ...@@ -15,6 +15,7 @@
#include <dune/vtk/types.hh> #include <dune/vtk/types.hh>
#include <dune/vtk/gridcreatorinterface.hh> #include <dune/vtk/gridcreatorinterface.hh>
#include <dune/vtk/gridfunctions/lagrangegridfunction.hh>
#include <dune/vtk/utility/lagrangepoints.hh> #include <dune/vtk/utility/lagrangepoints.hh>
namespace Dune namespace Dune
...@@ -58,11 +59,11 @@ namespace Dune ...@@ -58,11 +59,11 @@ namespace Dune
class LocalFunction; class LocalFunction;
public: public:
using Super::factory; using LocalGeometry = MultiLinearGeometry<typename Element::Geometry::ctype,Element::dimension,Element::dimension>;
LagrangeGridCreator (GridFactory<GridType>& factory) public:
: Super(factory) using Super::Super;
{} using Super::factory;
/// Implementation of the interface function `insertVertices()` /// Implementation of the interface function `insertVertices()`
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points, void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
...@@ -163,22 +164,41 @@ namespace Dune ...@@ -163,22 +164,41 @@ namespace Dune
* The returned LocalParametrization is a mapping `GlobalCoordinate(LocalCoordinate)` * The returned LocalParametrization is a mapping `GlobalCoordinate(LocalCoordinate)`
* where `LocalCoordinate is w.r.t. the local coordinate system in the passed element * where `LocalCoordinate is w.r.t. the local coordinate system in the passed element
* and `GlobalCoordinate` a world coordinate in the parametrized grid. * 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 LocalParametrization localParametrization (Element const& element) const
{ {
assert(!nodes_.empty() && !parametrization_.empty()); VTK_ASSERT(!nodes_.empty() && !parametrization_.empty());
unsigned int insertionIndex = factory().insertionIndex(element); unsigned int insertionIndex = factory().insertionIndex(element);
auto const& localParam = parametrization_.at(insertionIndex); 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 // collect indices of vertices
std::vector<unsigned int> indices(element.subEntities(GridType::dimension)); std::vector<unsigned int> indices(element.subEntities(GridType::dimension));
for (unsigned int i = 0; i < element.subEntities(GridType::dimension); ++i) for (unsigned int i = 0; i < element.subEntities(GridType::dimension); ++i)
...@@ -188,31 +208,41 @@ namespace Dune ...@@ -188,31 +208,41 @@ namespace Dune
std::vector<unsigned int> permutation(indices.size()); std::vector<unsigned int> permutation(indices.size());
for (std::size_t i = 0; i < indices.size(); ++i) { for (std::size_t i = 0; i < indices.size(); ++i) {
auto it = std::find(localParam.corners.begin(), localParam.corners.end(), indices[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); 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 /// Determine lagrange order from number of points
template <class LocalParam> int order (GeometryType type, std::size_t nNodes) const
int order (LocalParam const& localParam) const
{ {
GeometryType type = localParam.type; for (int o = 1; o <= int(nNodes); ++o)
int nNodes = localParam.nodes.size();
for (int o = 1; o <= nNodes; ++o)
if (numLagrangePoints(type.id(), type.dim(), o) == std::size_t(nNodes)) if (numLagrangePoints(type.id(), type.dim(), o) == std::size_t(nNodes))
return o; return o;
return 1; 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 /// Determine lagrange order from number of points from the first element parametrization
int order () const int order () const
{ {
assert(!parametrization_.empty()); assert(!parametrization_.empty());
return order(parametrization_.front()); auto const& localParam = parametrization_.front();
return order(localParam);
} }
public: public:
...@@ -274,6 +304,16 @@ namespace Dune ...@@ -274,6 +304,16 @@ namespace Dune
Parametrization parametrization_; 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> template <class Grid>
class LagrangeGridCreator<Grid>::LocalParametrization class LagrangeGridCreator<Grid>::LocalParametrization
...@@ -300,16 +340,11 @@ namespace Dune ...@@ -300,16 +340,11 @@ namespace Dune
} }
/// Construct a local element parametrization for elements with permuted corners /// Construct a local element parametrization for elements with permuted corners
template <class Nodes, class LocalParam, class Permutation> template <class Nodes, class LocalParam, class LG>
LocalParametrization (Nodes const& nodes, LocalParam const& param, int order, Permutation const& permutation) LocalParametrization (Nodes const& nodes, LocalParam const& param, int order, LG&& localGeometry)
: LocalParametrization(nodes, param, order) : LocalParametrization(nodes, param, order)
{ {
auto refElem = referenceElement<ctype,Grid::dimension>(param.type); localGeometry_.emplace(std::forward<LG>(localGeometry));
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);
} }
/// Evaluate the local parametrization in local coordinates /// Evaluate the local parametrization in local coordinates
......
...@@ -27,9 +27,9 @@ namespace Dune ...@@ -27,9 +27,9 @@ namespace Dune
// The GridFactory must support insertion of global vertex IDs // The GridFactory must support insertion of global vertex IDs
static_assert(Std::is_detected<HasInsertVertex, GridFactory<Grid>, GlobalCoordinate, VertexId>{}, ""); static_assert(Std::is_detected<HasInsertVertex, GridFactory<Grid>, GlobalCoordinate, VertexId>{}, "");
ParallelGridCreator (GridFactory<Grid>& factory) public:
: Super(factory)
{} using Super::Super;
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points, void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
std::vector<std::uint64_t> const& point_ids) std::vector<std::uint64_t> const& point_ids)
...@@ -48,5 +48,10 @@ namespace Dune ...@@ -48,5 +48,10 @@ namespace Dune
} }
}; };
// deduction guides
template <class Grid>
ParallelGridCreator(GridFactory<Grid>&)
-> ParallelGridCreator<Grid>;
} // end namespace Vtk } // end namespace Vtk
} // end namespace Dune } // end namespace Dune
...@@ -21,9 +21,9 @@ namespace Dune ...@@ -21,9 +21,9 @@ namespace Dune
using Super = GridCreatorInterface<Grid, Self>; using Super = GridCreatorInterface<Grid, Self>;
using GlobalCoordinate = typename Super::GlobalCoordinate; using GlobalCoordinate = typename Super::GlobalCoordinate;
SerialGridCreator (GridFactory<Grid>& factory) public:
: Super(factory)
{} using Super::Super;
void insertVerticesImpl (std::vector<GlobalCoordinate> const& points, void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
std::vector<std::uint64_t> const& /*point_ids*/) std::vector<std::uint64_t> const& /*point_ids*/)
...@@ -72,5 +72,10 @@ namespace Dune ...@@ -72,5 +72,10 @@ namespace Dune
std::vector<std::int64_t> shift_; std::vector<std::int64_t> shift_;
}; };
// deduction guides
template <class Grid>
SerialGridCreator(GridFactory<Grid>&)
-> SerialGridCreator<Grid>;
} // end namespace Vtk } // end namespace Vtk
} // end namespace Dune } // end namespace Dune
#install headers
install(FILES
common.hh
continuousgridfunction.hh
lagrangegridfunction.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vtk/gridfunctions)
#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
#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