diff --git a/dune/vtk/filereader.hh b/dune/vtk/filereader.hh
index 10a66d98c6c5b5e18f554cfc77521a6dc3b1e767..b473ec5e7c4a3384ba34c4426215d607a1f5a655 100644
--- a/dune/vtk/filereader.hh
+++ b/dune/vtk/filereader.hh
@@ -4,23 +4,67 @@
 
 namespace Dune
 {
-  template <class Grid>
+  template <class Grid, class GridReaderImp>
   class FileReader
   {
+  private:
+    // type of underlying implementation, for internal use only
+    using Implementation = GridReaderImp;
+
+    //! \brief An accessor class to call protected members of reader implementations.
+    struct Accessor : public Implementation
+    {
+      template <class... Args>
+      static std::unique_ptr<Grid> readImp (Args&&... args)
+      {
+        return Implementation::readImpl(std::forward<Args>(args)...);
+      }
+
+      template <class... Args>
+      static void readFactoryImpl (Args&&... args)
+      {
+        return Implementation::readFactoryImpl(std::forward<Args>(args)...);
+      }
+    };
+
   public:
-    /// Virtual destructor
-    virtual ~FileReader() = default;
+    //! Reads the grid from a file with filename and returns a unique_ptr to the created grid.
+    //! Redirects to concrete implementation of derivated class.
+    template <class... Args>
+    static std::unique_ptr<Grid> read (const std::string &filename, Args&&... args)
+    {
+      return Accessor::readImpl(filename, std::forward<Args>(args)...);
+    }
+
+    //! Reads the grid from a file with filename into a grid-factory.
+    //! Redirects to concrete implementation of derivated class.
+    template <class... Args>
+    static void read (GridFactory<Grid> &factory, const std::string &filename, Args&&... args)
+    {
+      Accessor::readFactoryImpl(factory, filename, std::forward<Args>(args)...);
+    }
 
-    /// Reads the grid from a file with `filename` and returns a unique_ptr to the created grid.
-    virtual std::unique_ptr<Grid> read(std::string const& filename)
+  protected: // default implementations
+
+    // Default implementation, redirects to factory read implementation.
+    template <class... Args>
+    static std::unique_ptr<Grid> readImpl (const std::string &filename, Args&&... args)
     {
       GridFactory<Grid> factory;
-      read(factory, filename);
+      read(factory, filename, std::forward<Args>(args)...);
+
       return std::unique_ptr<Grid>{ factory.createGrid() };
     }
 
-    /// Reads the grid from a file with `filename` into a grid-factory.
-    virtual void read(GridFactory<Grid>& factory, std::string const& filename) = 0;
+    // Default implementation for reading into grid-factory: produces a runtime-error.
+    template <class... Args>
+    static void readFactoryImpl (GridFactory<Grid> &/*factory*/, const std::string &/*filename*/,
+                                 Args&&... /*args*/)
+    {
+      DUNE_THROW(NotImplemented,
+        "GridReader using a factory argument not implemented for concrete reader implementation.");
+    }
   };
 
+
 } // end namespace Dune
diff --git a/dune/vtk/vtkfunction.hh b/dune/vtk/vtkfunction.hh
index 827629f326219f9c44ac754b247f4347b5daa6f9..1c6820b5ee4a615c0bc5a9e31cb6acfc2b3a5ed5 100644
--- a/dune/vtk/vtkfunction.hh
+++ b/dune/vtk/vtkfunction.hh
@@ -7,6 +7,7 @@
 
 namespace Dune
 {
+  /// An abstract base class for LocalFunctions
   template <class GridView>
   class VTKLocalFunctionInterface
   {
@@ -14,16 +15,16 @@ namespace Dune
     using Entity = typename GridView::template Codim<0>::Entity;
     using LocalCoordinate = typename Entity::Geometry::LocalCoordinate;
 
-    //! Bind the function to the grid entity
+    /// Bind the function to the grid entity
     virtual void bind (Entity const& entity) = 0;
 
-    //! Unbind from the currently bound entity
+    /// Unbind from the currently bound entity
     virtual void unbind () = 0;
 
-    //! Evaluate single component comp in the entity at local coordinates xi
+    /// Evaluate single component comp in the entity at local coordinates xi
     virtual double evaluate (int comp, LocalCoordinate const& xi) const = 0;
 
-    //! Virtual destructor
+    /// Virtual destructor
     virtual ~VTKLocalFunctionInterface () = default;
   };
 
@@ -64,14 +65,17 @@ namespace Dune
       }
 
     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 y[comp];
       }
 
-      double evaluateImpl (int /*comp*/, LocalCoordinate const& xi, std::false_type) const
+      // Return the scalar values
+      double evaluateImpl (int comp, LocalCoordinate const& xi, std::false_type) const
       {
+        assert(comp == 0);
         return this->get()(xi);
       }
     };
@@ -97,18 +101,19 @@ namespace Dune
 
     VTKLocalFunction () = default;
 
-    //! Bind the function to the grid entity
+    /// Bind the function to the grid entity
     void bind (Entity const& entity)
     {
       this->asInterface().bind(entity);
     }
 
-    //! Unbind from the currently bound 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);
@@ -117,14 +122,15 @@ namespace Dune
 
   // ---------------------------------------------------------------------------
 
+  /// An abstract base class for GlobalFunctions
   template <class GridView>
   class VTKFunctionInterface
   {
   public:
-    //! Create a local function
+    /// Create a local function
     virtual VTKLocalFunction<GridView> makelocalFunction () const = 0;
 
-    //! Virtual destructor
+    /// Virtual destructor
     virtual ~VTKFunctionInterface () = default;
   };
 
@@ -156,27 +162,35 @@ namespace Dune
     using Super = Functions::TypeErasureBase<VTKFunctionInterface<GridView>,
                                              VTKFunctionImpl<GridView>::template Model>;
 
+    template <class F>
+    using IsGridFunction = decltype(localFunction(std::declval<F>()));
+
   public:
     template <class F, disableCopyMove<VTKFunction, F> = 0>
     VTKFunction (F&& f, std::string name, int ncomps = 1)
       : Super(std::forward<F>(f))
       , name_(std::move(name))
       , ncomps_(ncomps)
-    {}
+    {
+      static_assert(Std::is_detected<IsGridFunction,F>::value,
+        "Requires A GridFunction to be passed to the VTKFunction.");
+    }
 
     VTKFunction () = default;
 
-    //! Bind the function to the grid entity
+    /// Create a LocalFunction
     friend VTKLocalFunction<GridView> localFunction (VTKFunction const& self)
     {
       return self.asInterface().makelocalFunction();
     }
 
+    /// Return the number of components of the Range
     int ncomps () const
     {
-      return ncomps_;
+      return ncomps_; // TODO: detect automatically. Is this always possible?
     }
 
+    /// Return a name associated with the function
     std::string name () const
     {
       return name_;
diff --git a/dune/vtk/vtkreader.hh b/dune/vtk/vtkreader.hh
index dacbc42f9e6e7a7de6cbe45297b970da1bb2e9eb..6e327c24007a0a8eb200cfad2789ceacf5f2023f 100644
--- a/dune/vtk/vtkreader.hh
+++ b/dune/vtk/vtkreader.hh
@@ -9,10 +9,17 @@
 namespace Dune
 {
   /// File-Reader for Vtk .vtu files
+  /**
+   * Reads .vtu files and constructs a grid from the cells stored in the file
+   * Additionally, stored data can be read.
+   *
+   * Assumption on the file structure: Each XML tag must be on a separate line.
+   **/
   template <class Grid>
   class VtkReader
-      : public FileReader<Grid>
+      : public FileReader<Grid, VtkReader<Grid>>
   {
+    // Sections visited during the xml parsing
     enum Sections {
       NO_SECTION = 0, VTK_FILE, UNSTRUCTURED_GRID, PIECE, POINT_DATA, PD_DATA_ARRAY, CELL_DATA, CD_DATA_ARRAY,
       POINTS, POINTS_DATA_ARRAY, CELLS, CELLS_DATA_ARRAY, APPENDED_DATA, XML_NAME, XML_NAME_ASSIGN, XML_VALUE
@@ -22,42 +29,51 @@ namespace Dune
     using GlobalCoordinate = typename Entity::Geometry::GlobalCoordinate;
 
   public:
-    /// Constructor
-    VtkReader ()
-      : buffer_(block_size)
+    /// Constructor. Stores a pointer to the GridFactory.
+    VtkReader (GridFactory<Grid>& factory)
+      : factory_(&factory)
     {}
 
-    /// read the file
-    virtual void read (GridFactory<Grid>& factory, std::string const& filename) override;
-    using FileReader<Grid>::read;
+    /// Read the grid from file with `filename` into the GridFactory `factory`
+    void readFromFile (std::string const& filename);
+
+    /// Implementation of \ref FileReader interface
+    static void readFactoryImpl (GridFactory<Grid>& factory, std::string const& filename)
+    {
+      VtkReader reader{factory};
+      reader.readFromFile(filename);
+    }
 
   private:
-    Sections readCellData (std::ifstream&,
-                           GridFactory<Grid>& /*factory*/,
+    // Read values stored on the cells with name `name`
+    template <class T>
+    Sections readCellData (std::ifstream& /*input*/,
+                           std::vector<T>& /*values*/,
                            std::string /*name*/,
                            Vtk::DataTypes /*type*/,
                            std::size_t /*nComponents*/,
                            std::string /*format*/,
-                           std::size_t /*offset*/)
+                           std::uint64_t /*offset*/)
     {
       /* does not read anything */
       return CD_DATA_ARRAY;
     }
 
-    Sections readPointData (std::ifstream&,
-                            GridFactory<Grid>& /*factory*/,
+    template <class T>
+    Sections readPointData (std::ifstream& /*input*/,
+                            std::vector<T>& /*values*/,
                             std::string /*name*/,
                             Vtk::DataTypes /*type*/,
                             std::size_t /*nComponents*/,
                             std::string /*format*/,
-                            std::size_t /*offset*/)
+                            std::uint64_t /*offset*/)
     {
       /* does not read anything */
       return PD_DATA_ARRAY;
     }
 
+    // Read vertex coordinates from `input` stream and store in into `factory`
     Sections readPoints (std::ifstream& input,
-                         GridFactory<Grid>& factory,
                          std::string name,
                          Vtk::DataTypes type,
                          std::size_t nComponents,
@@ -65,11 +81,10 @@ namespace Dune
                          std::uint64_t offset);
 
     template <class T>
-    void readPointsAppended (std::ifstream& input,
-                             GridFactory<Grid>& factory);
+    void readPointsAppended (std::ifstream& input);
 
+    // Read cell type, cell offsets and connectivity from `input` stream
     Sections readCells (std::ifstream& input,
-                        GridFactory<Grid>& factory,
                         std::string name,
                         Vtk::DataTypes type,
                         std::string format,
@@ -77,13 +92,15 @@ namespace Dune
 
     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);
 
-    inline bool isSection (std::string line,
-                           std::string key,
-                           Sections current,
-                           Sections parent = NO_SECTION)
+    // Test whether line belongs to section
+    bool isSection (std::string line,
+                    std::string key,
+                    Sections current,
+                    Sections parent = NO_SECTION)
     {
       bool result = line.substr(1, key.length()) == key;
       if (result && current != parent)
@@ -91,27 +108,33 @@ namespace Dune
       return result;
     }
 
-    std::map<std::string, std::string> parseXml(std::string const& line);
+    // Read attributes from current xml tag
+    std::map<std::string, std::string> parseXml(std::string const& line, bool& closed);
 
-    void createGrid(GridFactory<Grid>& factory) const;
+    // Construct a grid using the GridFactory `factory` and the read vectors
+    // \ref vec_types, \ref vec_offsets, and \ref vec_connectivity
+    void createGrid() const;
 
   private:
+    GridFactory<Grid>* factory_;
+
+    /// Data format, i.e. ASCII, BINARY or COMPRESSED. Read from xml attributes.
     Vtk::FormatTypes format_;
 
-    std::vector<std::uint8_t> vec_types;
-    std::vector<std::int64_t> vec_offsets;
-    std::vector<std::int64_t> vec_connectivity;
+    // Temporary data to construct the grid elements
+    std::vector<std::uint8_t> vec_types; //< VTK cell type ID
+    std::vector<std::int64_t> vec_offsets; //< offset of vertices of cell
+    std::vector<std::int64_t> vec_connectivity; //< vertex indices of cell
 
-    std::size_t numCells_;
-    std::size_t numVertices_;
-    std::size_t numData_;
+    std::size_t numCells_; //< Number of cells in the grid
+    std::size_t numVertices_; // Number of vertices in the grid
 
+    // offset information for appended data
     // map Name -> {DataType,Offset}
     std::map<std::string, std::pair<Vtk::DataTypes,std::uint64_t>> offsets_;
-    std::uint64_t offset0_;
 
-    std::size_t const block_size = 1024*32;
-    std::vector<unsigned char> buffer_;
+    /// Offset of beginning of appended data
+    std::uint64_t offset0_;
   };
 
 } // end namespace Dune
diff --git a/dune/vtk/vtkreader.impl.hh b/dune/vtk/vtkreader.impl.hh
index 4739e91f72d3b50d4b64ffd04fb2db3fe3d4b03d..6dcba2daa3b3841d7f751ea43ac0b1d90689ef93 100644
--- a/dune/vtk/vtkreader.impl.hh
+++ b/dune/vtk/vtkreader.impl.hh
@@ -2,19 +2,20 @@
 #include <fstream>
 #include <iterator>
 #include <string>
-//#include <regex>
 
 #ifdef HAVE_ZLIB
 #include <zlib.h>
 #endif
 
+#include <dune/common/classname.hh>
+
 #include "utility/filesystem.hh"
 #include "utility/string.hh"
 
 namespace Dune {
 
 template <class Grid>
-void VtkReader<Grid>::read(GridFactory<Grid>& factory, std::string const& filename)
+void VtkReader<Grid>::readFromFile (std::string const& filename)
 {
   // check whether file exists!
   if (!filesystem::exists(filename))
@@ -26,7 +27,7 @@ void VtkReader<Grid>::read(GridFactory<Grid>& factory, std::string const& filena
               data_format = "";
   Vtk::DataTypes data_type = Vtk::UNKNOWN;
   std::size_t data_components = 0;
-  std::size_t data_offset = 0;
+  std::uint64_t data_offset = 0;
 
   std::string file_type = "",
               byte_order = "",
@@ -40,15 +41,17 @@ void VtkReader<Grid>::read(GridFactory<Grid>& factory, std::string const& filena
     ltrim(line);
 
     if (isSection(line, "VTKFile", section)) {
-
-      auto attr = parseXml(line);
+      bool closed = false;
+      auto attr = parseXml(line, closed);
       file_type = attr["type"];
       if (file_type != "UnstructuredGrid")
         DUNE_THROW(NotImplemented, "Only UnstructuredGrid format implemented. Found: " << file_type);
       if (!attr["version"].empty())
-        version = std::stod(attr["version"]);
-      byte_order = attr["byte_order"];
-      header_type = attr["header_type"];
+        assert(std::stod(attr["version"]) == 1.0);
+      if (!attr["header_type"].empty())
+        assert(attr["byte_order"] == "LittleEndian");
+      if (!attr["header_type"].empty())
+        assert(attr["header_type"] == "UInt64");
       if (!attr["compressor"].empty())
         compressor = attr["compressor"];
       section = VTK_FILE;
@@ -60,7 +63,8 @@ void VtkReader<Grid>::read(GridFactory<Grid>& factory, std::string const& filena
     else if (isSection(line, "/UnstructuredGrid", section, UNSTRUCTURED_GRID))
       section = VTK_FILE;
     else if (isSection(line, "Piece", section, UNSTRUCTURED_GRID)) {
-      auto attr = parseXml(line);
+      bool closed = false;
+      auto attr = parseXml(line, closed);
       numVertices_ = std::stol(attr["NumberOfPoints"]);
       numCells_ = std::stol(attr["NumberOfCells"]);
       section = PIECE;
@@ -84,7 +88,8 @@ void VtkReader<Grid>::read(GridFactory<Grid>& factory, std::string const& filena
     else if (isSection(line, "/Cells", section, CELLS))
       section = PIECE;
     else if (line.substr(1,9) == "DataArray") {
-      auto attr = parseXml(line);
+      bool closed = false;
+      auto attr = parseXml(line, closed);
 
       data_type = Vtk::Map::datatype[attr["type"]];
       data_name = attr["Name"];
@@ -107,7 +112,8 @@ void VtkReader<Grid>::read(GridFactory<Grid>& factory, std::string const& filena
         assert(data_format == "appended");
       }
 
-      if (section == POINT_DATA)
+      if (closed) continue;
+      else if (section == POINT_DATA)
         section = PD_DATA_ARRAY;
       else if (section == POINTS)
         section = POINTS_DATA_ARRAY;
@@ -131,37 +137,52 @@ void VtkReader<Grid>::read(GridFactory<Grid>& factory, std::string const& filena
         DUNE_THROW(Exception, "Wrong section for </DataArray>");
     }
     else if (isSection(line, "AppendedData", section, VTK_FILE)) {
-      auto attr = parseXml(line);
+      bool closed = false;
+      auto attr = parseXml(line, closed);
       encoding = attr["encoding"];
       if (encoding != "raw")
         DUNE_THROW(NotImplemented, "Binary encoding != raw not implemented.");
-      offset0_ = input.tellg() + 1;
-      section = APPENDED_DATA;
+
+      // Find starting point of appended data
+      while (input.get() != '_')
+        /*do nothing*/;
+
+      offset0_ = input.tellg()+1;
+      if (offsets_["points"].first == Vtk::FLOAT32)
+        readPointsAppended<float>(input);
+      else
+        readPointsAppended<double>(input);
+
+      readCellsAppended(input);
+      section = NO_SECTION; // finish reading after appended section
     }
     else if (isSection(line, "/AppendedData", section, APPENDED_DATA))
       section = VTK_FILE;
 
     switch (section) {
       case PD_DATA_ARRAY:
-        section = readPointData(input, factory, data_name, data_type, data_components, data_format, data_offset);
+        if (data_type == Vtk::FLOAT32) {
+          std::vector<float> values;
+          section = readPointData(input, values, data_name, data_type, data_components, data_format, data_offset);
+        } else {
+          std::vector<double> values;
+          section = readPointData(input, values, data_name, data_type, data_components, data_format, data_offset);
+        }
         break;
       case POINTS_DATA_ARRAY:
-        section = readPoints(input, factory, data_name, data_type, data_components, data_format, data_offset);
+        section = readPoints(input, data_name, data_type, data_components, data_format, data_offset);
         break;
       case CD_DATA_ARRAY:
-        section = readCellData(input, factory, data_name, data_type, data_components, data_format, data_offset);
+        if (data_type == Vtk::FLOAT32) {
+          std::vector<float> values;
+          section = readCellData(input, values, data_name, data_type, data_components, data_format, data_offset);
+        } else {
+          std::vector<double> values;
+          section = readCellData(input, values, data_name, data_type, data_components, data_format, data_offset);
+        }
         break;
       case CELLS_DATA_ARRAY:
-        section = readCells(input, factory, data_name, data_type, data_format, data_offset);
-        break;
-      case APPENDED_DATA:
-        if (offsets_["points"].first == Vtk::FLOAT32)
-          readPointsAppended<float>(input, factory);
-        else
-          readPointsAppended<double>(input, factory);
-
-        readCellsAppended(input);
-        section = NO_SECTION; // finish reading after appended section
+        section = readCells(input, data_name, data_type, data_format, data_offset);
         break;
       default:
         // do nothing
@@ -175,16 +196,23 @@ void VtkReader<Grid>::read(GridFactory<Grid>& factory, std::string const& filena
   if (section != NO_SECTION)
     DUNE_THROW(IOError, "VTK-File is incomplete. It must end with </VTKFile>!");
 
-  createGrid(factory);
+  createGrid();
 }
 
 
+// @{ implementation detail
+/**
+ * Read ASCII data from `input` stream into vector `values`
+ * \param max_size  Upper bound for the number of values
+ * \param section   Current XML section you are reading in
+ * \param parent_section   XML Section to return when current `section` is finished.
+ **/
 template <class IStream, class T, class Sections>
-Sections read_data_array(IStream& input, std::vector<T>& values, std::size_t max_size,
-                         Sections section, Sections parent_section)
+Sections read_data_array (IStream& input, std::vector<T>& values, std::size_t max_size,
+                          Sections section, Sections parent_section)
 {
   values.reserve(max_size);
-  using S = std::conditional_t<sizeof(T) <= 8, std::uint16_t, T>;
+  using S = std::conditional_t<(sizeof(T) <= 1), std::uint16_t, T>;
 
   std::size_t idx = 0;
   std::string line;
@@ -201,13 +229,13 @@ Sections read_data_array(IStream& input, std::vector<T>& values, std::size_t max
 
   return section;
 }
+// @}
 
 
 template <class Grid>
 typename VtkReader<Grid>::Sections
-VtkReader<Grid>::readPoints(std::ifstream& input, GridFactory<Grid>& factory,
-                            std::string name, Vtk::DataTypes type,
-                            std::size_t nComponents, std::string format, std::uint64_t offset)
+VtkReader<Grid>::readPoints (std::ifstream& input, std::string name, Vtk::DataTypes type,
+                             std::size_t nComponents, std::string format, std::uint64_t offset)
 {
   if (format == "appended") {
     offsets_["points"] = {type, offset};
@@ -228,7 +256,7 @@ VtkReader<Grid>::readPoints(std::ifstream& input, GridFactory<Grid>& factory,
     for (std::size_t j = 0; j < p.size(); ++j)
       p[j] = point_values[idx++];
     idx += (3u - p.size());
-    factory.insertVertex(p);
+    factory_->insertVertex(p);
   }
 
   return sec;
@@ -237,7 +265,7 @@ VtkReader<Grid>::readPoints(std::ifstream& input, GridFactory<Grid>& factory,
 
 template <class Grid>
   template <class T>
-void VtkReader<Grid>::readPointsAppended(std::ifstream& input, GridFactory<Grid>& factory)
+void VtkReader<Grid>::readPointsAppended (std::ifstream& input)
 {
   assert(numVertices_ > 0);
   auto offset_data = offsets_["points"];
@@ -253,15 +281,15 @@ void VtkReader<Grid>::readPointsAppended(std::ifstream& input, GridFactory<Grid>
       p[j] = T(point_values[idx++]);
     idx += (3u - p.size());
 
-    factory.insertVertex(p);
+    factory_->insertVertex(p);
   }
 }
 
 
 template <class Grid>
 typename VtkReader<Grid>::Sections
-VtkReader<Grid>::readCells(std::ifstream& input, GridFactory<Grid>& factory,
-                           std::string name, Vtk::DataTypes type, std::string format, std::uint64_t offset)
+VtkReader<Grid>::readCells (std::ifstream& input, std::string name, Vtk::DataTypes type,
+                            std::string format, std::uint64_t offset)
 {
   if (format == "appended") {
     offsets_[name] = {type, offset};
@@ -292,7 +320,7 @@ VtkReader<Grid>::readCells(std::ifstream& input, GridFactory<Grid>& factory,
 
 
 template <class Grid>
-void VtkReader<Grid>::readCellsAppended(std::ifstream& input)
+void VtkReader<Grid>::readCellsAppended (std::ifstream& input)
 {
   assert(numCells_ > 0);
   auto types_data = offsets_["types"];
@@ -313,6 +341,14 @@ void VtkReader<Grid>::readCellsAppended(std::ifstream& input)
 }
 
 
+// @{ implementation detail
+/**
+ * Read compressed data into `buffer_in`, uncompress it and store the result in
+ * the concrete-data-type `buffer`
+ * \param bs     Size of the uncompressed data
+ * \param cbs    Size of the compressed data
+ * \param input  Stream to read from.
+ **/
 template <class T, class IStream>
 void read_compressed (T* buffer, unsigned char* buffer_in,
                       std::uint64_t bs, std::uint64_t cbs, IStream& input)
@@ -324,7 +360,8 @@ void read_compressed (T* buffer, unsigned char* buffer_in,
   Bytef* compressed_buffer = reinterpret_cast<Bytef*>(buffer_in);
   Bytef* uncompressed_buffer = reinterpret_cast<Bytef*>(buffer);
 
-  input.read((char*)(compressed_buffer), compressed_space);
+  input.readsome((char*)(compressed_buffer), compressed_space);
+  assert(input.gcount() == compressed_space);
 
   if (uncompress(uncompressed_buffer, &uncompressed_space, compressed_buffer, compressed_space) != Z_OK) {
     std::cerr << "Zlib error while uncompressing data.\n";
@@ -335,6 +372,7 @@ void read_compressed (T* buffer, unsigned char* buffer_in,
   std::abort();
 #endif
 }
+// @}
 
 
 template <class Grid>
@@ -367,6 +405,8 @@ void VtkReader<Grid>::readAppended (std::ifstream& input, std::vector<T>& values
   } else {
     input.read((char*)&size, sizeof(std::uint64_t));
   }
+  std::cout << "size = " << size << "\n";
+  assert(size > 0 && (size % sizeof(T)) == 0);
   values.resize(size / sizeof(T));
 
   if (format_ == Vtk::COMPRESSED) {
@@ -381,13 +421,14 @@ void VtkReader<Grid>::readAppended (std::ifstream& input, std::vector<T>& values
       read_compressed(values.data() + i*num_values, buffer_in.data(), bs, cbs[i], input);
     }
   } else {
-    input.read((char*)(values.data()), size);
+    input.readsome((char*)(values.data()), size);
+    assert(input.gcount() == size);
   }
 }
 
 
 template <class Grid>
-void VtkReader<Grid>::createGrid(GridFactory<Grid>& factory) const
+void VtkReader<Grid>::createGrid () const
 {
   assert(vec_types.size() == vec_offsets.size());
   std::size_t idx = 0;
@@ -408,14 +449,15 @@ void VtkReader<Grid>::createGrid(GridFactory<Grid>& factory) const
     for (std::size_t j = 0; j < nNodes; ++j)
       cell[j] = vtk_cell[cellType.localIndex(j)];
 
-    factory.insertElement(type,cell);
+    factory_->insertElement(type,cell);
   }
 }
 
 
 template <class Grid>
-std::map<std::string, std::string> VtkReader<Grid>::parseXml(std::string const& line)
+std::map<std::string, std::string> VtkReader<Grid>::parseXml (std::string const& line, bool& closed)
 {
+  closed = false;
   std::map<std::string, std::string> attr;
 
   Sections sec = NO_SECTION;
@@ -430,6 +472,8 @@ std::map<std::string, std::string> VtkReader<Grid>::parseXml(std::string const&
         name.clear();
         sec = XML_NAME;
         name.push_back(c);
+      } else if (c == '/') {
+        closed = true;
       }
       break;
     case XML_NAME:
diff --git a/dune/vtk/vtktypes.cc b/dune/vtk/vtktypes.cc
index 6351ba3ca0dfa8ca088aa43bd4b33b79e5457dbc..2e5acd4df4bf2853c048d145d57ee913f17a0ad3 100644
--- a/dune/vtk/vtktypes.cc
+++ b/dune/vtk/vtktypes.cc
@@ -5,7 +5,7 @@
 namespace Dune {
 namespace Vtk {
 
-std::map<size_t, GeometryType> Map::type = {
+std::map<std::uint8_t, GeometryType> Map::type = {
   { 1, GeometryTypes::vertex },
   { 3, GeometryTypes::line },
   { 5, GeometryTypes::triangle },
diff --git a/dune/vtk/vtktypes.hh b/dune/vtk/vtktypes.hh
index dbab2b27e028a428613b05e56d2e60135fe7de27..f0567b52c5c7e911f6c572d80a23affeb058e1d6 100644
--- a/dune/vtk/vtktypes.hh
+++ b/dune/vtk/vtktypes.hh
@@ -28,64 +28,6 @@ namespace Dune
       FLOAT64 = 64
     };
 
-    inline std::size_t size(DataTypes type)
-    {
-      switch (type) {
-        case INT8:    return 8;
-        case UINT8:   return 8;
-        case INT16:   return 16;
-        case UINT16:  return 16;
-        case INT32:   return 32;
-        case UINT32:  return 32;
-        case INT64:   return 64;
-        case UINT64:  return 64;
-        case FLOAT32: return 32;
-        case FLOAT64: return 64;
-        default:
-          std::abort();
-          return 0;
-      }
-    }
-
-    template <class T>
-    void cast_value(char const* data, DataTypes type, T& value)
-    {
-      switch (type) {
-        case INT8:
-          value = *reinterpret_cast<std::int8_t const*>(data);
-          break;
-        case UINT8:
-          value = *reinterpret_cast<std::uint8_t const*>(data);
-          break;
-        case INT16:
-          value = *reinterpret_cast<std::int16_t const*>(data);
-          break;
-        case UINT16:
-          value = *reinterpret_cast<std::uint16_t const*>(data);
-          break;
-        case INT32:
-          value = *reinterpret_cast<std::int32_t const*>(data);
-          break;
-        case UINT32:
-          value = *reinterpret_cast<std::uint32_t const*>(data);
-          break;
-        case INT64:
-          value = *reinterpret_cast<std::int64_t const*>(data);
-          break;
-        case UINT64:
-          value = *reinterpret_cast<std::uint64_t const*>(data);
-          break;
-        case FLOAT32:
-          value = *reinterpret_cast<float const*>(data);
-          break;
-        case FLOAT64:
-          value = *reinterpret_cast<double const*>(data);
-          break;
-        default:
-          std::abort();
-      }
-    }
-
     enum PositionTypes {
       VERTEX_DATA,
       CELL_DATA
@@ -98,8 +40,8 @@ namespace Dune
 
     struct Map
     {
-      static std::map<std::size_t, GeometryType> type;
-      static std::map<std::string, DataTypes> datatype;
+      static std::map<std::uint8_t, GeometryType> type; // VTK Cell type -> Dune::GeometryType
+      static std::map<std::string, DataTypes> datatype; // String -> DataTypes
     };
 
 
@@ -109,11 +51,13 @@ namespace Dune
     public:
       CellType (GeometryType const& t);
 
+      /// Return VTK Cell type
       std::uint8_t type () const
       {
         return type_;
       }
 
+      /// Return a permutation of Dune elemenr vertices to conform to VTK element numbering
       int localIndex (int idx) const
       {
         return permutation_[idx];
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
index 9021509147e5bc1f28fa9ce3ecfd40cbe79dc608..7b64d61ab2270847069d4c2f37122846c7203586 100644
--- a/lib/CMakeLists.txt
+++ b/lib/CMakeLists.txt
@@ -5,9 +5,7 @@ if(BUILD_SHARED_LIBS)
   set(_OBJECT_FLAG "")
 endif()
 
-dune_add_library(dunevtkwriter
+dune_add_library(dunevtk
   _DUNE_TARGET_OBJECTS:filesystem_
   _DUNE_TARGET_OBJECTS:vtktypes_
   ADD_LIBS ${DUNE_LIBS})
-
-# list(INSERT DUNE_DEFAULT_LIBS "dunevtkwriter")
\ No newline at end of file
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index c7fd6cd5bed29a168d6deaa51cf08442b00f3853..975fc79b486e14a55ac71f75306be1c7d8683fc5 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,7 +1,7 @@
 add_executable("vtkwriter" vtkwriter.cc)
 target_link_dune_default_libraries("vtkwriter")
-target_link_libraries("vtkwriter" dunevtkwriter)
+target_link_libraries("vtkwriter" dunevtk)
 
 add_executable("vtkreader" vtkreader.cc)
 target_link_dune_default_libraries("vtkreader")
-target_link_libraries("vtkreader" dunevtkwriter)
\ No newline at end of file
+target_link_libraries("vtkreader" dunevtk)
\ No newline at end of file
diff --git a/src/vtkreader.cc b/src/vtkreader.cc
index 6b0506ad666e1b2b1d9b97fbf36d3533bef62a61..8f8db803075e7322fe24673139b531d2e04b5a91 100644
--- a/src/vtkreader.cc
+++ b/src/vtkreader.cc
@@ -31,7 +31,7 @@ int main(int argc, char** argv)
   {
     FieldVector<double,dim> lowerLeft; lowerLeft = 0.0;
     FieldVector<double,dim> upperRight; upperRight = 1.0;
-    auto numElements = filledArray<dim,unsigned int>(1);
+    auto numElements = filledArray<dim,unsigned int>(4);
     auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);
     auto& grid = *gridPtr;
 
@@ -44,8 +44,7 @@ int main(int argc, char** argv)
   }
 
   {
-    VtkReader<GridType> vtkReader{};
-    auto gridPtr = vtkReader.read("test_ascii_float32.vtu");
+    auto gridPtr = VtkReader<GridType>::read("test_ascii_float32.vtu");
     auto& grid = *gridPtr;
 
     VtkWriter<GridView> vtkWriter(grid.leafGridView());
@@ -53,8 +52,7 @@ int main(int argc, char** argv)
   }
 
   {
-    VtkReader<GridType> vtkReader{};
-    auto gridPtr = vtkReader.read("test_binary_float32.vtu");
+    auto gridPtr = VtkReader<GridType>::read("test_binary_float32.vtu");
     auto& grid = *gridPtr;
 
     VtkWriter<GridView> vtkWriter(grid.leafGridView());
@@ -62,11 +60,20 @@ int main(int argc, char** argv)
   }
 
   {
-    VtkReader<GridType> vtkReader{};
-    auto gridPtr = vtkReader.read("test_compressed_float64.vtu");
+    auto gridPtr = VtkReader<GridType>::read("test_compressed_float64.vtu");
     auto& grid = *gridPtr;
 
     VtkWriter<GridView> vtkWriter(grid.leafGridView());
     vtkWriter.write("test_ascii_float64_3.vtu", Vtk::ASCII, Vtk::FLOAT64);
   }
+
+  {
+    using GridType2d = UGGrid<2>;
+    using GridView2d = typename GridType2d::LeafGridView;
+    auto gridPtr = VtkReader<GridType2d>::read("paraview_2d.vtu");
+    auto& grid = *gridPtr;
+
+    VtkWriter<GridView2d> vtkWriter(grid.leafGridView());
+    vtkWriter.write("paraview_2d_ascii.vtu", Vtk::ASCII, Vtk::FLOAT64);
+  }
 }
diff --git a/src/vtkwriter.cc b/src/vtkwriter.cc
index 1969bd1dee5b22670ea90fa81bb77d2e4b2fcc72..9275963856b08d47f419506121bfed414581ca68 100644
--- a/src/vtkwriter.cc
+++ b/src/vtkwriter.cc
@@ -36,13 +36,13 @@ int main(int argc, char** argv)
 #if GRID_TYPE == 1
   using GridType = YaspGrid<dim>;
   FieldVector<double,dim> upperRight; upperRight = 1.0;
-  auto numElements = filledArray<dim,int>(1);
+  auto numElements = filledArray<dim,int>(4);
   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>(2);
+  auto numElements = filledArray<dim,unsigned int>(4);
   auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);
   auto& grid = *gridPtr;
 #endif