data_out_base.cc

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2022 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
 
// TODO: Do neighbors for dx and povray smooth triangles
 
//--------------------------------------------------------------------
// Remarks on the implementations
//
// Variable names: in most functions, variable names have been
// standardized in the following way:
//
// n1, n2, ni Number of points in coordinate direction 1, 2, i
//    will be 1 if i>=dim
//
// i1, i2, ii Loop variable running up to ni
//
// d1, d2, di Multiplicators for ii to find positions in the
//    array of nodes.
//--------------------------------------------------------------------
 
#include <deal.II/base/data_out_base.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/mpi_large_count.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/utilities.h>
 
#include <deal.II/numerics/data_component_interpretation.h>
 
#include <algorithm>
#include <cmath>
#include <cstring>
#include <ctime>
#include <fstream>
#include <iomanip>
#include <memory>
#include <set>
#include <sstream>
 
// we use std::uint32_t and std::uint8_t below, which are declared here:
#include <cstdint>
#include <vector>
 
#ifdef DEAL_II_WITH_ZLIB
#  include <zlib.h>
#endif
 
#ifdef DEAL_II_WITH_HDF5
#  include <hdf5.h>
#endif
 
DEAL_II_NAMESPACE_OPEN
 
 
// we need the following exception from a global function, so can't declare it
// in the usual way inside a class
namespace
{
  DeclException2(ExcUnexpectedInput,
                 std::string,
                 std::string,
                 << "Unexpected input: expected line\n  <" << arg1
                 << ">\nbut got\n  <" << arg2 << ">");
}
 
 
namespace
{
#ifdef DEAL_II_WITH_ZLIB
  int
  get_zlib_compression_level(
    const DataOutBase::VtkFlags::ZlibCompressionLevel level)
  {
    switch (level)
      {
        case (DataOutBase::VtkFlags::no_compression):
          return Z_NO_COMPRESSION;
        case (DataOutBase::VtkFlags::best_speed):
          return Z_BEST_SPEED;
        case (DataOutBase::VtkFlags::best_compression):
          return Z_BEST_COMPRESSION;
        case (DataOutBase::VtkFlags::default_compression):
          return Z_DEFAULT_COMPRESSION;
        default:
          Assert(false, ExcNotImplemented());
          return Z_NO_COMPRESSION;
      }
  }
 
  template <typename T>
  void
  write_compressed_block(const std::vector<T> &       data,
                         const DataOutBase::VtkFlags &flags,
                         std::ostream &               output_stream)
  {
    if (data.size() != 0)
      {
        const std::size_t uncompressed_size = (data.size() * sizeof(T));
 
        // While zlib's compress2 uses unsigned long (which is 64bits
        // on Linux), the vtu compression header stores the block size
        // as an std::uint32_t (see below). While we could implement
        // writing several smaller blocks, we haven't done that. Let's
        // trigger an error for the user instead:
        AssertThrow(uncompressed_size <=
                      std::numeric_limits<std::uint32_t>::max(),
                    ExcNotImplemented());
 
        // allocate a buffer for compressing data and do so
        auto compressed_data_length = compressBound(uncompressed_size);
        AssertThrow(compressed_data_length <=
                      std::numeric_limits<std::uint32_t>::max(),
                    ExcNotImplemented());
 
        std::vector<unsigned char> compressed_data(compressed_data_length);
 
        int err =
          compress2(&compressed_data[0],
                    &compressed_data_length,
                    reinterpret_cast<const Bytef *>(data.data()),
                    uncompressed_size,
                    get_zlib_compression_level(flags.compression_level));
        (void)err;
        Assert(err == Z_OK, ExcInternalError());
 
        // Discard the unnecessary bytes
        compressed_data.resize(compressed_data_length);
 
        // now encode the compression header
        const std::uint32_t compression_header[4] = {
          1,                                             /* number of blocks */
          static_cast<std::uint32_t>(uncompressed_size), /* size of block */
          static_cast<std::uint32_t>(
            uncompressed_size), /* size of last block */
          static_cast<std::uint32_t>(
            compressed_data_length)}; /* list of compressed sizes of blocks */
 
        const auto header_start =
          reinterpret_cast<const unsigned char *>(&compression_header[0]);
 
        output_stream << Utilities::encode_base64(
                           {header_start,
                            header_start + 4 * sizeof(std::uint32_t)})
                      << Utilities::encode_base64(compressed_data);
      }
  }
#endif
} // namespace
 
 
// some declarations of functions and locally used classes
namespace DataOutBase
{
  namespace
  {
    class SvgCell
    {
    public:
      // Center of the cell (three-dimensional)
      Point<3> center;
 
      Point<3> vertices[4];
 
      float depth;
 
      Point<2> projected_vertices[4];
 
      // Center of the cell (projected, two-dimensional)
      Point<2> projected_center;
 
      bool
      operator<(const SvgCell &) const;
    };
 
    bool
    SvgCell::operator<(const SvgCell &e) const
    {
      // note the "wrong" order in which we sort the elements
      return depth > e.depth;
    }
 
 
 
    class EpsCell2d
    {
    public:
      Point<2> vertices[4];
 
      float color_value;
 
      float depth;
 
      bool
      operator<(const EpsCell2d &) const;
    };
 
    bool
    EpsCell2d::operator<(const EpsCell2d &e) const
    {
      // note the "wrong" order in which we sort the elements
      return depth > e.depth;
    }
 
 
 
    template <int dim, int spacedim, typename Number = double>
    void
    write_gmv_reorder_data_vectors(
      const std::vector<Patch<dim, spacedim>> &patches,
      Table<2, Number> &                       data_vectors)
    {
      // If there is nothing to write, just return
      if (patches.size() == 0)
        return;
 
      // unlike in the main function, we don't have here the data_names field,
      // so we initialize it with the number of data sets in the first patch.
      // the equivalence of these two definitions is checked in the main
      // function.
 
      // we have to take care, however, whether the points are appended to the
      // end of the patch.data table
      const unsigned int n_data_sets = patches[0].points_are_available ?
                                         (patches[0].data.n_rows() - spacedim) :
                                         patches[0].data.n_rows();
 
      Assert(data_vectors.size()[0] == n_data_sets, ExcInternalError());
 
      // loop over all patches
      unsigned int next_value = 0;
      for (const auto &patch : patches)
        {
          const unsigned int n_subdivisions = patch.n_subdivisions;
          (void)n_subdivisions;
 
          Assert((patch.data.n_rows() == n_data_sets &&
                  !patch.points_are_available) ||
                   (patch.data.n_rows() == n_data_sets + spacedim &&
                    patch.points_are_available),
                 ExcDimensionMismatch(patch.points_are_available ?
                                        (n_data_sets + spacedim) :
                                        n_data_sets,
                                      patch.data.n_rows()));
          Assert(patch.reference_cell != ReferenceCells::get_hypercube<dim>() ||
                   (n_data_sets == 0) ||
                   (patch.data.n_cols() ==
                    Utilities::fixed_power<dim>(n_subdivisions + 1)),
                 ExcInvalidDatasetSize(patch.data.n_cols(),
                                       n_subdivisions + 1));
 
          for (unsigned int i = 0; i < patch.data.n_cols(); ++i, ++next_value)
            for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
              data_vectors[data_set][next_value] = patch.data(data_set, i);
        }
 
      for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
        Assert(data_vectors[data_set].size() == next_value, ExcInternalError());
    }
  } // namespace
 
 
 
  DataOutFilter::DataOutFilter()
    : flags(false, true)
    , node_dim(numbers::invalid_unsigned_int)
    , num_cells(numbers::invalid_unsigned_int)
  {}
 
 
 
  DataOutFilter::DataOutFilter(const DataOutBase::DataOutFilterFlags &flags)
    : flags(flags)
    , node_dim(numbers::invalid_unsigned_int)
    , num_cells(numbers::invalid_unsigned_int)
  {}
 
 
 
  template <int dim>
  void
  DataOutFilter::write_point(const unsigned int index, const Point<dim> &p)
  {
    node_dim = dim;
 
    Point<3> int_pt;
    for (unsigned int d = 0; d < dim; ++d)
      int_pt(d) = p(d);
 
    const Map3DPoint::const_iterator it = existing_points.find(int_pt);
    unsigned int                     internal_ind;
 
    // If the point isn't in the set, or we're not filtering duplicate points,
    // add it
    if (it == existing_points.end() || !flags.filter_duplicate_vertices)
      {
        internal_ind = existing_points.size();
        existing_points.insert(std::make_pair(int_pt, internal_ind));
      }
    else
      {
        internal_ind = it->second;
      }
    // Now add the index to the list of filtered points
    filtered_points[index] = internal_ind;
  }
 
 
 
  void
  DataOutFilter::internal_add_cell(const unsigned int cell_index,
                                   const unsigned int pt_index)
  {
    filtered_cells[cell_index] = filtered_points[pt_index];
 
    // (Re)-initialize counter at any first call to this method.
    if (cell_index == 0)
      num_cells = 1;
  }
 
 
 
  void
  DataOutFilter::fill_node_data(std::vector<double> &node_data) const
  {
    node_data.resize(existing_points.size() * node_dim);
 
    for (const auto &existing_point : existing_points)
      {
        for (unsigned int d = 0; d < node_dim; ++d)
          node_data[node_dim * existing_point.second + d] =
            existing_point.first(d);
      }
  }
 
 
 
  void
  DataOutFilter::fill_cell_data(const unsigned int         local_node_offset,
                                std::vector<unsigned int> &cell_data) const
  {
    cell_data.resize(filtered_cells.size());
 
    for (const auto &filtered_cell : filtered_cells)
      {
        cell_data[filtered_cell.first] =
          filtered_cell.second + local_node_offset;
      }
  }
 
 
 
  std::string
  DataOutFilter::get_data_set_name(const unsigned int set_num) const
  {
    return data_set_names.at(set_num);
  }
 
 
 
  unsigned int
  DataOutFilter::get_data_set_dim(const unsigned int set_num) const
  {
    return data_set_dims.at(set_num);
  }
 
 
 
  const double *
  DataOutFilter::get_data_set(const unsigned int set_num) const
  {
    return data_sets[set_num].data();
  }
 
 
 
  unsigned int
  DataOutFilter::n_nodes() const
  {
    return existing_points.size();
  }
 
 
 
  unsigned int
  DataOutFilter::n_cells() const
  {
    return num_cells;
  }
 
 
 
  unsigned int
  DataOutFilter::n_data_sets() const
  {
    return data_set_names.size();
  }
 
 
 
  void
  DataOutFilter::flush_points()
  {}
 
 
 
  void
  DataOutFilter::flush_cells()
  {}
 
 
 
  template <int dim>
  void
  DataOutFilter::write_cell(const unsigned int index,
                            const unsigned int start,
                            const unsigned int d1,
                            const unsigned int d2,
                            const unsigned int d3)
  {
    ++num_cells;
 
    const unsigned int base_entry =
      index * GeometryInfo<dim>::vertices_per_cell;
 
    internal_add_cell(base_entry + 0, start);
    if (dim >= 1)
      {
        internal_add_cell(base_entry + 1, start + d1);
        if (dim >= 2)
          {
            internal_add_cell(base_entry + 2, start + d2 + d1);
            internal_add_cell(base_entry + 3, start + d2);
            if (dim >= 3)
              {
                internal_add_cell(base_entry + 4, start + d3);
                internal_add_cell(base_entry + 5, start + d3 + d1);
                internal_add_cell(base_entry + 6, start + d3 + d2 + d1);
                internal_add_cell(base_entry + 7, start + d3 + d2);
              }
          }
      }
  }
 
 
 
  void
  DataOutFilter::write_cell_single(const unsigned int   index,
                                   const unsigned int   start,
                                   const unsigned int   n_points,
                                   const ReferenceCell &reference_cell)
  {
    ++num_cells;
 
    const unsigned int base_entry = index * n_points;
 
    static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
 
    for (unsigned int i = 0; i < n_points; ++i)
      internal_add_cell(base_entry + i,
                        start + (reference_cell == ReferenceCells::Pyramid ?
                                   table[i] :
                                   i));
  }
 
 
 
  void
  DataOutFilter::write_data_set(const std::string &     name,
                                const unsigned int      dimension,
                                const unsigned int      set_num,
                                const Table<2, double> &data_vectors)
  {
    unsigned int new_dim;
 
    // HDF5/XDMF output only supports 1D or 3D output, so force rearrangement if
    // needed
    if (flags.xdmf_hdf5_output && dimension != 1)
      new_dim = 3;
    else
      new_dim = dimension;
 
    // Record the data set name, dimension, and allocate space for it
    data_set_names.push_back(name);
    data_set_dims.push_back(new_dim);
    data_sets.emplace_back(new_dim * existing_points.size());
 
    // TODO: averaging, min/max, etc for merged vertices
    for (unsigned int i = 0; i < filtered_points.size(); ++i)
      {
        const unsigned int r = filtered_points[i];
 
        for (unsigned int d = 0; d < new_dim; ++d)
          {
            if (d < dimension)
              data_sets.back()[r * new_dim + d] = data_vectors(set_num + d, i);
            else
              data_sets.back()[r * new_dim + d] = 0;
          }
      }
  }
} // namespace DataOutBase
 
 
 
//----------------------------------------------------------------------//
// Auxiliary data
//----------------------------------------------------------------------//
 
namespace
{
  const char *gmv_cell_type[4] = {"", "line 2", "quad 4", "hex 8"};
 
  const char *ucd_cell_type[4] = {"pt", "line", "quad", "hex"};
 
  const char *tecplot_cell_type[4] = {"", "lineseg", "quadrilateral", "brick"};
 
#ifdef DEAL_II_HAVE_TECPLOT
  const unsigned int tecplot_binary_cell_type[4] = {0, 0, 1, 3};
#endif
 
  // Define cell id using VTK nomenclature for linear, quadratic and
  // high-order Lagrange cells
  enum vtk_linear_cell_type
  {
    VTK_VERTEX     = 1,
    VTK_LINE       = 3,
    VTK_TRIANGLE   = 5,
    VTK_QUAD       = 9,
    VTK_TETRA      = 10,
    VTK_HEXAHEDRON = 12,
    VTK_WEDGE      = 13,
    VTK_PYRAMID    = 14
  };
 
  enum vtk_quadratic_cell_type
  {
    VTK_QUADRATIC_EDGE       = 21,
    VTK_QUADRATIC_TRIANGLE   = 22,
    VTK_QUADRATIC_QUAD       = 23,
    VTK_QUADRATIC_TETRA      = 24,
    VTK_QUADRATIC_HEXAHEDRON = 25,
    VTK_QUADRATIC_WEDGE      = 26,
    VTK_QUADRATIC_PYRAMID    = 27
  };
 
  enum vtk_lagrange_cell_type
  {
    VTK_LAGRANGE_CURVE         = 68,
    VTK_LAGRANGE_TRIANGLE      = 69,
    VTK_LAGRANGE_QUADRILATERAL = 70,
    VTK_LAGRANGE_TETRAHEDRON   = 71,
    VTK_LAGRANGE_HEXAHEDRON    = 72,
    VTK_LAGRANGE_WEDGE         = 73,
    VTK_LAGRANGE_PYRAMID       = 74
  };
 
  template <int dim, int spacedim>
  std::array<unsigned int, 3>
  extract_vtk_patch_info(const DataOutBase::Patch<dim, spacedim> &patch,
                         const bool write_higher_order_cells)
  {
    std::array<unsigned int, 3> vtk_cell_id{};
 
    if (write_higher_order_cells)
      {
        if (patch.reference_cell == ReferenceCells::get_hypercube<dim>())
          {
            const std::array<unsigned int, 4> cell_type_by_dim{
              {VTK_VERTEX,
               VTK_LAGRANGE_CURVE,
               VTK_LAGRANGE_QUADRILATERAL,
               VTK_LAGRANGE_HEXAHEDRON}};
            vtk_cell_id[0] = cell_type_by_dim[dim];
            vtk_cell_id[1] = 1;
          }
        else if (patch.reference_cell == ReferenceCells::Triangle)
          {
            vtk_cell_id[0] = VTK_LAGRANGE_TRIANGLE;
            vtk_cell_id[1] = 1;
          }
        else
          {
            Assert(false, ExcNotImplemented());
          }
      }
    else if (patch.reference_cell == ReferenceCells::Triangle &&
             patch.data.n_cols() == 3)
      {
        vtk_cell_id[0] = VTK_TRIANGLE;
        vtk_cell_id[1] = 1;
      }
    else if (patch.reference_cell == ReferenceCells::Triangle &&
             patch.data.n_cols() == 6)
      {
        vtk_cell_id[0] = VTK_QUADRATIC_TRIANGLE;
        vtk_cell_id[1] = 1;
      }
    else if (patch.reference_cell == ReferenceCells::Tetrahedron &&
             patch.data.n_cols() == 4)
      {
        vtk_cell_id[0] = VTK_TETRA;
        vtk_cell_id[1] = 1;
      }
    else if (patch.reference_cell == ReferenceCells::Tetrahedron &&
             patch.data.n_cols() == 10)
      {
        vtk_cell_id[0] = VTK_QUADRATIC_TETRA;
        vtk_cell_id[1] = 1;
      }
    else if (patch.reference_cell == ReferenceCells::Wedge &&
             patch.data.n_cols() == 6)
      {
        vtk_cell_id[0] = VTK_WEDGE;
        vtk_cell_id[1] = 1;
      }
    else if (patch.reference_cell == ReferenceCells::Pyramid &&
             patch.data.n_cols() == 5)
      {
        vtk_cell_id[0] = VTK_PYRAMID;
        vtk_cell_id[1] = 1;
      }
    else if (patch.reference_cell == ReferenceCells::get_hypercube<dim>())
      {
        const std::array<unsigned int, 4> cell_type_by_dim{
          {VTK_VERTEX, VTK_LINE, VTK_QUAD, VTK_HEXAHEDRON}};
        vtk_cell_id[0] = cell_type_by_dim[dim];
        vtk_cell_id[1] = Utilities::pow(patch.n_subdivisions, dim);
      }
    else
      {
        Assert(false, ExcNotImplemented());
      }
 
    if (patch.reference_cell != ReferenceCells::get_hypercube<dim>() ||
        write_higher_order_cells)
      vtk_cell_id[2] = patch.data.n_cols();
    else
      vtk_cell_id[2] = GeometryInfo<dim>::vertices_per_cell;
 
    return vtk_cell_id;
  }
 
  //----------------------------------------------------------------------//
  // Auxiliary functions
  //----------------------------------------------------------------------//
 
  // For a given patch that corresponds to a hypercube cell, compute the
  // location of a node interpolating the corner nodes linearly
  // at the point lattice_location/n_subdivisions where lattice_location
  // is a dim-dimensional integer vector. If the points are
  // saved in the patch.data member, return the saved point instead.
  template <int dim, int spacedim>
  inline Point<spacedim>
  get_equispaced_location(
    const DataOutBase::Patch<dim, spacedim> &  patch,
    const std::initializer_list<unsigned int> &lattice_location,
    const unsigned int                         n_subdivisions)
  {
    // This function only makes sense when called on hypercube cells
    Assert(patch.reference_cell == ReferenceCells::get_hypercube<dim>(),
           ExcInternalError());
 
    Assert(lattice_location.size() == dim, ExcInternalError());
 
    const unsigned int xstep = (dim > 0 ? *(lattice_location.begin() + 0) : 0);
    const unsigned int ystep = (dim > 1 ? *(lattice_location.begin() + 1) : 0);
    const unsigned int zstep = (dim > 2 ? *(lattice_location.begin() + 2) : 0);
 
    // If the patch stores the locations of nodes (rather than of only the
    // vertices), then obtain the location by direct lookup.
    if (patch.points_are_available)
      {
        Assert(n_subdivisions == patch.n_subdivisions, ExcNotImplemented());
 
        unsigned int point_no = 0;
        switch (dim)
          {
            case 3:
              AssertIndexRange(zstep, n_subdivisions + 1);
              point_no += (n_subdivisions + 1) * (n_subdivisions + 1) * zstep;
              DEAL_II_FALLTHROUGH;
            case 2:
              AssertIndexRange(ystep, n_subdivisions + 1);
              point_no += (n_subdivisions + 1) * ystep;
              DEAL_II_FALLTHROUGH;
            case 1:
              AssertIndexRange(xstep, n_subdivisions + 1);
              point_no += xstep;
              DEAL_II_FALLTHROUGH;
            case 0:
              // break here for dim<=3
              break;
 
            default:
              Assert(false, ExcNotImplemented());
          }
        Point<spacedim> node;
        for (unsigned int d = 0; d < spacedim; ++d)
          node[d] = patch.data(patch.data.size(0) - spacedim + d, point_no);
        return node;
      }
    else
      // The patch does not store node locations, so we have to interpolate
      // between its vertices:
      {
        if (dim == 0)
          return patch.vertices[0];
        else
          {
            // perform a dim-linear interpolation
            const double stepsize = 1. / n_subdivisions;
            const double xfrac    = xstep * stepsize;
 
            Point<spacedim> node =
              (patch.vertices[1] * xfrac) + (patch.vertices[0] * (1 - xfrac));
            if (dim > 1)
              {
                const double yfrac = ystep * stepsize;
                node *= 1 - yfrac;
                node += ((patch.vertices[3] * xfrac) +
                         (patch.vertices[2] * (1 - xfrac))) *
                        yfrac;
                if (dim > 2)
                  {
                    const double zfrac = zstep * stepsize;
                    node *= (1 - zfrac);
                    node += (((patch.vertices[5] * xfrac) +
                              (patch.vertices[4] * (1 - xfrac))) *
                               (1 - yfrac) +
                             ((patch.vertices[7] * xfrac) +
                              (patch.vertices[6] * (1 - xfrac))) *
                               yfrac) *
                            zfrac;
                  }
              }
            return node;
          }
      }
  }
 
  // For a given patch, compute the nodes for arbitrary (non-hypercube) cells.
  // If the points are saved in the patch.data member, return the saved point
  // instead.
  template <int dim, int spacedim>
  inline Point<spacedim>
  get_node_location(const DataOutBase::Patch<dim, spacedim> &patch,
                    const unsigned int                       node_index)
  {
    // Due to a historical accident, we are using a different indexing
    // for pyramids in this file than we do where we create patches.
    // So translate if necessary.
    unsigned int point_no_actual = node_index;
    if (patch.reference_cell == ReferenceCells::Pyramid)
      {
        AssertDimension(patch.n_subdivisions, 1);
 
        static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
        point_no_actual                                = table[node_index];
      }
 
    // If the patch stores the locations of nodes (rather than of only the
    // vertices), then obtain the location by direct lookup.
    if (patch.points_are_available)
      {
        Point<spacedim> node;
        for (unsigned int d = 0; d < spacedim; ++d)
          node[d] =
            patch.data(patch.data.size(0) - spacedim + d, point_no_actual);
        return node;
      }
    else
      // The patch does not store node locations, so we have to interpolate
      // between its vertices. This isn't currently implemented for anything
      // other than one subdivision, but would go here.
      //
      // For n_subdivisions==1, the locations are simply those of vertices, so
      // get the information from there.
      {
        AssertDimension(patch.n_subdivisions, 1);
 
        return patch.vertices[point_no_actual];
      }
  }
 
 
 
  int
  vtk_point_index_from_ijk(const unsigned i,
                           const unsigned j,
                           const unsigned,
                           const std::array<unsigned, 2> &order)
  {
    const bool ibdy = (i == 0 || i == order[0]);
    const bool jbdy = (j == 0 || j == order[1]);
    // How many boundaries do we lie on at once?
    const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0);
 
    if (nbdy == 2) // Vertex DOF
      { // ijk is a corner node. Return the proper index (somewhere in [0,3]):
        return (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0));
      }
 
    int offset = 4;
    if (nbdy == 1) // Edge DOF
      {
        if (!ibdy)
          { // On i axis
            return (i - 1) + (j != 0u ? order[0] - 1 + order[1] - 1 : 0) +
                   offset;
          }
 
        if (!jbdy)
          { // On j axis
            return (j - 1) +
                   (i != 0u ? order[0] - 1 :
                              2 * (order[0] - 1) + order[1] - 1) +
                   offset;
          }
      }
 
    offset += 2 * (order[0] - 1 + order[1] - 1);
    // nbdy == 0: Face DOF
    return offset + (i - 1) + (order[0] - 1) * ((j - 1));
  }
 
 
 
  int
  vtk_point_index_from_ijk(const unsigned                 i,
                           const unsigned                 j,
                           const unsigned                 k,
                           const std::array<unsigned, 3> &order)
  {
    const bool ibdy = (i == 0 || i == order[0]);
    const bool jbdy = (j == 0 || j == order[1]);
    const bool kbdy = (k == 0 || k == order[2]);
    // How many boundaries do we lie on at once?
    const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0);
 
    if (nbdy == 3) // Vertex DOF
      { // ijk is a corner node. Return the proper index (somewhere in [0,7]):
        return (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0)) +
               (k != 0u ? 4 : 0);
      }
 
    int offset = 8;
    if (nbdy == 2) // Edge DOF
      {
        if (!ibdy)
          { // On i axis
            return (i - 1) + (j != 0u ? order[0] - 1 + order[1] - 1 : 0) +
                   (k != 0u ? 2 * (order[0] - 1 + order[1] - 1) : 0) + offset;
          }
        if (!jbdy)
          { // On j axis
            return (j - 1) +
                   (i != 0u ? order[0] - 1 :
                              2 * (order[0] - 1) + order[1] - 1) +
                   (k != 0u ? 2 * (order[0] - 1 + order[1] - 1) : 0) + offset;
          }
        // !kbdy, On k axis
        offset += 4 * (order[0] - 1) + 4 * (order[1] - 1);
        return (k - 1) +
               (order[2] - 1) *
                 (i != 0u ? (j != 0u ? 3 : 1) : (j != 0u ? 2 : 0)) +
               offset;
      }
 
    offset += 4 * (order[0] - 1 + order[1] - 1 + order[2] - 1);
    if (nbdy == 1) // Face DOF
      {
        if (ibdy) // On i-normal face
          {
            return (j - 1) + ((order[1] - 1) * (k - 1)) +
                   (i != 0u ? (order[1] - 1) * (order[2] - 1) : 0) + offset;
          }
        offset += 2 * (order[1] - 1) * (order[2] - 1);
        if (jbdy) // On j-normal face
          {
            return (i - 1) + ((order[0] - 1) * (k - 1)) +
                   (j != 0u ? (order[2] - 1) * (order[0] - 1) : 0) + offset;
          }
        offset += 2 * (order[2] - 1) * (order[0] - 1);
        // kbdy, On k-normal face
        return (i - 1) + ((order[0] - 1) * (j - 1)) +
               (k != 0u ? (order[0] - 1) * (order[1] - 1) : 0) + offset;
      }
 
    // nbdy == 0: Body DOF
    offset +=
      2 * ((order[1] - 1) * (order[2] - 1) + (order[2] - 1) * (order[0] - 1) +
           (order[0] - 1) * (order[1] - 1));
    return offset + (i - 1) +
           (order[0] - 1) * ((j - 1) + (order[1] - 1) * ((k - 1)));
  }
 
 
 
  int
  vtk_point_index_from_ijk(const unsigned,
                           const unsigned,
                           const unsigned,
                           const std::array<unsigned, 0> &)
  {
    Assert(false, ExcNotImplemented());
    return 0;
  }
 
 
 
  int
  vtk_point_index_from_ijk(const unsigned,
                           const unsigned,
                           const unsigned,
                           const std::array<unsigned, 1> &)
  {
    Assert(false, ExcNotImplemented());
    return 0;
  }
 
 
 
  template <int dim, int spacedim>
  static void
  compute_sizes(const std::vector<DataOutBase::Patch<dim, spacedim>> &patches,
                unsigned int &                                        n_nodes,
                unsigned int &                                        n_cells)
  {
    n_nodes = 0;
    n_cells = 0;
    for (const auto &patch : patches)
      {
        Assert(patch.reference_cell != ReferenceCells::Invalid,
               ExcMessage(
                 "The reference cell for this patch is set to 'Invalid', "
                 "but that is clearly not a valid choice. Did you forget "
                 "to set the reference cell for the patch?"));
 
        // The following formula doesn't hold for non-tensor products.
        if (patch.reference_cell == ReferenceCells::get_hypercube<dim>())
          {
            n_nodes += Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
            n_cells += Utilities::fixed_power<dim>(patch.n_subdivisions);
          }
        else
          {
            Assert(patch.n_subdivisions == 1, ExcNotImplemented());
            n_nodes += patch.reference_cell.n_vertices();
            n_cells += 1;
          }
      }
  }
 
 
 
  template <int dim, int spacedim>
  static void
  compute_sizes(const std::vector<DataOutBase::Patch<dim, spacedim>> &patches,
                const bool    write_higher_order_cells,
                unsigned int &n_nodes,
                unsigned int &n_cells,
                unsigned int &n_points_and_n_cells)
  {
    n_nodes              = 0;
    n_cells              = 0;
    n_points_and_n_cells = 0;
 
    for (const auto &patch : patches)
      {
        // The following formulas don't hold for non-tensor products.
        if (patch.reference_cell == ReferenceCells::get_hypercube<dim>())
          {
            n_nodes += Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
 
            if (write_higher_order_cells)
              {
                n_cells += 1;
                n_points_and_n_cells +=
                  1 + Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
              }
            else
              {
                n_cells += Utilities::fixed_power<dim>(patch.n_subdivisions);
                n_points_and_n_cells +=
                  Utilities::fixed_power<dim>(patch.n_subdivisions) *
                  (1 + GeometryInfo<dim>::vertices_per_cell);
              }
          }
        else
          {
            n_nodes += patch.data.n_cols();
            n_cells += 1;
            n_points_and_n_cells += patch.data.n_cols() + 1;
          }
      }
  }
 
  template <typename FlagsType>
  class StreamBase
  {
  public:
    /*
     * Constructor. Stores a reference to the output stream for immediate use.
     */
    StreamBase(std::ostream &stream, const FlagsType &flags)
      : selected_component(numbers::invalid_unsigned_int)
      , stream(stream)
      , flags(flags)
    {}
 
    template <int dim>
    void
    write_point(const unsigned int, const Point<dim> &)
    {
      Assert(false,
             ExcMessage("The derived class you are using needs to "
                        "reimplement this function if you want to call "
                        "it."));
    }
 
    void
    flush_points()
    {}
 
    template <int dim>
    void
    write_cell(const unsigned int /*index*/,
               const unsigned int /*start*/,
               const unsigned int /*x_offset*/,
               const unsigned int /*y_offset*/,
               const unsigned int /*z_offset*/)
    {
      Assert(false,
             ExcMessage("The derived class you are using needs to "
                        "reimplement this function if you want to call "
                        "it."));
    }
 
    void
    write_cell_single(const unsigned int   index,
                      const unsigned int   start,
                      const unsigned int   n_points,
                      const ReferenceCell &reference_cell)
    {
      (void)index;
      (void)start;
      (void)n_points;
      (void)reference_cell;
 
      Assert(false,
             ExcMessage("The derived class you are using needs to "
                        "reimplement this function if you want to call "
                        "it."));
    }
 
    void
    flush_cells()
    {}
 
    template <typename T>
    std::ostream &
    operator<<(const T &t)
    {
      stream << t;
      return stream;
    }
 
    unsigned int selected_component;
 
  protected:
    std::ostream &stream;
 
    const FlagsType flags;
  };
 
  class DXStream : public StreamBase<DataOutBase::DXFlags>
  {
  public:
    DXStream(std::ostream &stream, const DataOutBase::DXFlags &flags);
 
    template <int dim>
    void
    write_point(const unsigned int index, const Point<dim> &);
 
    template <int dim>
    void
    write_cell(const unsigned int index,
               const unsigned int start,
               const unsigned int x_offset,
               const unsigned int y_offset,
               const unsigned int z_offset);
 
    template <typename data>
    void
    write_dataset(const unsigned int index, const std::vector<data> &values);
  };
 
  class GmvStream : public StreamBase<DataOutBase::GmvFlags>
  {
  public:
    GmvStream(std::ostream &stream, const DataOutBase::GmvFlags &flags);
 
    template <int dim>
    void
    write_point(const unsigned int index, const Point<dim> &);
 
    template <int dim>
    void
    write_cell(const unsigned int index,
               const unsigned int start,
               const unsigned int x_offset,
               const unsigned int y_offset,
               const unsigned int z_offset);
  };
 
  class TecplotStream : public StreamBase<DataOutBase::TecplotFlags>
  {
  public:
    TecplotStream(std::ostream &stream, const DataOutBase::TecplotFlags &flags);
 
    template <int dim>
    void
    write_point(const unsigned int index, const Point<dim> &);
 
    template <int dim>
    void
    write_cell(const unsigned int index,
               const unsigned int start,
               const unsigned int x_offset,
               const unsigned int y_offset,
               const unsigned int z_offset);
  };
 
  class UcdStream : public StreamBase<DataOutBase::UcdFlags>
  {
  public:
    UcdStream(std::ostream &stream, const DataOutBase::UcdFlags &flags);
 
    template <int dim>
    void
    write_point(const unsigned int index, const Point<dim> &);
 
    template <int dim>
    void
    write_cell(const unsigned int index,
               const unsigned int start,
               const unsigned int x_offset,
               const unsigned int y_offset,
               const unsigned int z_offset);
 
    template <typename data>
    void
    write_dataset(const unsigned int index, const std::vector<data> &values);
  };
 
  class VtkStream : public StreamBase<DataOutBase::VtkFlags>
  {
  public:
    VtkStream(std::ostream &stream, const DataOutBase::VtkFlags &flags);
 
    template <int dim>
    void
    write_point(const unsigned int index, const Point<dim> &);
 
    template <int dim>
    void
    write_cell(const unsigned int index,
               const unsigned int start,
               const unsigned int x_offset,
               const unsigned int y_offset,
               const unsigned int z_offset);
 
    void
    write_cell_single(const unsigned int   index,
                      const unsigned int   start,
                      const unsigned int   n_points,
                      const ReferenceCell &reference_cell);
 
    template <int dim>
    void
    write_high_order_cell(const unsigned int           index,
                          const unsigned int           start,
                          const std::vector<unsigned> &connectivity);
  };
 
 
  class VtuStream : public StreamBase<DataOutBase::VtkFlags>
  {
  public:
    VtuStream(std::ostream &stream, const DataOutBase::VtkFlags &flags);
 
    template <int dim>
    void
    write_point(const unsigned int index, const Point<dim> &);
 
    void
    flush_points();
 
    template <int dim>
    void
    write_cell(const unsigned int index,
               const unsigned int start,
               const unsigned int x_offset,
               const unsigned int y_offset,
               const unsigned int z_offset);
 
    void
    write_cell_single(const unsigned int   index,
                      const unsigned int   start,
                      const unsigned int   n_points,
                      const ReferenceCell &reference_cell);
 
    template <int dim>
    void
    write_high_order_cell(const unsigned int           index,
                          const unsigned int           start,
                          const std::vector<unsigned> &connectivity);
 
    void
    flush_cells();
 
    template <typename T>
    std::ostream &
    operator<<(const T &);
 
    template <typename T>
    std::ostream &
    operator<<(const std::vector<T> &);
 
  private:
    std::vector<float>   vertices;
    std::vector<int32_t> cells;
  };
 
 
  //----------------------------------------------------------------------//
 
  DXStream::DXStream(std::ostream &out, const DataOutBase::DXFlags &f)
    : StreamBase<DataOutBase::DXFlags>(out, f)
  {}
 
 
  template <int dim>
  void
  DXStream::write_point(const unsigned int, const Point<dim> &p)
  {
    if (flags.coordinates_binary)
      {
        float data[dim];
        for (unsigned int d = 0; d < dim; ++d)
          data[d] = p(d);
        stream.write(reinterpret_cast<const char *>(data), dim * sizeof(*data));
      }
    else
      {
        for (unsigned int d = 0; d < dim; ++d)
          stream << p(d) << '\t';
        stream << '\n';
      }
  }
 
 
 
  // Separate these out to avoid an internal compiler error with intel 17
  namespace DataOutBaseImplementation
  {
    template <int dim>
    std::array<unsigned int, GeometryInfo<dim>::vertices_per_cell>
    set_node_numbers(const unsigned int /*start*/,
                     const unsigned int /*d1*/,
                     const unsigned int /*d2*/,
                     const unsigned int /*d3*/)
    {
      Assert(false, ExcInternalError());
      return {};
    }
 
 
 
    template <>
    std::array<unsigned int, GeometryInfo<1>::vertices_per_cell>
    set_node_numbers<1>(const unsigned int start,
                        const unsigned int d1,
                        const unsigned int /*d2*/,
                        const unsigned int /*d3*/)
 
    {
      std::array<unsigned int, GeometryInfo<1>::vertices_per_cell> nodes;
      nodes[0] = start;
      nodes[1] = start + d1;
      return nodes;
    }
 
 
 
    template <>
    std::array<unsigned int, GeometryInfo<2>::vertices_per_cell>
    set_node_numbers<2>(const unsigned int start,
                        const unsigned int d1,
                        const unsigned int d2,
                        const unsigned int /*d3*/)
 
    {
      std::array<unsigned int, GeometryInfo<2>::vertices_per_cell> nodes;
      nodes[0] = start;
      nodes[1] = start + d1;
      nodes[2] = start + d2;
      nodes[3] = start + d2 + d1;
      return nodes;
    }
 
 
 
    template <>
    std::array<unsigned int, GeometryInfo<3>::vertices_per_cell>
    set_node_numbers<3>(const unsigned int start,
                        const unsigned int d1,
                        const unsigned int d2,
                        const unsigned int d3)
    {
      std::array<unsigned int, GeometryInfo<3>::vertices_per_cell> nodes;
      nodes[0] = start;
      nodes[1] = start + d1;
      nodes[2] = start + d2;
      nodes[3] = start + d2 + d1;
      nodes[4] = start + d3;
      nodes[5] = start + d3 + d1;
      nodes[6] = start + d3 + d2;
      nodes[7] = start + d3 + d2 + d1;
      return nodes;
    }
  } // namespace DataOutBaseImplementation
 
 
 
  template <int dim>
  void
  DXStream::write_cell(unsigned int,
                       unsigned int start,
                       unsigned int d1,
                       unsigned int d2,
                       unsigned int d3)
  {
    const auto nodes =
      DataOutBaseImplementation::set_node_numbers<dim>(start, d1, d2, d3);
 
    if (flags.int_binary)
      {
        std::array<unsigned int, GeometryInfo<dim>::vertices_per_cell> temp;
        for (unsigned int i = 0; i < nodes.size(); ++i)
          temp[i] = nodes[GeometryInfo<dim>::dx_to_deal[i]];
        stream.write(reinterpret_cast<const char *>(temp.data()),
                     temp.size() * sizeof(temp[0]));
      }
    else
      {
        for (unsigned int i = 0; i < nodes.size() - 1; ++i)
          stream << nodes[GeometryInfo<dim>::dx_to_deal[i]] << '\t';
        stream << nodes[GeometryInfo<dim>::dx_to_deal[nodes.size() - 1]]
               << '\n';
      }
  }
 
 
 
  template <typename data>
  inline void
  DXStream::write_dataset(const unsigned int, const std::vector<data> &values)
  {
    if (flags.data_binary)
      {
        stream.write(reinterpret_cast<const char *>(values.data()),
                     values.size() * sizeof(data));
      }
    else
      {
        for (unsigned int i = 0; i < values.size(); ++i)
          stream << '\t' << values[i];
        stream << '\n';
      }
  }
 
 
 
  //----------------------------------------------------------------------//
 
  GmvStream::GmvStream(std::ostream &out, const DataOutBase::GmvFlags &f)
    : StreamBase<DataOutBase::GmvFlags>(out, f)
  {}
 
 
  template <int dim>
  void
  GmvStream::write_point(const unsigned int, const Point<dim> &p)
  {
    Assert(selected_component != numbers::invalid_unsigned_int,
           ExcNotInitialized());
    stream << p(selected_component) << ' ';
  }
 
 
 
  template <int dim>
  void
  GmvStream::write_cell(unsigned int,
                        unsigned int s,
                        unsigned int d1,
                        unsigned int d2,
                        unsigned int d3)
  {
    // Vertices are numbered starting with one.
    const unsigned int start = s + 1;
    stream << gmv_cell_type[dim] << '\n';
 
    stream << start;
    if (dim >= 1)
      {
        stream << '\t' << start + d1;
        if (dim >= 2)
          {
            stream << '\t' << start + d2 + d1 << '\t' << start + d2;
            if (dim >= 3)
              {
                stream << '\t' << start + d3 << '\t' << start + d3 + d1 << '\t'
                       << start + d3 + d2 + d1 << '\t' << start + d3 + d2;
              }
          }
      }
    stream << '\n';
  }
 
 
 
  TecplotStream::TecplotStream(std::ostream &                   out,
                               const DataOutBase::TecplotFlags &f)
    : StreamBase<DataOutBase::TecplotFlags>(out, f)
  {}
 
 
  template <int dim>
  void
  TecplotStream::write_point(const unsigned int, const Point<dim> &p)
  {
    Assert(selected_component != numbers::invalid_unsigned_int,
           ExcNotInitialized());
    stream << p(selected_component) << '\n';
  }
 
 
 
  template <int dim>
  void
  TecplotStream::write_cell(unsigned int,
                            unsigned int s,
                            unsigned int d1,
                            unsigned int d2,
                            unsigned int d3)
  {
    const unsigned int start = s + 1;
 
    stream << start;
    if (dim >= 1)
      {
        stream << '\t' << start + d1;
        if (dim >= 2)
          {
            stream << '\t' << start + d2 + d1 << '\t' << start + d2;
            if (dim >= 3)
              {
                stream << '\t' << start + d3 << '\t' << start + d3 + d1 << '\t'
                       << start + d3 + d2 + d1 << '\t' << start + d3 + d2;
              }
          }
      }
    stream << '\n';
  }
 
 
 
  UcdStream::UcdStream(std::ostream &out, const DataOutBase::UcdFlags &f)
    : StreamBase<DataOutBase::UcdFlags>(out, f)
  {}
 
 
  template <int dim>
  void
  UcdStream::write_point(const unsigned int index, const Point<dim> &p)
  {
    stream << index + 1 << "   ";
    // write out coordinates
    for (unsigned int i = 0; i < dim; ++i)
      stream << p(i) << ' ';
    // fill with zeroes
    for (unsigned int i = dim; i < 3; ++i)
      stream << "0 ";
    stream << '\n';
  }
 
 
 
  template <int dim>
  void
  UcdStream::write_cell(unsigned int index,
                        unsigned int start,
                        unsigned int d1,
                        unsigned int d2,
                        unsigned int d3)
  {
    const auto nodes =
      DataOutBaseImplementation::set_node_numbers<dim>(start, d1, d2, d3);
 
    // Write out all cells and remember that all indices must be shifted by one.
    stream << index + 1 << "\t0 " << ucd_cell_type[dim];
    for (unsigned int i = 0; i < nodes.size(); ++i)
      stream << '\t' << nodes[GeometryInfo<dim>::ucd_to_deal[i]] + 1;
    stream << '\n';
  }
 
 
 
  template <typename data>
  inline void
  UcdStream::write_dataset(const unsigned int       index,
                           const std::vector<data> &values)
  {
    stream << index + 1;
    for (unsigned int i = 0; i < values.size(); ++i)
      stream << '\t' << values[i];
    stream << '\n';
  }
 
 
 
  //----------------------------------------------------------------------//
 
  VtkStream::VtkStream(std::ostream &out, const DataOutBase::VtkFlags &f)
    : StreamBase<DataOutBase::VtkFlags>(out, f)
  {}
 
 
  template <int dim>
  void
  VtkStream::write_point(const unsigned int, const Point<dim> &p)
  {
    // write out coordinates
    stream << p;
    // fill with zeroes
    for (unsigned int i = dim; i < 3; ++i)
      stream << " 0";
    stream << '\n';
  }
 
 
 
  template <int dim>
  void
  VtkStream::write_cell(unsigned int,
                        unsigned int start,
                        unsigned int d1,
                        unsigned int d2,
                        unsigned int d3)
  {
    stream << GeometryInfo<dim>::vertices_per_cell << '\t' << start;
 
    if (dim >= 1)
      stream << '\t' << start + d1;
    {
      if (dim >= 2)
        {
          stream << '\t' << start + d2 + d1 << '\t' << start + d2;
          if (dim >= 3)
            {
              stream << '\t' << start + d3 << '\t' << start + d3 + d1 << '\t'
                     << start + d3 + d2 + d1 << '\t' << start + d3 + d2;
            }
        }
    }
    stream << '\n';
  }
 
  void
  VtkStream::write_cell_single(const unsigned int   index,
                               const unsigned int   start,
                               const unsigned int   n_points,
                               const ReferenceCell &reference_cell)
  {
    (void)index;
 
    static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
 
    stream << '\t' << n_points;
    for (unsigned int i = 0; i < n_points; ++i)
      stream << '\t'
             << start +
                  (reference_cell == ReferenceCells::Pyramid ? table[i] : i);
    stream << '\n';
  }
 
  template <int dim>
  void
  VtkStream::write_high_order_cell(const unsigned int,
                                   const unsigned int           start,
                                   const std::vector<unsigned> &connectivity)
  {
    stream << connectivity.size();
    for (const auto &c : connectivity)
      stream << '\t' << start + c;
    stream << '\n';
  }
 
  VtuStream::VtuStream(std::ostream &out, const DataOutBase::VtkFlags &f)
    : StreamBase<DataOutBase::VtkFlags>(out, f)
  {}
 
 
  template <int dim>
  void
  VtuStream::write_point(const unsigned int, const Point<dim> &p)
  {
#if !defined(DEAL_II_WITH_ZLIB)
    // write out coordinates
    stream << p;
    // fill with zeroes
    for (unsigned int i = dim; i < 3; ++i)
      stream << " 0";
    stream << '\n';
#else
    // if we want to compress, then first collect all the data in an array
    for (unsigned int i = 0; i < dim; ++i)
      vertices.push_back(p[i]);
    for (unsigned int i = dim; i < 3; ++i)
      vertices.push_back(0);
#endif
  }
 
 
  void
  VtuStream::flush_points()
  {
#ifdef DEAL_II_WITH_ZLIB
    // compress the data we have in memory and write them to the stream. then
    // release the data
    *this << vertices << '\n';
    vertices.clear();
#endif
  }
 
 
  template <int dim>
  void
  VtuStream::write_cell(unsigned int,
                        unsigned int start,
                        unsigned int d1,
                        unsigned int d2,
                        unsigned int d3)
  {
#if !defined(DEAL_II_WITH_ZLIB)
    stream << start;
    if (dim >= 1)
      {
        stream << '\t' << start + d1;
        if (dim >= 2)
          {
            stream << '\t' << start + d2 + d1 << '\t' << start + d2;
            if (dim >= 3)
              {
                stream << '\t' << start + d3 << '\t' << start + d3 + d1 << '\t'
                       << start + d3 + d2 + d1 << '\t' << start + d3 + d2;
              }
          }
      }
    stream << '\n';
#else
    cells.push_back(start);
    if (dim >= 1)
      {
        cells.push_back(start + d1);
        if (dim >= 2)
          {
            cells.push_back(start + d2 + d1);
            cells.push_back(start + d2);
            if (dim >= 3)
              {
                cells.push_back(start + d3);
                cells.push_back(start + d3 + d1);
                cells.push_back(start + d3 + d2 + d1);
                cells.push_back(start + d3 + d2);
              }
          }
      }
#endif
  }
 
  void
  VtuStream::write_cell_single(const unsigned int   index,
                               const unsigned int   start,
                               const unsigned int   n_points,
                               const ReferenceCell &reference_cell)
  {
    (void)index;
 
    static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
 
#if !defined(DEAL_II_WITH_ZLIB)
    for (unsigned int i = 0; i < n_points; ++i)
      stream << '\t'
             << start +
                  (reference_cell == ReferenceCells::Pyramid ? table[i] : i);
    stream << '\n';
#else
    for (unsigned int i = 0; i < n_points; ++i)
      cells.push_back(
        start + (reference_cell == ReferenceCells::Pyramid ? table[i] : i));
#endif
  }
 
  template <int dim>
  void
  VtuStream::write_high_order_cell(const unsigned int,
                                   const unsigned int           start,
                                   const std::vector<unsigned> &connectivity)
  {
#if !defined(DEAL_II_WITH_ZLIB)
    for (const auto &c : connectivity)
      stream << '\t' << start + c;
    stream << '\n';
#else
    for (const auto &c : connectivity)
      cells.push_back(start + c);
#endif
  }
 
  void
  VtuStream::flush_cells()
  {
#ifdef DEAL_II_WITH_ZLIB
    // compress the data we have in memory and write them to the stream. then
    // release the data
    *this << cells << '\n';
    cells.clear();
#endif
  }
 
 
  template <typename T>
  std::ostream &
  VtuStream::operator<<(const std::vector<T> &data)
  {
#ifdef DEAL_II_WITH_ZLIB
    // compress the data we have in memory and write them to the stream. then
    // release the data
    write_compressed_block(data, flags, stream);
#else
    for (unsigned int i = 0; i < data.size(); ++i)
      stream << data[i] << ' ';
#endif
 
    return stream;
  }
} // namespace
 
 
 
namespace DataOutBase
{
  const unsigned int Deal_II_IntermediateFlags::format_version = 4;
 
 
  template <int dim, int spacedim>
  const unsigned int Patch<dim, spacedim>::space_dim;
 
 
  template <int dim, int spacedim>
  const unsigned int Patch<dim, spacedim>::no_neighbor;
 
 
  template <int dim, int spacedim>
  Patch<dim, spacedim>::Patch()
    : patch_index(no_neighbor)
    , n_subdivisions(1)
    , points_are_available(false)
    , reference_cell(ReferenceCells::Invalid)
  // all the other data has a constructor of its own, except for the "neighbors"
  // field, which we set to invalid values.
  {
    for (unsigned int i : GeometryInfo<dim>::face_indices())
      neighbors[i] = no_neighbor;
 
    AssertIndexRange(dim, spacedim + 1);
    Assert(spacedim <= 3, ExcNotImplemented());
  }
 
 
 
  template <int dim, int spacedim>
  bool
  Patch<dim, spacedim>::operator==(const Patch &patch) const
  {
    if (reference_cell != patch.reference_cell)
      return false;
 
    // TODO: make tolerance relative
    const double epsilon = 3e-16;
    for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
      if (vertices[i].distance(patch.vertices[i]) > epsilon)
        return false;
 
    for (unsigned int i : GeometryInfo<dim>::face_indices())
      if (neighbors[i] != patch.neighbors[i])
        return false;
 
    if (patch_index != patch.patch_index)
      return false;
 
    if (n_subdivisions != patch.n_subdivisions)
      return false;
 
    if (points_are_available != patch.points_are_available)
      return false;
 
    if (data.n_rows() != patch.data.n_rows())
      return false;
 
    if (data.n_cols() != patch.data.n_cols())
      return false;
 
    for (unsigned int i = 0; i < data.n_rows(); ++i)
      for (unsigned int j = 0; j < data.n_cols(); ++j)
        if (data[i][j] != patch.data[i][j])
          return false;
 
    return true;
  }
 
 
 
  template <int dim, int spacedim>
  std::size_t
  Patch<dim, spacedim>::memory_consumption() const
  {
    return (sizeof(vertices) / sizeof(vertices[0]) *
              MemoryConsumption::memory_consumption(vertices[0]) +
            sizeof(neighbors) / sizeof(neighbors[0]) *
              MemoryConsumption::memory_consumption(neighbors[0]) +
            MemoryConsumption::memory_consumption(patch_index) +
            MemoryConsumption::memory_consumption(n_subdivisions) +
            MemoryConsumption::memory_consumption(data) +
            MemoryConsumption::memory_consumption(points_are_available) +
            sizeof(reference_cell));
  }
 
 
 
  template <int dim, int spacedim>
  void
  Patch<dim, spacedim>::swap(Patch<dim, spacedim> &other_patch)
  {
    std::swap(vertices, other_patch.vertices);
    std::swap(neighbors, other_patch.neighbors);
    std::swap(patch_index, other_patch.patch_index);
    std::swap(n_subdivisions, other_patch.n_subdivisions);
    data.swap(other_patch.data);
    std::swap(points_are_available, other_patch.points_are_available);
    std::swap(reference_cell, other_patch.reference_cell);
  }
 
 
 
  template <int spacedim>
  const unsigned int Patch<0, spacedim>::space_dim;
 
 
  template <int spacedim>
  const unsigned int Patch<0, spacedim>::no_neighbor;
 
 
  template <int spacedim>
  unsigned int Patch<0, spacedim>::neighbors[1] = {
    Patch<0, spacedim>::no_neighbor};
 
  template <int spacedim>
  const unsigned int Patch<0, spacedim>::n_subdivisions = 1;
 
  template <int spacedim>
  const ReferenceCell Patch<0, spacedim>::reference_cell =
    ReferenceCells::Vertex;
 
  template <int spacedim>
  Patch<0, spacedim>::Patch()
    : patch_index(no_neighbor)
    , points_are_available(false)
  {
    Assert(spacedim <= 3, ExcNotImplemented());
  }
 
 
 
  template <int spacedim>
  bool
  Patch<0, spacedim>::operator==(const Patch &patch) const
  {
    const unsigned int dim = 0;
 
    // TODO: make tolerance relative
    const double epsilon = 3e-16;
    for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
      if (vertices[i].distance(patch.vertices[i]) > epsilon)
        return false;
 
    if (patch_index != patch.patch_index)
      return false;
 
    if (points_are_available != patch.points_are_available)
      return false;
 
    if (data.n_rows() != patch.data.n_rows())
      return false;
 
    if (data.n_cols() != patch.data.n_cols())
      return false;
 
    for (unsigned int i = 0; i < data.n_rows(); ++i)
      for (unsigned int j = 0; j < data.n_cols(); ++j)
        if (data[i][j] != patch.data[i][j])
          return false;
 
    return true;
  }
 
 
 
  template <int spacedim>
  std::size_t
  Patch<0, spacedim>::memory_consumption() const
  {
    return (sizeof(vertices) / sizeof(vertices[0]) *
              MemoryConsumption::memory_consumption(vertices[0]) +
            MemoryConsumption::memory_consumption(data) +
            MemoryConsumption::memory_consumption(points_are_available));
  }
 
 
 
  template <int spacedim>
  void
  Patch<0, spacedim>::swap(Patch<0, spacedim> &other_patch)
  {
    std::swap(vertices, other_patch.vertices);
    std::swap(patch_index, other_patch.patch_index);
    data.swap(other_patch.data);
    std::swap(points_are_available, other_patch.points_are_available);
  }
 
 
 
  UcdFlags::UcdFlags(const bool write_preamble)
    : write_preamble(write_preamble)
  {}
 
 
 
  GnuplotFlags::GnuplotFlags()
  {
    space_dimension_labels.emplace_back("x");
    space_dimension_labels.emplace_back("y");
    space_dimension_labels.emplace_back("z");
  }
 
 
 
  GnuplotFlags::GnuplotFlags(const std::vector<std::string> &labels)
    : space_dimension_labels(labels)
  {}
 
 
 
  std::size_t
  GnuplotFlags::memory_consumption() const
  {
    return MemoryConsumption::memory_consumption(space_dimension_labels);
  }
 
 
 
  PovrayFlags::PovrayFlags(const bool smooth,
                           const bool bicubic_patch,
                           const bool external_data)
    : smooth(smooth)
    , bicubic_patch(bicubic_patch)
    , external_data(external_data)
  {}
 
 
  DataOutFilterFlags::DataOutFilterFlags(const bool filter_duplicate_vertices,
                                         const bool xdmf_hdf5_output)
    : filter_duplicate_vertices(filter_duplicate_vertices)
    , xdmf_hdf5_output(xdmf_hdf5_output)
  {}
 
 
  void
  DataOutFilterFlags::declare_parameters(ParameterHandler &prm)
  {
    prm.declare_entry(
      "Filter duplicate vertices",
      "false",
      Patterns::Bool(),
      "Whether to remove duplicate vertex values. deal.II duplicates "
      "vertices once for each adjacent cell so that it can output "
      "discontinuous quantities for which there may be more than one "
      "value for each vertex position. Setting this flag to "
      "'true' will merge all of these values by selecting a "
      "random one and outputting this as 'the' value for the vertex. "
      "As long as the data to be output corresponds to continuous "
      "fields, merging vertices has no effect. On the other hand, "
      "if the data to be output corresponds to discontinuous fields "
      "(either because you are using a discontinuous finite element, "
      "or because you are using a DataPostprocessor that yields "
      "discontinuous data, or because the data to be output has been "
      "produced by entirely different means), then the data in the "
      "output file no longer faithfully represents the underlying data "
      "because the discontinuous field has been replaced by a "
      "continuous one. Note also that the filtering can not occur "
      "on processor boundaries. Thus, a filtered discontinuous field "
      "looks like a continuous field inside of a subdomain, "
      "but like a discontinuous field at the subdomain boundary."
      "\n\n"
      "In any case, filtering results in drastically smaller output "
      "files (smaller by about a factor of 2^dim).");
    prm.declare_entry(
      "XDMF HDF5 output",
      "false",
      Patterns::Bool(),
      "Whether the data will be used in an XDMF/HDF5 combination.");
  }
 
 
 
  void
  DataOutFilterFlags::parse_parameters(const ParameterHandler &prm)
  {
    filter_duplicate_vertices = prm.get_bool("Filter duplicate vertices");
    xdmf_hdf5_output          = prm.get_bool("XDMF HDF5 output");
  }
 
 
 
  DXFlags::DXFlags(const bool write_neighbors,
                   const bool int_binary,
                   const bool coordinates_binary,
                   const bool data_binary)
    : write_neighbors(write_neighbors)
    , int_binary(int_binary)
    , coordinates_binary(coordinates_binary)
    , data_binary(data_binary)
    , data_double(false)
  {}
 
 
  void
  DXFlags::declare_parameters(ParameterHandler &prm)
  {
    prm.declare_entry("Write neighbors",
                      "true",
                      Patterns::Bool(),
                      "A boolean field indicating whether neighborship "
                      "information between cells is to be written to the "
                      "OpenDX output file");
    prm.declare_entry("Integer format",
                      "ascii",
                      Patterns::Selection("ascii|32|64"),
                      "Output format of integer numbers, which is "
                      "either a text representation (ascii) or binary integer "
                      "values of 32 or 64 bits length");
    prm.declare_entry("Coordinates format",
                      "ascii",
                      Patterns::Selection("ascii|32|64"),
                      "Output format of vertex coordinates, which is "
                      "either a text representation (ascii) or binary "
                      "floating point values of 32 or 64 bits length");
    prm.declare_entry("Data format",
                      "ascii",
                      Patterns::Selection("ascii|32|64"),
                      "Output format of data values, which is "
                      "either a text representation (ascii) or binary "
                      "floating point values of 32 or 64 bits length");
  }
 
 
 
  void
  DXFlags::parse_parameters(const ParameterHandler &prm)
  {
    write_neighbors = prm.get_bool("Write neighbors");
    // TODO:[GK] Read the new  parameters
  }
 
 
 
  void
  UcdFlags::declare_parameters(ParameterHandler &prm)
  {
    prm.declare_entry("Write preamble",
                      "true",
                      Patterns::Bool(),
                      "A flag indicating whether a comment should be "
                      "written to the beginning of the output file "
                      "indicating date and time of creation as well "
                      "as the creating program");
  }
 
 
 
  void
  UcdFlags::parse_parameters(const ParameterHandler &prm)
  {
    write_preamble = prm.get_bool("Write preamble");
  }
 
 
 
  SvgFlags::SvgFlags(const unsigned int height_vector,
                     const int          azimuth_angle,
                     const int          polar_angle,
                     const unsigned int line_thickness,
                     const bool         margin,
                     const bool         draw_colorbar)
    : height(4000)
    , width(0)
    , height_vector(height_vector)
    , azimuth_angle(azimuth_angle)
    , polar_angle(polar_angle)
    , line_thickness(line_thickness)
    , margin(margin)
    , draw_colorbar(draw_colorbar)
  {}
 
 
 
  void
  PovrayFlags::declare_parameters(ParameterHandler &prm)
  {
    prm.declare_entry("Use smooth triangles",
                      "false",
                      Patterns::Bool(),
                      "A flag indicating whether POVRAY should use smoothed "
                      "triangles instead of the usual ones");
    prm.declare_entry("Use bicubic patches",
                      "false",
                      Patterns::Bool(),
                      "Whether POVRAY should use bicubic patches");
    prm.declare_entry("Include external file",
                      "true",
                      Patterns::Bool(),
                      "Whether camera and lighting information should "
                      "be put into an external file \"data.inc\" or into "
                      "the POVRAY input file");
  }
 
 
 
  void
  PovrayFlags::parse_parameters(const ParameterHandler &prm)
  {
    smooth        = prm.get_bool("Use smooth triangles");
    bicubic_patch = prm.get_bool("Use bicubic patches");
    external_data = prm.get_bool("Include external file");
  }
 
 
 
  EpsFlags::EpsFlags(const unsigned int  height_vector,
                     const unsigned int  color_vector,
                     const SizeType      size_type,
                     const unsigned int  size,
                     const double        line_width,
                     const double        azimut_angle,
                     const double        turn_angle,
                     const double        z_scaling,
                     const bool          draw_mesh,
                     const bool          draw_cells,
                     const bool          shade_cells,
                     const ColorFunction color_function)
    : height_vector(height_vector)
    , color_vector(color_vector)
    , size_type(size_type)
    , size(size)
    , line_width(line_width)
    , azimut_angle(azimut_angle)
    , turn_angle(turn_angle)
    , z_scaling(z_scaling)
    , draw_mesh(draw_mesh)
    , draw_cells(draw_cells)
    , shade_cells(shade_cells)
    , color_function(color_function)
  {}
 
 
 
  EpsFlags::RgbValues
  EpsFlags::default_color_function(const double x,
                                   const double xmin,
                                   const double xmax)
  {
    RgbValues rgb_values = {0, 0, 0};
 
    // A difficult color scale:
    //     xmin          = black  (1)
    // 3/4*xmin+1/4*xmax = blue   (2)
    // 1/2*xmin+1/2*xmax = green  (3)
    // 1/4*xmin+3/4*xmax = red    (4)
    //              xmax = white  (5)
    // Makes the following color functions:
    //
    // red      green    blue
    //       __
    //      /      /\  /  /\    /
    // ____/    __/  \/  /  \__/
 
    //     { 0                                (1) - (3)
    // r = { ( 4*x-2*xmin+2*xmax)/(xmax-xmin) (3) - (4)
    //     { 1                                (4) - (5)
    //
    //     { 0                                (1) - (2)
    // g = { ( 4*x-3*xmin-  xmax)/(xmax-xmin) (2) - (3)
    //     { (-4*x+  xmin+3*xmax)/(xmax-xmin) (3) - (4)
    //     { ( 4*x-  xmin-3*xmax)/(xmax-xmin) (4) - (5)
    //
    //     { ( 4*x-4*xmin       )/(xmax-xmin) (1) - (2)
    // b = { (-4*x+2*xmin+2*xmax)/(xmax-xmin) (2) - (3)
    //     { 0                                (3) - (4)
    //     { ( 4*x-  xmin-3*xmax)/(xmax-xmin) (4) - (5)
 
    double sum    = xmax + xmin;
    double sum13  = xmin + 3 * xmax;
    double sum22  = 2 * xmin + 2 * xmax;
    double sum31  = 3 * xmin + xmax;
    double dif    = xmax - xmin;
    double rezdif = 1.0 / dif;
 
    int where;
 
    if (x < (sum31) / 4)
      where = 0;
    else if (x < (sum22) / 4)
      where = 1;
    else if (x < (sum13) / 4)
      where = 2;
    else
      where = 3;
 
    if (dif != 0)
      {
        switch (where)
          {
            case 0:
              rgb_values.red   = 0;
              rgb_values.green = 0;
              rgb_values.blue  = (x - xmin) * 4. * rezdif;
              break;
            case 1:
              rgb_values.red   = 0;
              rgb_values.green = (4 * x - 3 * xmin - xmax) * rezdif;
              rgb_values.blue  = (sum22 - 4. * x) * rezdif;
              break;
            case 2:
              rgb_values.red   = (4 * x - 2 * sum) * rezdif;
              rgb_values.green = (xmin + 3 * xmax - 4 * x) * rezdif;
              rgb_values.blue  = 0;
              break;
            case 3:
              rgb_values.red   = 1;
              rgb_values.green = (4 * x - xmin - 3 * xmax) * rezdif;
              rgb_values.blue  = (4. * x - sum13) * rezdif;
              break;
            default:
              break;
          }
      }
    else // White
      rgb_values.red = rgb_values.green = rgb_values.blue = 1;
 
    return rgb_values;
  }
 
 
 
  EpsFlags::RgbValues
  EpsFlags::grey_scale_color_function(const double x,
                                      const double xmin,
                                      const double xmax)
  {
    EpsFlags::RgbValues rgb_values;
    rgb_values.red = rgb_values.blue = rgb_values.green =
      (x - xmin) / (xmax - xmin);
    return rgb_values;
  }
 
 
 
  EpsFlags::RgbValues
  EpsFlags::reverse_grey_scale_color_function(const double x,
                                              const double xmin,
                                              const double xmax)
  {
    EpsFlags::RgbValues rgb_values;
    rgb_values.red = rgb_values.blue = rgb_values.green =
      1 - (x - xmin) / (xmax - xmin);
    return rgb_values;
  }
 
 
 
  void
  EpsFlags::declare_parameters(ParameterHandler &prm)
  {
    prm.declare_entry("Index of vector for height",
                      "0",
                      Patterns::Integer(),
                      "Number of the input vector that is to be used to "
                      "generate height information");
    prm.declare_entry("Index of vector for color",
                      "0",
                      Patterns::Integer(),
                      "Number of the input vector that is to be used to "
                      "generate color information");
    prm.declare_entry("Scale to width or height",
                      "width",
                      Patterns::Selection("width|height"),
                      "Whether width or height should be scaled to match "
                      "the given size");
    prm.declare_entry("Size (width or height) in eps units",
                      "300",
                      Patterns::Integer(),
                      "The size (width or height) to which the eps output "
                      "file is to be scaled");
    prm.declare_entry("Line widths in eps units",
                      "0.5",
                      Patterns::Double(),
                      "The width in which the postscript renderer is to "
                      "plot lines");
    prm.declare_entry("Azimut angle",
                      "60",
                      Patterns::Double(0, 180),
                      "Angle of the viewing position against the vertical "
                      "axis");
    prm.declare_entry("Turn angle",
                      "30",
                      Patterns::Double(0, 360),
                      "Angle of the viewing direction against the y-axis");
    prm.declare_entry("Scaling for z-axis",
                      "1",
                      Patterns::Double(),
                      "Scaling for the z-direction relative to the scaling "
                      "used in x- and y-directions");
    prm.declare_entry("Draw mesh lines",
                      "true",
                      Patterns::Bool(),
                      "Whether the mesh lines, or only the surface should be "
                      "drawn");
    prm.declare_entry("Fill interior of cells",
                      "true",
                      Patterns::Bool(),
                      "Whether only the mesh lines, or also the interior of "
                      "cells should be plotted. If this flag is false, then "
                      "one can see through the mesh");
    prm.declare_entry("Color shading of interior of cells",
                      "true",
                      Patterns::Bool(),
                      "Whether the interior of cells shall be shaded");
    prm.declare_entry("Color function",
                      "default",
                      Patterns::Selection(
                        "default|grey scale|reverse grey scale"),
                      "Name of a color function used to colorize mesh lines "
                      "and/or cell interiors");
  }
 
 
 
  void
  EpsFlags::parse_parameters(const ParameterHandler &prm)
  {
    height_vector = prm.get_integer("Index of vector for height");
    color_vector  = prm.get_integer("Index of vector for color");
    if (prm.get("Scale to width or height") == "width")
      size_type = width;
    else
      size_type = height;
    size         = prm.get_integer("Size (width or height) in eps units");
    line_width   = prm.get_double("Line widths in eps units");
    azimut_angle = prm.get_double("Azimut angle");
    turn_angle   = prm.get_double("Turn angle");
    z_scaling    = prm.get_double("Scaling for z-axis");
    draw_mesh    = prm.get_bool("Draw mesh lines");
    draw_cells   = prm.get_bool("Fill interior of cells");
    shade_cells  = prm.get_bool("Color shading of interior of cells");
    if (prm.get("Color function") == "default")
      color_function = &default_color_function;
    else if (prm.get("Color function") == "grey scale")
      color_function = &grey_scale_color_function;
    else if (prm.get("Color function") == "reverse grey scale")
      color_function = &reverse_grey_scale_color_function;
    else
      // we shouldn't get here, since the parameter object should already have
      // checked that the given value is valid
      Assert(false, ExcInternalError());
  }
 
 
 
  TecplotFlags::TecplotFlags(const char *zone_name, const double solution_time)
    : zone_name(zone_name)
    , solution_time(solution_time)
  {}
 
 
 
  std::size_t
  TecplotFlags::memory_consumption() const
  {
    return sizeof(*this) + MemoryConsumption::memory_consumption(zone_name);
  }
 
 
 
  VtkFlags::VtkFlags(const double                         time,
                     const unsigned int                   cycle,
                     const bool                           print_date_and_time,
                     const VtkFlags::ZlibCompressionLevel compression_level,
                     const bool write_higher_order_cells,
                     const std::map<std::string, std::string> &physical_units)
    : time(time)
    , cycle(cycle)
    , print_date_and_time(print_date_and_time)
    , compression_level(compression_level)
    , write_higher_order_cells(write_higher_order_cells)
    , physical_units(physical_units)
  {}
 
 
 
  OutputFormat
  parse_output_format(const std::string &format_name)
  {
    if (format_name == "none")
      return none;
 
    if (format_name == "dx")
      return dx;
 
    if (format_name == "ucd")
      return ucd;
 
    if (format_name == "gnuplot")
      return gnuplot;
 
    if (format_name == "povray")
      return povray;
 
    if (format_name == "eps")
      return eps;
 
    if (format_name == "gmv")
      return gmv;
 
    if (format_name == "tecplot")
      return tecplot;
 
    if (format_name == "tecplot_binary")
      return tecplot_binary;
 
    if (format_name == "vtk")
      return vtk;
 
    if (format_name == "vtu")
      return vtu;
 
    if (format_name == "deal.II intermediate")
      return deal_II_intermediate;
 
    if (format_name == "hdf5")
      return hdf5;
 
    AssertThrow(false,
                ExcMessage("The given file format name is not recognized: <" +
                           format_name + ">"));
 
    // return something invalid
    return OutputFormat(-1);
  }
 
 
 
  std::string
  get_output_format_names()
  {
    return "none|dx|ucd|gnuplot|povray|eps|gmv|tecplot|tecplot_binary|vtk|vtu|hdf5|svg|deal.II intermediate";
  }
 
 
 
  std::string
  default_suffix(const OutputFormat output_format)
  {
    switch (output_format)
      {
        case none:
          return "";
        case dx:
          return ".dx";
        case ucd:
          return ".inp";
        case gnuplot:
          return ".gnuplot";
        case povray:
          return ".pov";
        case eps:
          return ".eps";
        case gmv:
          return ".gmv";
        case tecplot:
          return ".dat";
        case tecplot_binary:
          return ".plt";
        case vtk:
          return ".vtk";
        case vtu:
          return ".vtu";
        case deal_II_intermediate:
          return ".d2";
        case hdf5:
          return ".h5";
        case svg:
          return ".svg";
        default:
          Assert(false, ExcNotImplemented());
          return "";
      }
  }
 
 
  //----------------------------------------------------------------------//
 
  template <int dim, int spacedim, typename StreamType>
  void
  write_nodes(const std::vector<Patch<dim, spacedim>> &patches, StreamType &out)
  {
    Assert(dim <= 3, ExcNotImplemented());
    unsigned int count = 0;
 
    static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
 
    for (const auto &patch : patches)
      {
        // special treatment of non-hypercube cells
        if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
          {
            for (unsigned int point_no = 0; point_no < patch.data.n_cols();
                 ++point_no)
              out.write_point(count++,
                              get_node_location(patch,
                                                (patch.reference_cell ==
                                                     ReferenceCells::Pyramid ?
                                                   table[point_no] :
                                                   point_no)));
          }
        else
          {
            const unsigned int n_subdivisions = patch.n_subdivisions;
            const unsigned int n              = n_subdivisions + 1;
 
            switch (dim)
              {
                case 0:
                  out.write_point(count++,
                                  get_equispaced_location(patch,
                                                          {},
                                                          n_subdivisions));
                  break;
                case 1:
                  for (unsigned int i1 = 0; i1 < n; ++i1)
                    out.write_point(count++,
                                    get_equispaced_location(patch,
                                                            {i1},
                                                            n_subdivisions));
                  break;
                case 2:
                  for (unsigned int i2 = 0; i2 < n; ++i2)
                    for (unsigned int i1 = 0; i1 < n; ++i1)
                      out.write_point(count++,
                                      get_equispaced_location(patch,
                                                              {i1, i2},
                                                              n_subdivisions));
                  break;
                case 3:
                  for (unsigned int i3 = 0; i3 < n; ++i3)
                    for (unsigned int i2 = 0; i2 < n; ++i2)
                      for (unsigned int i1 = 0; i1 < n; ++i1)
                        out.write_point(count++,
                                        get_equispaced_location(
                                          patch, {i1, i2, i3}, n_subdivisions));
                  break;
 
                default:
                  Assert(false, ExcInternalError());
              }
          }
      }
    out.flush_points();
  }
 
  template <int dim, int spacedim, typename StreamType>
  void
  write_cells(const std::vector<Patch<dim, spacedim>> &patches, StreamType &out)
  {
    Assert(dim <= 3, ExcNotImplemented());
    unsigned int count                 = 0;
    unsigned int first_vertex_of_patch = 0;
    for (const auto &patch : patches)
      {
        // special treatment of simplices since they are not subdivided
        if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
          {
            out.write_cell_single(count++,
                                  first_vertex_of_patch,
                                  patch.data.n_cols(),
                                  patch.reference_cell);
            first_vertex_of_patch += patch.data.n_cols();
          }
        else
          {
            const unsigned int n_subdivisions = patch.n_subdivisions;
            const unsigned int n              = n_subdivisions + 1;
            // Length of loops in all dimensions
            const unsigned int n1 = (dim > 0) ? n_subdivisions : 1;
            const unsigned int n2 = (dim > 1) ? n_subdivisions : 1;
            const unsigned int n3 = (dim > 2) ? n_subdivisions : 1;
            // Offsets of outer loops
            const unsigned int d1 = 1;
            const unsigned int d2 = n;
            const unsigned int d3 = n * n;
            for (unsigned int i3 = 0; i3 < n3; ++i3)
              for (unsigned int i2 = 0; i2 < n2; ++i2)
                for (unsigned int i1 = 0; i1 < n1; ++i1)
                  {
                    const unsigned int offset =
                      first_vertex_of_patch + i3 * d3 + i2 * d2 + i1 * d1;
                    // First write line in x direction
                    out.template write_cell<dim>(count++, offset, d1, d2, d3);
                  }
            // finally update the number of the first vertex of this patch
            first_vertex_of_patch +=
              Utilities::fixed_power<dim>(n_subdivisions + 1);
          }
      }
 
    out.flush_cells();
  }
 
  template <int dim, int spacedim, typename StreamType>
  void
  write_high_order_cells(const std::vector<Patch<dim, spacedim>> &patches,
                         StreamType &                             out)
  {
    Assert(dim <= 3 && dim > 1, ExcNotImplemented());
    unsigned int first_vertex_of_patch = 0;
    unsigned int count                 = 0;
    // Array to hold all the node numbers of a cell
    std::vector<unsigned> connectivity;
    // Array to hold cell order in each dimension
    std::array<unsigned, dim> cell_order;
 
    for (const auto &patch : patches)
      {
        if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
          {
            connectivity.resize(patch.data.n_cols());
 
            for (unsigned int i = 0; i < patch.data.n_cols(); ++i)
              connectivity[i] = i;
 
            out.template write_high_order_cell<dim>(count++,
                                                    first_vertex_of_patch,
                                                    connectivity);
 
            first_vertex_of_patch += patch.data.n_cols();
          }
        else
          {
            const unsigned int n_subdivisions = patch.n_subdivisions;
            const unsigned int n              = n_subdivisions + 1;
 
            cell_order.fill(n_subdivisions);
            connectivity.resize(Utilities::fixed_power<dim>(n));
 
            // Length of loops in all dimensons
            const unsigned int n1 = (dim > 0) ? n_subdivisions : 0;
            const unsigned int n2 = (dim > 1) ? n_subdivisions : 0;
            const unsigned int n3 = (dim > 2) ? n_subdivisions : 0;
            // Offsets of outer loops
            const unsigned int d1 = 1;
            const unsigned int d2 = n;
            const unsigned int d3 = n * n;
            for (unsigned int i3 = 0; i3 <= n3; ++i3)
              for (unsigned int i2 = 0; i2 <= n2; ++i2)
                for (unsigned int i1 = 0; i1 <= n1; ++i1)
                  {
                    const unsigned int local_index =
                      i3 * d3 + i2 * d2 + i1 * d1;
                    const unsigned int connectivity_index =
                      vtk_point_index_from_ijk(i1, i2, i3, cell_order);
                    connectivity[connectivity_index] = local_index;
                  }
 
            out.template write_high_order_cell<dim>(count++,
                                                    first_vertex_of_patch,
                                                    connectivity);
 
            // finally update the number of the first vertex of this patch
            first_vertex_of_patch += Utilities::fixed_power<dim>(n);
          }
      }
 
    out.flush_cells();
  }
 
 
  template <int dim, int spacedim, class StreamType>
  void
  write_data(const std::vector<Patch<dim, spacedim>> &patches,
             unsigned int                             n_data_sets,
             const bool                               double_precision,
             StreamType &                             out)
  {
    Assert(dim <= 3, ExcNotImplemented());
    unsigned int count = 0;
 
    for (const auto &patch : patches)
      {
        const unsigned int n_subdivisions = patch.n_subdivisions;
        const unsigned int n              = n_subdivisions + 1;
        // Length of loops in all dimensions
        Assert((patch.data.n_rows() == n_data_sets &&
                !patch.points_are_available) ||
                 (patch.data.n_rows() == n_data_sets + spacedim &&
                  patch.points_are_available),
               ExcDimensionMismatch(patch.points_are_available ?
                                      (n_data_sets + spacedim) :
                                      n_data_sets,
                                    patch.data.n_rows()));
        Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(n),
               ExcInvalidDatasetSize(patch.data.n_cols(), n));
 
        std::vector<float>  floats(n_data_sets);
        std::vector<double> doubles(n_data_sets);
 
        // Data is already in lexicographic ordering
        for (unsigned int i = 0; i < Utilities::fixed_power<dim>(n);
             ++i, ++count)
          if (double_precision)
            {
              for (unsigned int data_set = 0; data_set < n_data_sets;
                   ++data_set)
                doubles[data_set] = patch.data(data_set, i);
              out.write_dataset(count, doubles);
            }
          else
            {
              for (unsigned int data_set = 0; data_set < n_data_sets;
                   ++data_set)
                floats[data_set] = patch.data(data_set, i);
              out.write_dataset(count, floats);
            }
      }
  }
 
 
 
  namespace
  {
    Point<2>
    svg_project_point(Point<3> point,
                      Point<3> camera_position,
                      Point<3> camera_direction,
                      Point<3> camera_horizontal,
                      float    camera_focus)
    {
      Point<3> camera_vertical;
      camera_vertical[0] = camera_horizontal[1] * camera_direction[2] -
                           camera_horizontal[2] * camera_direction[1];
      camera_vertical[1] = camera_horizontal[2] * camera_direction[0] -
                           camera_horizontal[0] * camera_direction[2];
      camera_vertical[2] = camera_horizontal[0] * camera_direction[1] -
                           camera_horizontal[1] * camera_direction[0];
 
      float phi;
      phi = camera_focus;
      phi /= (point[0] - camera_position[0]) * camera_direction[0] +
             (point[1] - camera_position[1]) * camera_direction[1] +
             (point[2] - camera_position[2]) * camera_direction[2];
 
      Point<3> projection;
      projection[0] =
        camera_position[0] + phi * (point[0] - camera_position[0]);
      projection[1] =
        camera_position[1] + phi * (point[1] - camera_position[1]);
      projection[2] =
        camera_position[2] + phi * (point[2] - camera_position[2]);
 
      Point<2> projection_decomposition;
      projection_decomposition[0] = (projection[0] - camera_position[0] -
                                     camera_focus * camera_direction[0]) *
                                    camera_horizontal[0];
      projection_decomposition[0] += (projection[1] - camera_position[1] -
                                      camera_focus * camera_direction[1]) *
                                     camera_horizontal[1];
      projection_decomposition[0] += (projection[2] - camera_position[2] -
                                      camera_focus * camera_direction[2]) *
                                     camera_horizontal[2];
 
      projection_decomposition[1] = (projection[0] - camera_position[0] -
                                     camera_focus * camera_direction[0]) *
                                    camera_vertical[0];
      projection_decomposition[1] += (projection[1] - camera_position[1] -
                                      camera_focus * camera_direction[1]) *
                                     camera_vertical[1];
      projection_decomposition[1] += (projection[2] - camera_position[2] -
                                      camera_focus * camera_direction[2]) *
                                     camera_vertical[2];
 
      return projection_decomposition;
    }
 
 
    Point<6>
    svg_get_gradient_parameters(Point<3> points[])
    {
      Point<3> v_min, v_max, v_inter;
 
      // Use the Bubblesort algorithm to sort the points with respect to the
      // third coordinate
      for (int i = 0; i < 2; ++i)
        {
          for (int j = 0; j < 2 - i; ++j)
            {
              if (points[j][2] > points[j + 1][2])
                {
                  Point<3> temp = points[j];
                  points[j]     = points[j + 1];
                  points[j + 1] = temp;
                }
            }
        }
 
      // save the related three-dimensional vectors v_min, v_inter, and v_max
      v_min   = points[0];
      v_inter = points[1];
      v_max   = points[2];
 
      Point<2> A[2];
      Point<2> b, gradient;
 
      // determine the plane offset c
      A[0][0] = v_max[0] - v_min[0];
      A[0][1] = v_inter[0] - v_min[0];
      A[1][0] = v_max[1] - v_min[1];
      A[1][1] = v_inter[1] - v_min[1];
 
      b[0] = -v_min[0];
      b[1] = -v_min[1];
 
      double x, sum;
      bool   col_change = false;
 
      if (A[0][0] == 0)
        {
          col_change = true;
 
          A[0][0] = A[0][1];
          A[0][1] = 0;
 
          double temp = A[1][0];
          A[1][0]     = A[1][1];
          A[1][1]     = temp;
        }
 
      for (unsigned int k = 0; k < 1; ++k)
        {
          for (unsigned int i = k + 1; i < 2; ++i)
            {
              x = A[i][k] / A[k][k];
 
              for (unsigned int j = k + 1; j < 2; ++j)
                A[i][j] = A[i][j] - A[k][j] * x;
 
              b[i] = b[i] - b[k] * x;
            }
        }
 
      b[1] = b[1] / A[1][1];
 
      for (int i = 0; i >= 0; i--)
        {
          sum = b[i];
 
          for (unsigned int j = i + 1; j < 2; ++j)
            sum = sum - A[i][j] * b[j];
 
          b[i] = sum / A[i][i];
        }
 
      if (col_change)
        {
          double temp = b[0];
          b[0]        = b[1];
          b[1]        = temp;
        }
 
      double c = b[0] * (v_max[2] - v_min[2]) + b[1] * (v_inter[2] - v_min[2]) +
                 v_min[2];
 
      // Determine the first entry of the gradient (phi, cf. documentation)
      A[0][0] = v_max[0] - v_min[0];
      A[0][1] = v_inter[0] - v_min[0];
      A[1][0] = v_max[1] - v_min[1];
      A[1][1] = v_inter[1] - v_min[1];
 
      b[0] = 1.0 - v_min[0];
      b[1] = -v_min[1];
 
      col_change = false;
 
      if (A[0][0] == 0)
        {
          col_change = true;
 
          A[0][0] = A[0][1];
          A[0][1] = 0;
 
          double temp = A[1][0];
          A[1][0]     = A[1][1];
          A[1][1]     = temp;
        }
 
      for (unsigned int k = 0; k < 1; ++k)
        {
          for (unsigned int i = k + 1; i < 2; ++i)
            {
              x = A[i][k] / A[k][k];
 
              for (unsigned int j = k + 1; j < 2; ++j)
                A[i][j] = A[i][j] - A[k][j] * x;
 
              b[i] = b[i] - b[k] * x;
            }
        }
 
      b[1] = b[1] / A[1][1];
 
      for (int i = 0; i >= 0; i--)
        {
          sum = b[i];
 
          for (unsigned int j = i + 1; j < 2; ++j)
            sum = sum - A[i][j] * b[j];
 
          b[i] = sum / A[i][i];
        }
 
      if (col_change)
        {
          double temp = b[0];
          b[0]        = b[1];
          b[1]        = temp;
        }
 
      gradient[0] = b[0] * (v_max[2] - v_min[2]) +
                    b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
 
      // determine the second entry of the gradient
      A[0][0] = v_max[0] - v_min[0];
      A[0][1] = v_inter[0] - v_min[0];
      A[1][0] = v_max[1] - v_min[1];
      A[1][1] = v_inter[1] - v_min[1];
 
      b[0] = -v_min[0];
      b[1] = 1.0 - v_min[1];
 
      col_change = false;
 
      if (A[0][0] == 0)
        {
          col_change = true;
 
          A[0][0] = A[0][1];
          A[0][1] = 0;
 
          double temp = A[1][0];
          A[1][0]     = A[1][1];
          A[1][1]     = temp;
        }
 
      for (unsigned int k = 0; k < 1; ++k)
        {
          for (unsigned int i = k + 1; i < 2; ++i)
            {
              x = A[i][k] / A[k][k];
 
              for (unsigned int j = k + 1; j < 2; ++j)
                A[i][j] = A[i][j] - A[k][j] * x;
 
              b[i] = b[i] - b[k] * x;
            }
        }
 
      b[1] = b[1] / A[1][1];
 
      for (int i = 0; i >= 0; i--)
        {
          sum = b[i];
 
          for (unsigned int j = i + 1; j < 2; ++j)
            sum = sum - A[i][j] * b[j];
 
          b[i] = sum / A[i][i];
        }
 
      if (col_change)
        {
          double temp = b[0];
          b[0]        = b[1];
          b[1]        = temp;
        }
 
      gradient[1] = b[0] * (v_max[2] - v_min[2]) +
                    b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
 
      // normalize the gradient
      gradient /= gradient.norm();
 
      const double lambda = -gradient[0] * (v_min[0] - v_max[0]) -
                            gradient[1] * (v_min[1] - v_max[1]);
 
      Point<6> gradient_parameters;
 
      gradient_parameters[0] = v_min[0];
      gradient_parameters[1] = v_min[1];
 
      gradient_parameters[2] = v_min[0] + lambda * gradient[0];
      gradient_parameters[3] = v_min[1] + lambda * gradient[1];
 
      gradient_parameters[4] = v_min[2];
      gradient_parameters[5] = v_max[2];
 
      return gradient_parameters;
    }
  } // namespace
 
 
 
  template <int dim, int spacedim>
  void
  write_ucd(
    const std::vector<Patch<dim, spacedim>> &patches,
    const std::vector<std::string> &         data_names,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>> &,
    const UcdFlags &flags,
    std::ostream &  out)
  {
    // Note that while in theory dim==0 should be implemented, this is not
    // tested, therefore currently not allowed.
    AssertThrow(dim > 0, ExcNotImplemented());
 
    AssertThrow(out.fail() == false, ExcIO());
 
#ifndef DEAL_II_WITH_MPI
    // verify that there are indeed patches to be written out. most of the
    // times, people just forget to call build_patches when there are no
    // patches, so a warning is in order. that said, the assertion is disabled
    // if we support MPI since then it can happen that on the coarsest mesh, a
    // processor simply has no cells it actually owns, and in that case it is
    // legit if there are no patches
    Assert(patches.size() > 0, ExcNoPatches());
#else
    if (patches.size() == 0)
      return;
#endif
 
    const unsigned int n_data_sets = data_names.size();
 
    UcdStream ucd_out(out, flags);
 
    // first count the number of cells and cells for later use
    unsigned int n_nodes;
    unsigned int n_cells;
    compute_sizes<dim, spacedim>(patches, n_nodes, n_cells);
    //---------------------
    // preamble
    if (flags.write_preamble)
      {
        out
          << "# This file was generated by the deal.II library." << '\n'
          << "# Date =  " << Utilities::System::get_date() << '\n'
          << "# Time =  " << Utilities::System::get_time() << '\n'
          << "#" << '\n'
          << "# For a description of the UCD format see the AVS Developer's guide."
          << '\n'
          << "#" << '\n';
      }
 
    // start with ucd data
    out << n_nodes << ' ' << n_cells << ' ' << n_data_sets << ' ' << 0
        << ' ' // no cell data at present
        << 0   // no model data
        << '\n';
 
    write_nodes(patches, ucd_out);
    out << '\n';
 
    write_cells(patches, ucd_out);
    out << '\n';
 
    //---------------------------
    // now write data
    if (n_data_sets != 0)
      {
        out << n_data_sets << "    "; // number of vectors
        for (unsigned int i = 0; i < n_data_sets; ++i)
          out << 1 << ' '; // number of components;
        // only 1 supported presently
        out << '\n';
 
        for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
          out << data_names[data_set]
              << ",dimensionless" // no units supported at present
              << '\n';
 
        write_data(patches, n_data_sets, true, ucd_out);
      }
    // make sure everything now gets to disk
    out.flush();
 
    // assert the stream is still ok
    AssertThrow(out.fail() == false, ExcIO());
  }
 
 
  template <int dim, int spacedim>
  void
  write_dx(
    const std::vector<Patch<dim, spacedim>> &patches,
    const std::vector<std::string> &         data_names,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>> &,
    const DXFlags &flags,
    std::ostream & out)
  {
    // Point output is currently not implemented.
    AssertThrow(dim > 0, ExcNotImplemented());
 
    AssertThrow(out.fail() == false, ExcIO());
 
#ifndef DEAL_II_WITH_MPI
    // verify that there are indeed patches to be written out. most of the
    // times, people just forget to call build_patches when there are no
    // patches, so a warning is in order. that said, the assertion is disabled
    // if we support MPI since then it can happen that on the coarsest mesh, a
    // processor simply has no cells it actually owns, and in that case it is
    // legit if there are no patches
    Assert(patches.size() > 0, ExcNoPatches());
#else
    if (patches.size() == 0)
      return;
#endif
    // Stream with special features for dx output
    DXStream dx_out(out, flags);
 
    // Variable counting the offset of binary data.
    unsigned int offset = 0;
 
    const unsigned int n_data_sets = data_names.size();
 
    // first count the number of cells and cells for later use
    unsigned int n_nodes;
    unsigned int n_cells;
    compute_sizes<dim, spacedim>(patches, n_nodes, n_cells);
    // start with vertices order is lexicographical, x varying fastest
    out << "object \"vertices\" class array type float rank 1 shape "
        << spacedim << " items " << n_nodes;
 
    if (flags.coordinates_binary)
      {
        out << " lsb ieee data 0" << '\n';
        offset += n_nodes * spacedim * sizeof(float);
      }
    else
      {
        out << " data follows" << '\n';
        write_nodes(patches, dx_out);
      }
 
    //-----------------------------
    // first write the coordinates of all vertices
 
    //---------------------------------------
    // write cells
    out << "object \"cells\" class array type int rank 1 shape "
        << GeometryInfo<dim>::vertices_per_cell << " items " << n_cells;
 
    if (flags.int_binary)
      {
        out << " lsb binary data " << offset << '\n';
        offset += n_cells * sizeof(int);
      }
    else
      {
        out << " data follows" << '\n';
        write_cells(patches, dx_out);
        out << '\n';
      }
 
 
    out << "attribute \"element type\" string \"";
    if (dim == 1)
      out << "lines";
    if (dim == 2)
      out << "quads";
    if (dim == 3)
      out << "cubes";
    out << "\"" << '\n' << "attribute \"ref\" string \"positions\"" << '\n';
 
    // TODO:[GK] Patches must be of same size!
    //---------------------------
    // write neighbor information
    if (flags.write_neighbors)
      {
        out << "object \"neighbors\" class array type int rank 1 shape "
            << GeometryInfo<dim>::faces_per_cell << " items " << n_cells
            << " data follows";
 
        for (const auto &patch : patches)
          {
            const unsigned int n               = patch.n_subdivisions;
            const unsigned int n1              = (dim > 0) ? n : 1;
            const unsigned int n2              = (dim > 1) ? n : 1;
            const unsigned int n3              = (dim > 2) ? n : 1;
            unsigned int       cells_per_patch = Utilities::fixed_power<dim>(n);
            unsigned int       dx              = 1;
            unsigned int       dy              = n;
            unsigned int       dz              = n * n;
 
            const unsigned int patch_start =
              patch.patch_index * cells_per_patch;
 
            for (unsigned int i3 = 0; i3 < n3; ++i3)
              for (unsigned int i2 = 0; i2 < n2; ++i2)
                for (unsigned int i1 = 0; i1 < n1; ++i1)
                  {
                    const unsigned int nx = i1 * dx;
                    const unsigned int ny = i2 * dy;
                    const unsigned int nz = i3 * dz;
 
                    // There are no neighbors for dim==0. Note that this case is
                    // caught by the AssertThrow at the beginning of this
                    // function anyway. This condition avoids compiler warnings.
                    if (dim < 1)
                      continue;
 
                    out << '\n';
                    // Direction -x Last cell in row of other patch
                    if (i1 == 0)
                      {
                        const unsigned int nn = patch.neighbors[0];
                        out << '\t';
                        if (nn != patch.no_neighbor)
                          out
                            << (nn * cells_per_patch + ny + nz + dx * (n - 1));
                        else
                          out << "-1";
                      }
                    else
                      {
                        out << '\t' << patch_start + nx - dx + ny + nz;
                      }
                    // Direction +x First cell in row of other patch
                    if (i1 == n - 1)
                      {
                        const unsigned int nn = patch.neighbors[1];
                        out << '\t';
                        if (nn != patch.no_neighbor)
                          out << (nn * cells_per_patch + ny + nz);
                        else
                          out << "-1";
                      }
                    else
                      {
                        out << '\t' << patch_start + nx + dx + ny + nz;
                      }
                    if (dim < 2)
                      continue;
                    // Direction -y
                    if (i2 == 0)
                      {
                        const unsigned int nn = patch.neighbors[2];
                        out << '\t';
                        if (nn != patch.no_neighbor)
                          out
                            << (nn * cells_per_patch + nx + nz + dy * (n - 1));
                        else
                          out << "-1";
                      }
                    else
                      {
                        out << '\t' << patch_start + nx + ny - dy + nz;
                      }
                    // Direction +y
                    if (i2 == n - 1)
                      {
                        const unsigned int nn = patch.neighbors[3];
                        out << '\t';
                        if (nn != patch.no_neighbor)
                          out << (nn * cells_per_patch + nx + nz);
                        else
                          out << "-1";
                      }
                    else
                      {
                        out << '\t' << patch_start + nx + ny + dy + nz;
                      }
                    if (dim < 3)
                      continue;
 
                    // Direction -z
                    if (i3 == 0)
                      {
                        const unsigned int nn = patch.neighbors[4];
                        out << '\t';
                        if (nn != patch.no_neighbor)
                          out
                            << (nn * cells_per_patch + nx + ny + dz * (n - 1));
                        else
                          out << "-1";
                      }
                    else
                      {
                        out << '\t' << patch_start + nx + ny + nz - dz;
                      }
                    // Direction +z
                    if (i3 == n - 1)
                      {
                        const unsigned int nn = patch.neighbors[5];
                        out << '\t';
                        if (nn != patch.no_neighbor)
                          out << (nn * cells_per_patch + nx + ny);
                        else
                          out << "-1";
                      }
                    else
                      {
                        out << '\t' << patch_start + nx + ny + nz + dz;
                      }
                  }
            out << '\n';
          }
      }
    //---------------------------
    // now write data
    if (n_data_sets != 0)
      {
        out << "object \"data\" class array type float rank 1 shape "
            << n_data_sets << " items " << n_nodes;
 
        if (flags.data_binary)
          {
            out << " lsb ieee data " << offset << '\n';
            offset += n_data_sets * n_nodes *
                      ((flags.data_double) ? sizeof(double) : sizeof(float));
          }
        else
          {
            out << " data follows" << '\n';
            write_data(patches, n_data_sets, flags.data_double, dx_out);
          }
 
        // loop over all patches
        out << "attribute \"dep\" string \"positions\"" << '\n';
      }
    else
      {
        out << "object \"data\" class constantarray type float rank 0 items "
            << n_nodes << " data follows" << '\n'
            << '0' << '\n';
      }
 
    // no model data
 
    out << "object \"deal data\" class field" << '\n'
        << "component \"positions\" value \"vertices\"" << '\n'
        << "component \"connections\" value \"cells\"" << '\n'
        << "component \"data\" value \"data\"" << '\n';
 
    if (flags.write_neighbors)
      out << "component \"neighbors\" value \"neighbors\"" << '\n';
 
    {
      out << "attribute \"created\" string \"" << Utilities::System::get_date()
          << ' ' << Utilities::System::get_time() << '"' << '\n';
    }
 
    out << "end" << '\n';
    // Write all binary data now
    if (flags.coordinates_binary)
      write_nodes(patches, dx_out);
    if (flags.int_binary)
      write_cells(patches, dx_out);
    if (flags.data_binary)
      write_data(patches, n_data_sets, flags.data_double, dx_out);
 
    // make sure everything now gets to disk
    out.flush();
 
    // assert the stream is still ok
    AssertThrow(out.fail() == false, ExcIO());
  }
 
 
 
  template <int dim, int spacedim>
  void
  write_gnuplot(
    const std::vector<Patch<dim, spacedim>> &patches,
    const std::vector<std::string> &         data_names,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>> &,
    const GnuplotFlags &flags,
    std::ostream &      out)
  {
    AssertThrow(out.fail() == false, ExcIO());
 
#ifndef DEAL_II_WITH_MPI
    // verify that there are indeed patches to be written out. most
    // of the times, people just forget to call build_patches when there
    // are no patches, so a warning is in order. that said, the
    // assertion is disabled if we support MPI since then it can
    // happen that on the coarsest mesh, a processor simply has no
    // cells it actually owns, and in that case it is legit if there
    // are no patches
    Assert(patches.size() > 0, ExcNoPatches());
#else
    if (patches.size() == 0)
      return;
#endif
 
    const unsigned int n_data_sets = data_names.size();
 
    // write preamble
    {
      out << "# This file was generated by the deal.II library." << '\n'
          << "# Date =  " << Utilities::System::get_date() << '\n'
          << "# Time =  " << Utilities::System::get_time() << '\n'
          << "#" << '\n'
          << "# For a description of the GNUPLOT format see the GNUPLOT manual."
          << '\n'
          << "#" << '\n'
          << "# ";
 
      AssertThrow(spacedim <= flags.space_dimension_labels.size(),
                  GnuplotFlags::ExcNotEnoughSpaceDimensionLabels());
      for (unsigned int spacedim_n = 0; spacedim_n < spacedim; ++spacedim_n)
        {
          out << '<' << flags.space_dimension_labels.at(spacedim_n) << "> ";
        }
 
      for (const auto &data_name : data_names)
        out << '<' << data_name << "> ";
      out << '\n';
    }
 
 
    // loop over all patches
    for (const auto &patch : patches)
      {
        const unsigned int n_subdivisions         = patch.n_subdivisions;
        const unsigned int n_points_per_direction = n_subdivisions + 1;
 
        Assert((patch.data.n_rows() == n_data_sets &&
                !patch.points_are_available) ||
                 (patch.data.n_rows() == n_data_sets + spacedim &&
                  patch.points_are_available),
               ExcDimensionMismatch(patch.points_are_available ?
                                      (n_data_sets + spacedim) :
                                      n_data_sets,
                                    patch.data.n_rows()));
 
        auto output_point_data =
          [&out, &patch, n_data_sets](const unsigned int point_index) mutable {
            for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
              out << patch.data(data_set, point_index) << ' ';
          };
 
        switch (dim)
          {
            case 0:
              {
                Assert(patch.reference_cell == ReferenceCells::Vertex,
                       ExcInternalError());
                Assert(patch.data.n_cols() == 1,
                       ExcInvalidDatasetSize(patch.data.n_cols(),
                                             n_subdivisions + 1));
 
 
                // compute coordinates for this patch point
                out << get_equispaced_location(patch, {}, n_subdivisions)
                    << ' ';
                output_point_data(0);
                out << '\n';
                out << '\n';
                break;
              }
 
            case 1:
              {
                Assert(patch.reference_cell == ReferenceCells::Line,
                       ExcInternalError());
                Assert(patch.data.n_cols() ==
                         Utilities::fixed_power<dim>(n_points_per_direction),
                       ExcInvalidDatasetSize(patch.data.n_cols(),
                                             n_subdivisions + 1));
 
                for (unsigned int i1 = 0; i1 < n_points_per_direction; ++i1)
                  {
                    // compute coordinates for this patch point
                    out << get_equispaced_location(patch, {i1}, n_subdivisions)
                        << ' ';
 
                    output_point_data(i1);
                    out << '\n';
                  }
                // end of patch
                out << '\n';
                out << '\n';
                break;
              }
 
            case 2:
              {
                if (patch.reference_cell == ReferenceCells::Quadrilateral)
                  {
                    Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(
                                                    n_points_per_direction),
                           ExcInvalidDatasetSize(patch.data.n_cols(),
                                                 n_subdivisions + 1));
 
                    for (unsigned int i2 = 0; i2 < n_points_per_direction; ++i2)
                      {
                        for (unsigned int i1 = 0; i1 < n_points_per_direction;
                             ++i1)
                          {
                            // compute coordinates for this patch point
                            out << get_equispaced_location(patch,
                                                           {i1, i2},
                                                           n_subdivisions)
                                << ' ';
 
                            output_point_data(i1 + i2 * n_points_per_direction);
                            out << '\n';
                          }
                        // end of row in patch
                        out << '\n';
                      }
                  }
                else if (patch.reference_cell == ReferenceCells::Triangle)
                  {
                    Assert(n_subdivisions == 1, ExcNotImplemented());
 
                    Assert(patch.data.n_cols() == 3, ExcInternalError());
 
                    // Gnuplot can only plot surfaces if each facet of the
                    // surface is a bilinear patch, or a subdivided bilinear
                    // patch with equally many points along each row of the
                    // subdivision. This is what the code above for
                    // quadrilaterals does. We emulate this by repeating the
                    // third point of a triangle twice so that there are two
                    // points for that row as well -- i.e., we write a 2x2
                    // bilinear patch where two of the points are collapsed onto
                    // one vertex.
                    //
                    // This also matches the example here:
                    // https://stackoverflow.com/questions/42784369/drawing-triangular-mesh-using-gnuplot
                    out << get_node_location(patch, 0) << ' ';
                    output_point_data(0);
                    out << '\n';
 
                    out << get_node_location(patch, 1) << ' ';
                    output_point_data(1);
                    out << '\n';
                    out << '\n'; // end of one row of points
 
                    out << get_node_location(patch, 2) << ' ';
                    output_point_data(2);
                    out << '\n';
 
                    out << get_node_location(patch, 2) << ' ';
                    output_point_data(2);
                    out << '\n';
                    out << '\n'; // end of the second row of points
                    out << '\n'; // end of the entire patch
                  }
                else
                  // There aren't any other reference cells in 2d than the
                  // quadrilateral and the triangle. So whatever we got here
                  // can't be any good
                  Assert(false, ExcInternalError());
                // end of patch
                out << '\n';
 
                break;
              }
 
            case 3:
              {
                if (patch.reference_cell == ReferenceCells::Hexahedron)
                  {
                    Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(
                                                    n_points_per_direction),
                           ExcInvalidDatasetSize(patch.data.n_cols(),
                                                 n_subdivisions + 1));
 
                    // for all grid points: draw lines into all positive
                    // coordinate directions if there is another grid point
                    // there
                    for (unsigned int i3 = 0; i3 < n_points_per_direction; ++i3)
                      for (unsigned int i2 = 0; i2 < n_points_per_direction;
                           ++i2)
                        for (unsigned int i1 = 0; i1 < n_points_per_direction;
                             ++i1)
                          {
                            // compute coordinates for this patch point
                            const Point<spacedim> this_point =
                              get_equispaced_location(patch,
                                                      {i1, i2, i3},
                                                      n_subdivisions);
                            // line into positive x-direction if possible
                            if (i1 < n_subdivisions)
                              {
                                // write point here and its data
                                out << this_point << ' ';
                                output_point_data(i1 +
                                                  i2 * n_points_per_direction +
                                                  i3 * n_points_per_direction *
                                                    n_points_per_direction);
                                out << '\n';
 
                                // write point there and its data
                                out << get_equispaced_location(patch,
                                                               {i1 + 1, i2, i3},
                                                               n_subdivisions)
                                    << ' ';
 
                                output_point_data((i1 + 1) +
                                                  i2 * n_points_per_direction +
                                                  i3 * n_points_per_direction *
                                                    n_points_per_direction);
                                out << '\n';
 
                                // end of line
                                out << '\n' << '\n';
                              }
 
                            // line into positive y-direction if possible
                            if (i2 < n_subdivisions)
                              {
                                // write point here and its data
                                out << this_point << ' ';
                                output_point_data(i1 +
                                                  i2 * n_points_per_direction +
                                                  i3 * n_points_per_direction *
                                                    n_points_per_direction);
                                out << '\n';
 
                                // write point there and its data
                                out << get_equispaced_location(patch,
                                                               {i1, i2 + 1, i3},
                                                               n_subdivisions)
                                    << ' ';
 
                                output_point_data(
                                  i1 + (i2 + 1) * n_points_per_direction +
                                  i3 * n_points_per_direction *
                                    n_points_per_direction);
                                out << '\n';
 
                                // end of line
                                out << '\n' << '\n';
                              }
 
                            // line into positive z-direction if possible
                            if (i3 < n_subdivisions)
                              {
                                // write point here and its data
                                out << this_point << ' ';
                                output_point_data(i1 +
                                                  i2 * n_points_per_direction +
                                                  i3 * n_points_per_direction *
                                                    n_points_per_direction);
                                out << '\n';
 
                                // write point there and its data
                                out << get_equispaced_location(patch,
                                                               {i1, i2, i3 + 1},
                                                               n_subdivisions)
                                    << ' ';
 
                                output_point_data(
                                  i1 + i2 * n_points_per_direction +
                                  (i3 + 1) * n_points_per_direction *
                                    n_points_per_direction);
                                out << '\n';
                                // end of line
                                out << '\n' << '\n';
                              }
                          }
                  }
                else if (patch.reference_cell == ReferenceCells::Tetrahedron)
                  {
                    Assert(n_subdivisions == 1, ExcNotImplemented());
 
                    // Draw the tetrahedron as a collection of two lines.
                    for (const unsigned int v : {0, 1, 2, 0, 3, 2})
                      {
                        out << get_node_location(patch, v) << ' ';
                        output_point_data(v);
                        out << '\n';
                      }
                    out << '\n'; // end of first line
 
                    for (const unsigned int v : {3, 1})
                      {
                        out << get_node_location(patch, v) << ' ';
                        output_point_data(v);
                        out << '\n';
                      }
                    out << '\n'; // end of second line
                  }
                else if (patch.reference_cell == ReferenceCells::Pyramid)
                  {
                    Assert(n_subdivisions == 1, ExcNotImplemented());
 
                    // Draw the pyramid as a collection of two lines.
                    for (const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
                      {
                        out << get_node_location(patch, v) << ' ';
                        output_point_data(v);
                        out << '\n';
                      }
                    out << '\n'; // end of first line
 
                    for (const unsigned int v : {2, 4, 3})
                      {
                        out << get_node_location(patch, v) << ' ';
                        output_point_data(v);
                        out << '\n';
                      }
                    out << '\n'; // end of second line
                  }
                else if (patch.reference_cell == ReferenceCells::Wedge)
                  {
                    Assert(n_subdivisions == 1, ExcNotImplemented());
 
                    // Draw the wedge as a collection of three
                    // lines. The first one wraps around the base,
                    // goes up to the top, and wraps around that. The
                    // second and third are just individual lines
                    // going from base to top.
                    for (const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
                      {
                        out << get_node_location(patch, v) << ' ';
                        output_point_data(v);
                        out << '\n';
                      }
                    out << '\n'; // end of first line
 
                    for (const unsigned int v : {1, 4})
                      {
                        out << get_node_location(patch, v) << ' ';
                        output_point_data(v);
                        out << '\n';
                      }
                    out << '\n'; // end of second line
 
                    for (const unsigned int v : {2, 5})
                      {
                        out << get_node_location(patch, v) << ' ';
                        output_point_data(v);
                        out << '\n';
                      }
                    out << '\n'; // end of second line
                  }
                else
                  // No other reference cells are currently implemented
                  Assert(false, ExcNotImplemented());
 
                break;
              }
 
            default:
              Assert(false, ExcNotImplemented());
          }
      }
    // make sure everything now gets to disk
    out.flush();
 
    AssertThrow(out.fail() == false, ExcIO());
  }
 
 
  namespace
  {
    template <int dim, int spacedim>
    void
    do_write_povray(const std::vector<Patch<dim, spacedim>> &,
                    const std::vector<std::string> &,
                    const PovrayFlags &,
                    std::ostream &)
    {
      Assert(false,
             ExcMessage("Writing files in POVRAY format is only supported "
                        "for two-dimensional meshes."));
    }
 
 
 
    void
    do_write_povray(const std::vector<Patch<2, 2>> &patches,
                    const std::vector<std::string> &data_names,
                    const PovrayFlags &             flags,
                    std::ostream &                  out)
    {
      AssertThrow(out.fail() == false, ExcIO());
 
#ifndef DEAL_II_WITH_MPI
      // verify that there are indeed patches to be written out. most
      // of the times, people just forget to call build_patches when there
      // are no patches, so a warning is in order. that said, the
      // assertion is disabled if we support MPI since then it can
      // happen that on the coarsest mesh, a processor simply has no cells it
      // actually owns, and in that case it is legit if there are no patches
      Assert(patches.size() > 0, ExcNoPatches());
#else
      if (patches.size() == 0)
        return;
#endif
      constexpr int dim = 2;
      (void)dim;
      constexpr int spacedim = 2;
 
      const unsigned int n_data_sets = data_names.size();
      (void)n_data_sets;
 
      // write preamble
      {
        out
          << "/* This file was generated by the deal.II library." << '\n'
          << "   Date =  " << Utilities::System::get_date() << '\n'
          << "   Time =  " << Utilities::System::get_time() << '\n'
          << '\n'
          << "   For a description of the POVRAY format see the POVRAY manual."
          << '\n'
          << "*/ " << '\n';
 
        // include files
        out << "#include \"colors.inc\" " << '\n'
            << "#include \"textures.inc\" " << '\n';
 
 
        // use external include file for textures, camera and light
        if (flags.external_data)
          out << "#include \"data.inc\" " << '\n';
        else // all definitions in data file
          {
            // camera
            out << '\n'
                << '\n'
                << "camera {" << '\n'
                << "  location <1,4,-7>" << '\n'
                << "  look_at <0,0,0>" << '\n'
                << "  angle 30" << '\n'
                << "}" << '\n';
 
            // light
            out << '\n'
                << "light_source {" << '\n'
                << "  <1,4,-7>" << '\n'
                << "  color Grey" << '\n'
                << "}" << '\n';
            out << '\n'
                << "light_source {" << '\n'
                << "  <0,20,0>" << '\n'
                << "  color White" << '\n'
                << "}" << '\n';
          }
      }
 
      // max. and min. height of solution
      Assert(patches.size() > 0, ExcNoPatches());
      double hmin = patches[0].data(0, 0);
      double hmax = patches[0].data(0, 0);
 
      for (const auto &patch : patches)
        {
          const unsigned int n_subdivisions = patch.n_subdivisions;
 
          Assert((patch.data.n_rows() == n_data_sets &&
                  !patch.points_are_available) ||
                   (patch.data.n_rows() == n_data_sets + spacedim &&
                    patch.points_are_available),
                 ExcDimensionMismatch(patch.points_are_available ?
                                        (n_data_sets + spacedim) :
                                        n_data_sets,
                                      patch.data.n_rows()));
          Assert(patch.data.n_cols() ==
                   Utilities::fixed_power<dim>(n_subdivisions + 1),
                 ExcInvalidDatasetSize(patch.data.n_cols(),
                                       n_subdivisions + 1));
 
          for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
            for (unsigned int j = 0; j < n_subdivisions + 1; ++j)
              {
                const int dl = i * (n_subdivisions + 1) + j;
                if (patch.data(0, dl) < hmin)
                  hmin = patch.data(0, dl);
                if (patch.data(0, dl) > hmax)
                  hmax = patch.data(0, dl);
              }
        }
 
      out << "#declare HMIN=" << hmin << ";" << '\n'
          << "#declare HMAX=" << hmax << ";" << '\n'
          << '\n';
 
      if (!flags.external_data)
        {
          // texture with scaled niveau lines 10 lines in the surface
          out << "#declare Tex=texture{" << '\n'
              << "  pigment {" << '\n'
              << "    gradient y" << '\n'
              << "    scale y*(HMAX-HMIN)*" << 0.1 << '\n'
              << "    color_map {" << '\n'
              << "      [0.00 color Light_Purple] " << '\n'
              << "      [0.95 color Light_Purple] " << '\n'
              << "      [1.00 color White]    " << '\n'
              << "} } }" << '\n'
              << '\n';
        }
 
      if (!flags.bicubic_patch)
        {
          // start of mesh header
          out << '\n' << "mesh {" << '\n';
        }
 
      // loop over all patches
      for (const auto &patch : patches)
        {
          const unsigned int n_subdivisions = patch.n_subdivisions;
          const unsigned int n              = n_subdivisions + 1;
          const unsigned int d1             = 1;
          const unsigned int d2             = n;
 
          Assert((patch.data.n_rows() == n_data_sets &&
                  !patch.points_are_available) ||
                   (patch.data.n_rows() == n_data_sets + spacedim &&
                    patch.points_are_available),
                 ExcDimensionMismatch(patch.points_are_available ?
                                        (n_data_sets + spacedim) :
                                        n_data_sets,
                                      patch.data.n_rows()));
          Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(n),
                 ExcInvalidDatasetSize(patch.data.n_cols(),
                                       n_subdivisions + 1));
 
 
          std::vector<Point<spacedim>> ver(n * n);
 
          for (unsigned int i2 = 0; i2 < n; ++i2)
            for (unsigned int i1 = 0; i1 < n; ++i1)
              {
                // compute coordinates for this patch point, storing in ver
                ver[i1 * d1 + i2 * d2] =
                  get_equispaced_location(patch, {i1, i2}, n_subdivisions);
              }
 
 
          if (!flags.bicubic_patch)
            {
              // approximate normal vectors in patch
              std::vector<Point<3>> nrml;
              // only if smooth triangles are used
              if (flags.smooth)
                {
                  nrml.resize(n * n);
                  // These are difference quotients of the surface
                  // mapping. We take them symmetric inside the
                  // patch and one-sided at the edges
                  Point<3> h1, h2;
                  // Now compute normals in every point
                  for (unsigned int i = 0; i < n; ++i)
                    for (unsigned int j = 0; j < n; ++j)
                      {
                        const unsigned int il = (i == 0) ? i : (i - 1);
                        const unsigned int ir =
                          (i == n_subdivisions) ? i : (i + 1);
                        const unsigned int jl = (j == 0) ? j : (j - 1);
                        const unsigned int jr =
                          (j == n_subdivisions) ? j : (j + 1);
 
                        h1(0) =
                          ver[ir * d1 + j * d2](0) - ver[il * d1 + j * d2](0);
                        h1(1) = patch.data(0, ir * d1 + j * d2) -
                                patch.data(0, il * d1 + j * d2);
                        h1(2) =
                          ver[ir * d1 + j * d2](1) - ver[il * d1 + j * d2](1);
 
                        h2(0) =
                          ver[i * d1 + jr * d2](0) - ver[i * d1 + jl * d2](0);
                        h2(1) = patch.data(0, i * d1 + jr * d2) -
                                patch.data(0, i * d1 + jl * d2);
                        h2(2) =
                          ver[i * d1 + jr * d2](1) - ver[i * d1 + jl * d2](1);
 
                        nrml[i * d1 + j * d2](0) =
                          h1(1) * h2(2) - h1(2) * h2(1);
                        nrml[i * d1 + j * d2](1) =
                          h1(2) * h2(0) - h1(0) * h2(2);
                        nrml[i * d1 + j * d2](2) =
                          h1(0) * h2(1) - h1(1) * h2(0);
 
                        // normalize Vector
                        double norm =
                          std::sqrt(std::pow(nrml[i * d1 + j * d2](0), 2.) +
                                    std::pow(nrml[i * d1 + j * d2](1), 2.) +
                                    std::pow(nrml[i * d1 + j * d2](2), 2.));
 
                        if (nrml[i * d1 + j * d2](1) < 0)
                          norm *= -1.;
 
                        for (unsigned int k = 0; k < 3; ++k)
                          nrml[i * d1 + j * d2](k) /= norm;
                      }
                }
 
              // setting up triangles
              for (unsigned int i = 0; i < n_subdivisions; ++i)
                for (unsigned int j = 0; j < n_subdivisions; ++j)
                  {
                    // down/left vertex of triangle
                    const int dl = i * d1 + j * d2;
                    if (flags.smooth)
                      {
                        // writing smooth_triangles
 
                        // down/right triangle
                        out << "smooth_triangle {" << '\n'
                            << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
                            << "," << ver[dl](1) << ">, <" << nrml[dl](0)
                            << ", " << nrml[dl](1) << ", " << nrml[dl](2)
                            << ">," << '\n';
                        out << " \t<" << ver[dl + d1](0) << ","
                            << patch.data(0, dl + d1) << "," << ver[dl + d1](1)
                            << ">, <" << nrml[dl + d1](0) << ", "
                            << nrml[dl + d1](1) << ", " << nrml[dl + d1](2)
                            << ">," << '\n';
                        out << "\t<" << ver[dl + d1 + d2](0) << ","
                            << patch.data(0, dl + d1 + d2) << ","
                            << ver[dl + d1 + d2](1) << ">, <"
                            << nrml[dl + d1 + d2](0) << ", "
                            << nrml[dl + d1 + d2](1) << ", "
                            << nrml[dl + d1 + d2](2) << ">}" << '\n';
 
                        // upper/left triangle
                        out << "smooth_triangle {" << '\n'
                            << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
                            << "," << ver[dl](1) << ">, <" << nrml[dl](0)
                            << ", " << nrml[dl](1) << ", " << nrml[dl](2)
                            << ">," << '\n';
                        out << "\t<" << ver[dl + d1 + d2](0) << ","
                            << patch.data(0, dl + d1 + d2) << ","
                            << ver[dl + d1 + d2](1) << ">, <"
                            << nrml[dl + d1 + d2](0) << ", "
                            << nrml[dl + d1 + d2](1) << ", "
                            << nrml[dl + d1 + d2](2) << ">," << '\n';
                        out << "\t<" << ver[dl + d2](0) << ","
                            << patch.data(0, dl + d2) << "," << ver[dl + d2](1)
                            << ">, <" << nrml[dl + d2](0) << ", "
                            << nrml[dl + d2](1) << ", " << nrml[dl + d2](2)
                            << ">}" << '\n';
                      }
                    else
                      {
                        // writing standard triangles down/right triangle
                        out << "triangle {" << '\n'
                            << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
                            << "," << ver[dl](1) << ">," << '\n';
                        out << "\t<" << ver[dl + d1](0) << ","
                            << patch.data(0, dl + d1) << "," << ver[dl + d1](1)
                            << ">," << '\n';
                        out << "\t<" << ver[dl + d1 + d2](0) << ","
                            << patch.data(0, dl + d1 + d2) << ","
                            << ver[dl + d1 + d2](1) << ">}" << '\n';
 
                        // upper/left triangle
                        out << "triangle {" << '\n'
                            << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
                            << "," << ver[dl](1) << ">," << '\n';
                        out << "\t<" << ver[dl + d1 + d2](0) << ","
                            << patch.data(0, dl + d1 + d2) << ","
                            << ver[dl + d1 + d2](1) << ">," << '\n';
                        out << "\t<" << ver[dl + d2](0) << ","
                            << patch.data(0, dl + d2) << "," << ver[dl + d2](1)
                            << ">}" << '\n';
                      }
                  }
            }
          else
            {
              // writing bicubic_patch
              Assert(n_subdivisions == 3,
                     ExcDimensionMismatch(n_subdivisions, 3));
              out << '\n'
                  << "bicubic_patch {" << '\n'
                  << "  type 0" << '\n'
                  << "  flatness 0" << '\n'
                  << "  u_steps 0" << '\n'
                  << "  v_steps 0" << '\n';
              for (int i = 0; i < 16; ++i)
                {
                  out << "\t<" << ver[i](0) << "," << patch.data(0, i) << ","
                      << ver[i](1) << ">";
                  if (i != 15)
                    out << ",";
                  out << '\n';
                }
              out << "  texture {Tex}" << '\n' << "}" << '\n';
            }
        }
 
      if (!flags.bicubic_patch)
        {
          // the end of the mesh
          out << "  texture {Tex}" << '\n' << "}" << '\n' << '\n';
        }
 
      // make sure everything now gets to disk
      out.flush();
 
      AssertThrow(out.fail() == false, ExcIO());
    }
  } // namespace
 
 
 
  template <int dim, int spacedim>
  void
  write_povray(
    const std::vector<Patch<dim, spacedim>> &patches,
    const std::vector<std::string> &         data_names,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>> &,
    const PovrayFlags &flags,
    std::ostream &     out)
  {
    do_write_povray(patches, data_names, flags, out);
  }
 
 
 
  template <int dim, int spacedim>
  void
  write_eps(
    const std::vector<Patch<dim, spacedim>> & /*patches*/,
    const std::vector<std::string> & /*data_names*/,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>> &,
    const EpsFlags & /*flags*/,
    std::ostream & /*out*/)
  {
    // not implemented, see the documentation of the function
    AssertThrow(dim == 2, ExcNotImplemented());
  }
 
 
  template <int spacedim>
  void
  write_eps(
    const std::vector<Patch<2, spacedim>> &patches,
    const std::vector<std::string> & /*data_names*/,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>> &,
    const EpsFlags &flags,
    std::ostream &  out)
  {
    AssertThrow(out.fail() == false, ExcIO());
 
#ifndef DEAL_II_WITH_MPI
    // verify that there are indeed patches to be written out. most of the
    // times, people just forget to call build_patches when there are no
    // patches, so a warning is in order. that said, the assertion is disabled
    // if we support MPI since then it can happen that on the coarsest mesh, a
    // processor simply has no cells it actually owns, and in that case it is
    // legit if there are no patches
    Assert(patches.size() > 0, ExcNoPatches());
#else
    if (patches.size() == 0)
      return;
#endif
 
    // set up an array of cells to be written later. this array holds the cells
    // of all the patches as projected to the plane perpendicular to the line of
    // sight.
    //
    // note that they are kept sorted by the set, where we chose the value of
    // the center point of the cell along the line of sight as value for sorting
    std::multiset<EpsCell2d> cells;
 
    // two variables in which we will store the minimum and maximum values of
    // the field to be used for colorization
    float min_color_value = std::numeric_limits<float>::max();
    float max_color_value = std::numeric_limits<float>::min();
 
    // Array for z-coordinates of points. The elevation determined by a function
    // if spacedim=2 or the z-cooridate of the grid point if spacedim=3
    double heights[4] = {0, 0, 0, 0};
 
    // compute the cells for output and enter them into the set above note that
    // since dim==2, we have exactly four vertices per patch and per cell
    for (const auto &patch : patches)
      {
        const unsigned int n_subdivisions = patch.n_subdivisions;
        const unsigned int n              = n_subdivisions + 1;
        const unsigned int d1             = 1;
        const unsigned int d2             = n;
 
        for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
          for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
            {
              Point<spacedim> points[4];
              points[0] =
                get_equispaced_location(patch, {i1, i2}, n_subdivisions);
              points[1] =
                get_equispaced_location(patch, {i1 + 1, i2}, n_subdivisions);
              points[2] =
                get_equispaced_location(patch, {i1, i2 + 1}, n_subdivisions);
              points[3] = get_equispaced_location(patch,
                                                  {i1 + 1, i2 + 1},
                                                  n_subdivisions);
 
              switch (spacedim)
                {
                  case 2:
                    Assert((flags.height_vector < patch.data.n_rows()) ||
                             patch.data.n_rows() == 0,
                           ExcIndexRange(flags.height_vector,
                                         0,
                                         patch.data.n_rows()));
                    heights[0] =
                      patch.data.n_rows() != 0 ?
                        patch.data(flags.height_vector, i1 * d1 + i2 * d2) *
                          flags.z_scaling :
                        0;
                    heights[1] = patch.data.n_rows() != 0 ?
                                   patch.data(flags.height_vector,
                                              (i1 + 1) * d1 + i2 * d2) *
                                     flags.z_scaling :
                                   0;
                    heights[2] = patch.data.n_rows() != 0 ?
                                   patch.data(flags.height_vector,
                                              i1 * d1 + (i2 + 1) * d2) *
                                     flags.z_scaling :
                                   0;
                    heights[3] = patch.data.n_rows() != 0 ?
                                   patch.data(flags.height_vector,
                                              (i1 + 1) * d1 + (i2 + 1) * d2) *
                                     flags.z_scaling :
                                   0;
 
                    break;
                  case 3:
                    // Copy z-coordinates into the height vector
                    for (unsigned int i = 0; i < 4; ++i)
                      heights[i] = points[i](2);
                    break;
                  default:
                    Assert(false, ExcNotImplemented());
                }
 
 
              // now compute the projection of the bilinear cell given by the
              // four vertices and their heights and write them to a proper cell
              // object. note that we only need the first two components of the
              // projected position for output, but we need the value along the
              // line of sight for sorting the cells for back-to- front-output
              //
              // this computation was first written by Stefan Nauber. please
              // no-one ask me why it works that way (or may be not), especially
              // not about the angles and the sign of the height field, I don't
              // know it.
              EpsCell2d    eps_cell;
              const double pi = numbers::PI;
              const double cx =
                             -std::cos(pi - flags.azimut_angle * 2 * pi / 360.),
                           cz = -std::cos(flags.turn_angle * 2 * pi / 360.),
                           sx =
                             std::sin(pi - flags.azimut_angle * 2 * pi / 360.),
                           sz = std::sin(flags.turn_angle * 2 * pi / 360.);
              for (unsigned int vertex = 0; vertex < 4; ++vertex)
                {
                  const double x = points[vertex](0), y = points[vertex](1),
                               z = -heights[vertex];
 
                  eps_cell.vertices[vertex](0) = -cz * x + sz * y;
                  eps_cell.vertices[vertex](1) =
                    -cx * sz * x - cx * cz * y - sx * z;
 
                  //      ( 1 0    0 )
                  // D1 = ( 0 cx -sx )
                  //      ( 0 sx  cx )
 
                  //      ( cy 0 sy )
                  // Dy = (  0 1  0 )
                  //      (-sy 0 cy )
 
                  //      ( cz -sz 0 )
                  // Dz = ( sz  cz 0 )
                  //      (  0   0 1 )
 
                  //       ( cz -sz 0 )( 1 0    0 )(x)   (
                  //       cz*x-sz*(cx*y-sx*z)+0*(sx*y+cx*z) )
                  // Dxz = ( sz  cz 0 )( 0 cx -sx )(y) = (
                  // sz*x+cz*(cx*y-sx*z)+0*(sx*y+cx*z) )
                  //       (  0   0 1 )( 0 sx  cx )(z)   (  0*x+
                  //       *(cx*y-sx*z)+1*(sx*y+cx*z) )
                }
 
              // compute coordinates of center of cell
              const Point<spacedim> center_point =
                (points[0] + points[1] + points[2] + points[3]) / 4;
              const double center_height =
                -(heights[0] + heights[1] + heights[2] + heights[3]) / 4;
 
              // compute the depth into the picture
              eps_cell.depth = -sx * sz * center_point(0) -
                               sx * cz * center_point(1) + cx * center_height;
 
              if (flags.draw_cells && flags.shade_cells)
                {
                  Assert((flags.color_vector < patch.data.n_rows()) ||
                           patch.data.n_rows() == 0,
                         ExcIndexRange(flags.color_vector,
                                       0,
                                       patch.data.n_rows()));
                  const double color_values[4] = {
                    patch.data.n_rows() != 0 ?
                      patch.data(flags.color_vector, i1 * d1 + i2 * d2) :
                      1,
 
                    patch.data.n_rows() != 0 ?
                      patch.data(flags.color_vector, (i1 + 1) * d1 + i2 * d2) :
                      1,
 
                    patch.data.n_rows() != 0 ?
                      patch.data(flags.color_vector, i1 * d1 + (i2 + 1) * d2) :
                      1,
 
                    patch.data.n_rows() != 0 ?
                      patch.data(flags.color_vector,
                                 (i1 + 1) * d1 + (i2 + 1) * d2) :
                      1};
 
                  // set color value to average of the value at the vertices
                  eps_cell.color_value = (color_values[0] + color_values[1] +
                                          color_values[3] + color_values[2]) /
                                         4;
 
                  // update bounds of color field
                  min_color_value =
                    std::min(min_color_value, eps_cell.color_value);
                  max_color_value =
                    std::max(max_color_value, eps_cell.color_value);
                }
 
              // finally add this cell
              cells.insert(eps_cell);
            }
      }
 
    // find out minimum and maximum x and y coordinates to compute offsets and
    // scaling factors
    double x_min = cells.begin()->vertices[0](0);
    double x_max = x_min;
    double y_min = cells.begin()->vertices[0](1);
    double y_max = y_min;
 
    for (const auto &cell : cells)
      for (const auto &vertex : cell.vertices)
        {
          x_min = std::min(x_min, vertex(0));
          x_max = std::max(x_max, vertex(0));
          y_min = std::min(y_min, vertex(1));
          y_max = std::max(y_max, vertex(1));
        }
 
    // scale in x-direction such that in the output 0 <= x <= 300. don't scale
    // in y-direction to preserve the shape of the triangulation
    const double scale =
      (flags.size /
       (flags.size_type == EpsFlags::width ? x_max - x_min : y_min - y_max));
 
    const Point<2> offset(x_min, y_min);
 
 
    // now write preamble
    {
      out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
          << "%%Title: deal.II Output" << '\n'
          << "%%Creator: the deal.II library" << '\n'
          << "%%Creation Date: " << Utilities::System::get_date() << " - "
          << Utilities::System::get_time() << '\n'
          << "%%BoundingBox: "
          // lower left corner
          << "0 0 "
          // upper right corner
          << static_cast<unsigned int>((x_max - x_min) * scale + 0.5) << ' '
          << static_cast<unsigned int>((y_max - y_min) * scale + 0.5) << '\n';
 
      // define some abbreviations to keep the output small:
      // m=move turtle to
      // l=define a line
      // s=set rgb color
      // sg=set gray value
      // lx=close the line and plot the line
      // lf=close the line and fill the interior
      out << "/m {moveto} bind def" << '\n'
          << "/l {lineto} bind def" << '\n'
          << "/s {setrgbcolor} bind def" << '\n'
          << "/sg {setgray} bind def" << '\n'
          << "/lx {lineto closepath stroke} bind def" << '\n'
          << "/lf {lineto closepath fill} bind def" << '\n';
 
      out << "%%EndProlog" << '\n' << '\n';
      // set fine lines
      out << flags.line_width << " setlinewidth" << '\n';
    }
 
    // check if min and max values for the color are actually different. If
    // that is not the case (such things happen, for example, in the very first
    // time step of a time dependent problem, if the initial values are zero),
    // all values are equal, and then we can draw everything in an arbitrary
    // color. Thus, change one of the two values arbitrarily
    if (max_color_value == min_color_value)
      max_color_value = min_color_value + 1;
 
    // now we've got all the information we need. write the cells. note: due to
    // the ordering, we traverse the list of cells back-to-front
    for (const auto &cell : cells)
      {
        if (flags.draw_cells)
          {
            if (flags.shade_cells)
              {
                const EpsFlags::RgbValues rgb_values =
                  (*flags.color_function)(cell.color_value,
                                          min_color_value,
                                          max_color_value);
 
                // write out color
                if (rgb_values.is_grey())
                  out << rgb_values.red << " sg ";
                else
                  out << rgb_values.red << ' ' << rgb_values.green << ' '
                      << rgb_values.blue << " s ";
              }
            else
              out << "1 sg ";
 
            out << (cell.vertices[0] - offset) * scale << " m "
                << (cell.vertices[1] - offset) * scale << " l "
                << (cell.vertices[3] - offset) * scale << " l "
                << (cell.vertices[2] - offset) * scale << " lf" << '\n';
          }
 
        if (flags.draw_mesh)
          out << "0 sg " // draw lines in black
              << (cell.vertices[0] - offset) * scale << " m "
              << (cell.vertices[1] - offset) * scale << " l "
              << (cell.vertices[3] - offset) * scale << " l "
              << (cell.vertices[2] - offset) * scale << " lx" << '\n';
      }
    out << "showpage" << '\n';
 
    out.flush();
 
    AssertThrow(out.fail() == false, ExcIO());
  }
 
 
 
  template <int dim, int spacedim>
  void
  write_gmv(
    const std::vector<Patch<dim, spacedim>> &patches,
    const std::vector<std::string> &         data_names,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>> &,
    const GmvFlags &flags,
    std::ostream &  out)
  {
    // The gmv format does not support cells that only consist of a single
    // point. It does support the output of point data using the keyword
    // 'tracers' instead of 'nodes' and 'cells', but this output format is
    // currently not implemented.
    AssertThrow(dim > 0, ExcNotImplemented());
 
    Assert(dim <= 3, ExcNotImplemented());
    AssertThrow(out.fail() == false, ExcIO());
 
#ifndef DEAL_II_WITH_MPI
    // verify that there are indeed patches to be written out. most of the
    // times, people just forget to call build_patches when there are no
    // patches, so a warning is in order. that said, the assertion is disabled
    // if we support MPI since then it can happen that on the coarsest mesh, a
    // processor simply has no cells it actually owns, and in that case it is
    // legit if there are no patches
    Assert(patches.size() > 0, ExcNoPatches());
#else
    if (patches.size() == 0)
      return;
#endif
 
    GmvStream          gmv_out(out, flags);
    const unsigned int n_data_sets = data_names.size();
    // check against # of data sets in first patch. checks against all other
    // patches are made in write_gmv_reorder_data_vectors
    Assert((patches[0].data.n_rows() == n_data_sets &&
            !patches[0].points_are_available) ||
             (patches[0].data.n_rows() == n_data_sets + spacedim &&
              patches[0].points_are_available),
           ExcDimensionMismatch(patches[0].points_are_available ?
                                  (n_data_sets + spacedim) :
                                  n_data_sets,
                                patches[0].data.n_rows()));
 
    //---------------------
    // preamble
    out << "gmvinput ascii" << '\n' << '\n';
 
    // first count the number of cells and cells for later use
    unsigned int n_nodes;
    unsigned int n_cells;
    compute_sizes<dim, spacedim>(patches, n_nodes, n_cells);
 
    // in gmv format the vertex coordinates and the data have an order that is a
    // bit unpleasant (first all x coordinates, then all y coordinate, ...;
    // first all data of variable 1, then variable 2, etc), so we have to copy
    // the data vectors a bit around
    //
    // note that we copy vectors when looping over the patches since we have to
    // write them one variable at a time and don't want to use more than one
    // loop
    //
    // this copying of data vectors can be done while we already output the
    // vertices, so do this on a separate task and when wanting to write out the
    // data, we wait for that task to finish
    Table<2, double> data_vectors(n_data_sets, n_nodes);
    void (*fun_ptr)(const std::vector<Patch<dim, spacedim>> &,
                    Table<2, double> &) =
      &write_gmv_reorder_data_vectors<dim, spacedim>;
    Threads::Task<> reorder_task =
      Threads::new_task(fun_ptr, patches, data_vectors);
 
    //-----------------------------
    // first make up a list of used vertices along with their coordinates
    //
    // note that we have to print 3 dimensions
    out << "nodes " << n_nodes << '\n';
    for (unsigned int d = 0; d < spacedim; ++d)
      {
        gmv_out.selected_component = d;
        write_nodes(patches, gmv_out);
        out << '\n';
      }
    gmv_out.selected_component = numbers::invalid_unsigned_int;
 
    for (unsigned int d = spacedim; d < 3; ++d)
      {
        for (unsigned int i = 0; i < n_nodes; ++i)
          out << "0 ";
        out << '\n';
      }
 
    //-------------------------------
    // now for the cells. note that vertices are counted from 1 onwards
    out << "cells " << n_cells << '\n';
    write_cells(patches, gmv_out);
 
    //-------------------------------------
    // data output.
    out << "variable" << '\n';
 
    // now write the data vectors to @p{out} first make sure that all data is in
    // place
    reorder_task.join();
 
    // then write data. the '1' means: node data (as opposed to cell data, which
    // we do not support explicitly here)
    for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
      {
        out << data_names[data_set] << " 1" << '\n';
        std::copy(data_vectors[data_set].begin(),
                  data_vectors[data_set].end(),
                  std::ostream_iterator<double>(out, " "));
        out << '\n' << '\n';
      }
 
 
 
    // end of variable section
    out << "endvars" << '\n';
 
    // end of output
    out << "endgmv" << '\n';
 
    // make sure everything now gets to disk
    out.flush();
 
    // assert the stream is still ok
    AssertThrow(out.fail() == false, ExcIO());
  }
 
 
 
  template <int dim, int spacedim>
  void
  write_tecplot(
    const std::vector<Patch<dim, spacedim>> &patches,
    const std::vector<std::string> &         data_names,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>> &,
    const TecplotFlags &flags,
    std::ostream &      out)
  {
    AssertThrow(out.fail() == false, ExcIO());
 
    // The FEBLOCK or FEPOINT formats of tecplot only allows full elements (e.g.
    // triangles), not single points. Other tecplot format allow point output,
    // but they are currently not implemented.
    AssertThrow(dim > 0, ExcNotImplemented());
 
#ifndef DEAL_II_WITH_MPI
    // verify that there are indeed patches to be written out. most of the
    // times, people just forget to call build_patches when there are no
    // patches, so a warning is in order. that said, the assertion is disabled
    // if we support MPI since then it can happen that on the coarsest mesh, a
    // processor simply has no cells it actually owns, and in that case it is
    // legit if there are no patches
    Assert(patches.size() > 0, ExcNoPatches());
#else
    if (patches.size() == 0)
      return;
#endif
 
    TecplotStream tecplot_out(out, flags);
 
    const unsigned int n_data_sets = data_names.size();
    // check against # of data sets in first patch. checks against all other
    // patches are made in write_gmv_reorder_data_vectors
    Assert((patches[0].data.n_rows() == n_data_sets &&
            !patches[0].points_are_available) ||
             (patches[0].data.n_rows() == n_data_sets + spacedim &&
              patches[0].points_are_available),
           ExcDimensionMismatch(patches[0].points_are_available ?
                                  (n_data_sets + spacedim) :
                                  n_data_sets,
                                patches[0].data.n_rows()));
 
    // first count the number of cells and cells for later use
    unsigned int n_nodes;
    unsigned int n_cells;
    compute_sizes<dim, spacedim>(patches, n_nodes, n_cells);
 
    //---------
    // preamble
    {
      out
        << "# This file was generated by the deal.II library." << '\n'
        << "# Date =  " << Utilities::System::get_date() << '\n'
        << "# Time =  " << Utilities::System::get_time() << '\n'
        << "#" << '\n'
        << "# For a description of the Tecplot format see the Tecplot documentation."
        << '\n'
        << "#" << '\n';
 
 
      out << "Variables=";
 
      switch (spacedim)
        {
          case 1:
            out << "\"x\"";
            break;
          case 2:
            out << "\"x\", \"y\"";
            break;
          case 3:
            out << "\"x\", \"y\", \"z\"";
            break;
          default:
            Assert(false, ExcNotImplemented());
        }
 
      for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
        out << ", \"" << data_names[data_set] << "\"";
 
      out << '\n';
 
      out << "zone ";
      if (flags.zone_name)
        out << "t=\"" << flags.zone_name << "\" ";
 
      if (flags.solution_time >= 0.0)
        out << "strandid=1, solutiontime=" << flags.solution_time << ", ";
 
      out << "f=feblock, n=" << n_nodes << ", e=" << n_cells
          << ", et=" << tecplot_cell_type[dim] << '\n';
    }
 
 
    // in Tecplot FEBLOCK format the vertex coordinates and the data have an
    // order that is a bit unpleasant (first all x coordinates, then all y
    // coordinate, ...; first all data of variable 1, then variable 2, etc), so
    // we have to copy the data vectors a bit around
    //
    // note that we copy vectors when looping over the patches since we have to
    // write them one variable at a time and don't want to use more than one
    // loop
    //
    // this copying of data vectors can be done while we already output the
    // vertices, so do this on a separate task and when wanting to write out the
    // data, we wait for that task to finish
 
    Table<2, double> data_vectors(n_data_sets, n_nodes);
 
    void (*fun_ptr)(const std::vector<Patch<dim, spacedim>> &,
                    Table<2, double> &) =
      &write_gmv_reorder_data_vectors<dim, spacedim>;
    Threads::Task<> reorder_task =
      Threads::new_task(fun_ptr, patches, data_vectors);
 
    //-----------------------------
    // first make up a list of used vertices along with their coordinates
 
 
    for (unsigned int d = 0; d < spacedim; ++d)
      {
        tecplot_out.selected_component = d;
        write_nodes(patches, tecplot_out);
        out << '\n';
      }
 
 
    //-------------------------------------
    // data output.
    //
    // now write the data vectors to @p{out} first make sure that all data is in
    // place
    reorder_task.join();
 
    // then write data.
    for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
      {
        std::copy(data_vectors[data_set].begin(),
                  data_vectors[data_set].end(),
                  std::ostream_iterator<double>(out, "\n"));
        out << '\n';
      }
 
    write_cells(patches, tecplot_out);
 
    // make sure everything now gets to disk
    out.flush();
 
    // assert the stream is still ok
    AssertThrow(out.fail() == false, ExcIO());
  }
 
 
 
  //---------------------------------------------------------------------------
  // Macros for handling Tecplot API data
 
#ifdef DEAL_II_HAVE_TECPLOT
 
  namespace
  {
    class TecplotMacros
    {
    public:
      TecplotMacros(const unsigned int n_nodes = 0,
                    const unsigned int n_vars  = 0,
                    const unsigned int n_cells = 0,
                    const unsigned int n_vert  = 0);
      ~TecplotMacros();
      float &
      nd(const unsigned int i, const unsigned int j);
      int &
                         cd(const unsigned int i, const unsigned int j);
      std::vector<float> nodalData;
      std::vector<int>   connData;
 
    private:
      unsigned int n_nodes;
      unsigned int n_vars;
      unsigned int n_cells;
      unsigned int n_vert;
    };
 
 
    inline TecplotMacros::TecplotMacros(const unsigned int n_nodes,
                                        const unsigned int n_vars,
                                        const unsigned int n_cells,
                                        const unsigned int n_vert)
      : n_nodes(n_nodes)
      , n_vars(n_vars)
      , n_cells(n_cells)
      , n_vert(n_vert)
    {
      nodalData.resize(n_nodes * n_vars);
      connData.resize(n_cells * n_vert);
    }
 
 
 
    inline TecplotMacros::~TecplotMacros()
    {}
 
 
 
    inline float &
    TecplotMacros::nd(const unsigned int i, const unsigned int j)
    {
      return nodalData[i * n_nodes + j];
    }
 
 
 
    inline int &
    TecplotMacros::cd(const unsigned int i, const unsigned int j)
    {
      return connData[i + j * n_vert];
    }
 
  } // namespace
 
 
#endif
  //---------------------------------------------------------------------------
 
 
 
  template <int dim, int spacedim>
  void
  write_tecplot_binary(
    const std::vector<Patch<dim, spacedim>> &patches,
    const std::vector<std::string> &         data_names,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>>
      &                 nonscalar_data_ranges,
    const TecplotFlags &flags,
    std::ostream &      out)
  {
    // The FEBLOCK or FEPOINT formats of tecplot only allows full elements (e.g.
    // triangles), not single points. Other tecplot format allow point output,
    // but they are currently not implemented.
    AssertThrow(dim > 0, ExcNotImplemented());
 
#ifndef DEAL_II_HAVE_TECPLOT
 
    // simply call the ASCII output function if the Tecplot API isn't present
    write_tecplot(patches, data_names, nonscalar_data_ranges, flags, out);
    return;
 
#else
 
    // Tecplot binary output only good for 2D & 3D
    if (dim == 1)
      {
        write_tecplot(patches, data_names, nonscalar_data_ranges, flags, out);
        return;
      }
 
    // if the user hasn't specified a file name we should call the ASCII
    // function and use the ostream @p{out} instead of doing something silly
    // later
    char *file_name = (char *)flags.tecplot_binary_file_name;
 
    if (file_name == nullptr)
      {
        // At least in debug mode we should tell users why they don't get
        // tecplot binary output
        Assert(false,
               ExcMessage("Specify the name of the tecplot_binary"
                          " file through the TecplotFlags interface."));
        write_tecplot(patches, data_names, nonscalar_data_ranges, flags, out);
        return;
      }
 
 
    AssertThrow(out.fail() == false, ExcIO());
 
#  ifndef DEAL_II_WITH_MPI
    // verify that there are indeed patches to be written out. most of the
    // times, people just forget to call build_patches when there are no
    // patches, so a warning is in order. that said, the assertion is disabled
    // if we support MPI since then it can happen that on the coarsest mesh, a
    // processor simply has no cells it actually owns, and in that case it is
    // legit if there are no patches
    Assert(patches.size() > 0, ExcNoPatches());
#  else
    if (patches.size() == 0)
      return;
#  endif
 
    const unsigned int n_data_sets = data_names.size();
    // check against # of data sets in first patch. checks against all other
    // patches are made in write_gmv_reorder_data_vectors
    Assert((patches[0].data.n_rows() == n_data_sets &&
            !patches[0].points_are_available) ||
             (patches[0].data.n_rows() == n_data_sets + spacedim &&
              patches[0].points_are_available),
           ExcDimensionMismatch(patches[0].points_are_available ?
                                  (n_data_sets + spacedim) :
                                  n_data_sets,
                                patches[0].data.n_rows()));
 
    // first count the number of cells and cells for later use
    unsigned int n_nodes;
    unsigned int n_cells;
    compute_sizes<dim, spacedim>(patches, n_nodes, n_cells);
    // local variables only needed to write Tecplot binary output files
    const unsigned int vars_per_node  = (spacedim + n_data_sets),
                       nodes_per_cell = GeometryInfo<dim>::vertices_per_cell;
 
    TecplotMacros tm(n_nodes, vars_per_node, n_cells, nodes_per_cell);
 
    int is_double = 0, tec_debug = 0, cell_type = tecplot_binary_cell_type[dim];
 
    std::string tec_var_names;
    switch (spacedim)
      {
        case 2:
          tec_var_names = "x y";
          break;
        case 3:
          tec_var_names = "x y z";
          break;
        default:
          Assert(false, ExcNotImplemented());
      }
 
    for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
      {
        tec_var_names += " ";
        tec_var_names += data_names[data_set];
      }
    // in Tecplot FEBLOCK format the vertex coordinates and the data have an
    // order that is a bit unpleasant (first all x coordinates, then all y
    // coordinate, ...; first all data of variable 1, then variable 2, etc), so
    // we have to copy the data vectors a bit around
    //
    // note that we copy vectors when looping over the patches since we have to
    // write them one variable at a time and don't want to use more than one
    // loop
    //
    // this copying of data vectors can be done while we already output the
    // vertices, so do this on a separate task and when wanting to write out the
    // data, we wait for that task to finish
    Table<2, double> data_vectors(n_data_sets, n_nodes);
 
    void (*fun_ptr)(const std::vector<Patch<dim, spacedim>> &,
                    Table<2, double> &) =
      &write_gmv_reorder_data_vectors<dim, spacedim>;
    Threads::Task<> reorder_task =
      Threads::new_task(fun_ptr, patches, data_vectors);
 
    //-----------------------------
    // first make up a list of used vertices along with their coordinates
    for (unsigned int d = 1; d <= spacedim; ++d)
      {
        unsigned int entry = 0;
 
        for (const auto &patch : patches)
          {
            const unsigned int n_subdivisions = patch.n_subdivisions;
 
            switch (dim)
              {
                case 2:
                  {
                    for (unsigned int j = 0; j < n_subdivisions + 1; ++j)
                      for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
                        {
                          const double x_frac = i * 1. / n_subdivisions,
                                       y_frac = j * 1. / n_subdivisions;
 
                          tm.nd((d - 1), entry) = static_cast<float>(
                            (((patch.vertices[1](d - 1) * x_frac) +
                              (patch.vertices[0](d - 1) * (1 - x_frac))) *
                               (1 - y_frac) +
                             ((patch.vertices[3](d - 1) * x_frac) +
                              (patch.vertices[2](d - 1) * (1 - x_frac))) *
                               y_frac));
                          entry++;
                        }
                    break;
                  }
 
                case 3:
                  {
                    for (unsigned int j = 0; j < n_subdivisions + 1; ++j)
                      for (unsigned int k = 0; k < n_subdivisions + 1; ++k)
                        for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
                          {
                            const double x_frac = i * 1. / n_subdivisions,
                                         y_frac = k * 1. / n_subdivisions,
                                         z_frac = j * 1. / n_subdivisions;
 
                            // compute coordinates for this patch point
                            tm.nd((d - 1), entry) = static_cast<float>(
                              ((((patch.vertices[1](d - 1) * x_frac) +
                                 (patch.vertices[0](d - 1) * (1 - x_frac))) *
                                  (1 - y_frac) +
                                ((patch.vertices[3](d - 1) * x_frac) +
                                 (patch.vertices[2](d - 1) * (1 - x_frac))) *
                                  y_frac) *
                                 (1 - z_frac) +
                               (((patch.vertices[5](d - 1) * x_frac) +
                                 (patch.vertices[4](d - 1) * (1 - x_frac))) *
                                  (1 - y_frac) +
                                ((patch.vertices[7](d - 1) * x_frac) +
                                 (patch.vertices[6](d - 1) * (1 - x_frac))) *
                                  y_frac) *
                                 z_frac));
                            entry++;
                          }
                    break;
                  }
 
                default:
                  Assert(false, ExcNotImplemented());
              }
          }
      }
 
 
    //-------------------------------------
    // data output.
    //
    reorder_task.join();
 
    // then write data.
    for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
      for (unsigned int entry = 0; entry < data_vectors[data_set].size();
           entry++)
        tm.nd((spacedim + data_set), entry) =
          static_cast<float>(data_vectors[data_set][entry]);
 
 
 
    //-------------------------------
    // now for the cells. note that vertices are counted from 1 onwards
    unsigned int first_vertex_of_patch = 0;
    unsigned int elem                  = 0;
 
    for (const auto &patch : patches)
      {
        const unsigned int n_subdivisions = patch.n_subdivisions;
        const unsigned int n              = n_subdivisions + 1;
        const unsigned int d1             = 1;
        const unsigned int d2             = n;
        const unsigned int d3             = n * n;
        // write out the cells making up this patch
        switch (dim)
          {
            case 2:
              {
                for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
                  for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
                    {
                      tm.cd(0, elem) =
                        first_vertex_of_patch + (i1)*d1 + (i2)*d2 + 1;
                      tm.cd(1, elem) =
                        first_vertex_of_patch + (i1 + 1) * d1 + (i2)*d2 + 1;
                      tm.cd(2, elem) = first_vertex_of_patch + (i1 + 1) * d1 +
                                       (i2 + 1) * d2 + 1;
                      tm.cd(3, elem) =
                        first_vertex_of_patch + (i1)*d1 + (i2 + 1) * d2 + 1;
 
                      elem++;
                    }
                break;
              }
 
            case 3:
              {
                for (unsigned int i3 = 0; i3 < n_subdivisions; ++i3)
                  for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
                    for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
                      {
                        // note: vertex indices start with 1!
 
 
                        tm.cd(0, elem) = first_vertex_of_patch + (i1)*d1 +
                                         (i2)*d2 + (i3)*d3 + 1;
                        tm.cd(1, elem) = first_vertex_of_patch + (i1 + 1) * d1 +
                                         (i2)*d2 + (i3)*d3 + 1;
                        tm.cd(2, elem) = first_vertex_of_patch + (i1 + 1) * d1 +
                                         (i2 + 1) * d2 + (i3)*d3 + 1;
                        tm.cd(3, elem) = first_vertex_of_patch + (i1)*d1 +
                                         (i2 + 1) * d2 + (i3)*d3 + 1;
                        tm.cd(4, elem) = first_vertex_of_patch + (i1)*d1 +
                                         (i2)*d2 + (i3 + 1) * d3 + 1;
                        tm.cd(5, elem) = first_vertex_of_patch + (i1 + 1) * d1 +
                                         (i2)*d2 + (i3 + 1) * d3 + 1;
                        tm.cd(6, elem) = first_vertex_of_patch + (i1 + 1) * d1 +
                                         (i2 + 1) * d2 + (i3 + 1) * d3 + 1;
                        tm.cd(7, elem) = first_vertex_of_patch + (i1)*d1 +
                                         (i2 + 1) * d2 + (i3 + 1) * d3 + 1;
 
                        elem++;
                      }
                break;
              }
 
            default:
              Assert(false, ExcNotImplemented());
          }
 
 
        // finally update the number of the first vertex of this patch
        first_vertex_of_patch += Utilities::fixed_power<dim>(n);
      }
 
 
    {
      int ierr = 0, num_nodes = static_cast<int>(n_nodes),
          num_cells = static_cast<int>(n_cells);
 
      char dot[2] = {'.', 0};
      // Unfortunately, TECINI takes a char *, but c_str() gives a const char *.
      // As we don't do anything else with tec_var_names following const_cast is
      // ok
      char *var_names = const_cast<char *>(tec_var_names.c_str());
      ierr = TECINI(nullptr, var_names, file_name, dot, &tec_debug, &is_double);
 
      Assert(ierr == 0, ExcErrorOpeningTecplotFile(file_name));
 
      char FEBLOCK[] = {'F', 'E', 'B', 'L', 'O', 'C', 'K', 0};
      ierr =
        TECZNE(nullptr, &num_nodes, &num_cells, &cell_type, FEBLOCK, nullptr);
 
      Assert(ierr == 0, ExcTecplotAPIError());
 
      int total = (vars_per_node * num_nodes);
 
      ierr = TECDAT(&total, tm.nodalData.data(), &is_double);
 
      Assert(ierr == 0, ExcTecplotAPIError());
 
      ierr = TECNOD(tm.connData.data());
 
      Assert(ierr == 0, ExcTecplotAPIError());
 
      ierr = TECEND();
 
      Assert(ierr == 0, ExcTecplotAPIError());
    }
#endif
  }
 
 
 
  template <int dim, int spacedim>
  void
  write_vtk(
    const std::vector<Patch<dim, spacedim>> &patches,
    const std::vector<std::string> &         data_names,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>>
      &             nonscalar_data_ranges,
    const VtkFlags &flags,
    std::ostream &  out)
  {
    AssertThrow(out.fail() == false, ExcIO());
 
#ifndef DEAL_II_WITH_MPI
    // verify that there are indeed patches to be written out. most of the
    // times, people just forget to call build_patches when there are no
    // patches, so a warning is in order. that said, the assertion is disabled
    // if we support MPI since then it can happen that on the coarsest mesh, a
    // processor simply has no cells it actually owns, and in that case it is
    // legit if there are no patches
    Assert(patches.size() > 0, ExcNoPatches());
#else
    if (patches.size() == 0)
      return;
#endif
 
    VtkStream vtk_out(out, flags);
 
    const unsigned int n_data_sets = data_names.size();
    // check against # of data sets in first patch.
    if (patches[0].points_are_available)
      {
        AssertDimension(n_data_sets + spacedim, patches[0].data.n_rows())
      }
    else
      {
        AssertDimension(n_data_sets, patches[0].data.n_rows())
      }
 
    //---------------------
    // preamble
    {
      out << "# vtk DataFile Version 3.0" << '\n'
          << "#This file was generated by the deal.II library";
      if (flags.print_date_and_time)
        {
          out << " on " << Utilities::System::get_date() << " at "
              << Utilities::System::get_time();
        }
      else
        out << '.';
      out << '\n' << "ASCII" << '\n';
      // now output the data header
      out << "DATASET UNSTRUCTURED_GRID\n" << '\n';
    }
 
    // if desired, output time and cycle of the simulation, following the
    // instructions at
    // http://www.visitusers.org/index.php?title=Time_and_Cycle_in_VTK_files
    {
      const unsigned int n_metadata =
        ((flags.cycle != std::numeric_limits<unsigned int>::min() ? 1 : 0) +
         (flags.time != std::numeric_limits<double>::min() ? 1 : 0));
      if (n_metadata > 0)
        {
          out << "FIELD FieldData " << n_metadata << '\n';
 
          if (flags.cycle != std::numeric_limits<unsigned int>::min())
            {
              out << "CYCLE 1 1 int\n" << flags.cycle << '\n';
            }
          if (flags.time != std::numeric_limits<double>::min())
            {
              out << "TIME 1 1 double\n" << flags.time << '\n';
            }
        }
    }
 
    // first count the number of cells and cells for later use
    unsigned int n_nodes, n_cells, n_points_and_n_cells;
    compute_sizes(patches,
                  flags.write_higher_order_cells,
                  n_nodes,
                  n_cells,
                  n_points_and_n_cells);
 
    // in gmv format the vertex coordinates and the data have an order that is a
    // bit unpleasant (first all x coordinates, then all y coordinate, ...;
    // first all data of variable 1, then variable 2, etc), so we have to copy
    // the data vectors a bit around
    //
    // note that we copy vectors when looping over the patches since we have to
    // write them one variable at a time and don't want to use more than one
    // loop
    //
    // this copying of data vectors can be done while we already output the
    // vertices, so do this on a separate task and when wanting to write out the
    // data, we wait for that task to finish
    Table<2, double> data_vectors(n_data_sets, n_nodes);
 
    void (*fun_ptr)(const std::vector<Patch<dim, spacedim>> &,
                    Table<2, double> &) =
      &write_gmv_reorder_data_vectors<dim, spacedim>;
    Threads::Task<> reorder_task =
      Threads::new_task(fun_ptr, patches, data_vectors);
 
    //-----------------------------
    // first make up a list of used vertices along with their coordinates
    //
    // note that we have to print d=1..3 dimensions
    out << "POINTS " << n_nodes << " double" << '\n';
    write_nodes(patches, vtk_out);
    out << '\n';
    //-------------------------------
    // now for the cells
    out << "CELLS " << n_cells << ' ' << n_points_and_n_cells << '\n';
    if (flags.write_higher_order_cells)
      write_high_order_cells(patches, vtk_out);
    else
      write_cells(patches, vtk_out);
    out << '\n';
    // next output the types of the cells. since all cells are the same, this is
    // simple
    out << "CELL_TYPES " << n_cells << '\n';
 
    // need to distinguish between linear cells, simplex cells (linear or
    // quadratic), and  high order cells
    for (const auto &patch : patches)
      {
        const auto vtk_cell_id =
          extract_vtk_patch_info(patch, flags.write_higher_order_cells);
 
        for (unsigned int i = 0; i < vtk_cell_id[1]; ++i)
          out << ' ' << vtk_cell_id[0];
      }
 
    out << '\n';
    //-------------------------------------
    // data output.
 
    // now write the data vectors to @p{out} first make sure that all data is in
    // place
    reorder_task.join();
 
    // then write data.  the 'POINT_DATA' means: node data (as opposed to cell
    // data, which we do not support explicitly here). all following data sets
    // are point data
    out << "POINT_DATA " << n_nodes << '\n';
 
    // when writing, first write out all vector data, then handle the scalar
    // data sets that have been left over
    std::vector<bool> data_set_written(n_data_sets, false);
    for (const auto &nonscalar_data_range : nonscalar_data_ranges)
      {
        AssertThrow(std::get<3>(nonscalar_data_range) !=
                      DataComponentInterpretation::component_is_part_of_tensor,
                    ExcNotImplemented());
 
        AssertThrow(std::get<1>(nonscalar_data_range) >=
                      std::get<0>(nonscalar_data_range),
                    ExcLowerRange(std::get<1>(nonscalar_data_range),
                                  std::get<0>(nonscalar_data_range)));
        AssertThrow(std::get<1>(nonscalar_data_range) < n_data_sets,
                    ExcIndexRange(std::get<1>(nonscalar_data_range),
                                  0,
                                  n_data_sets));
        AssertThrow(std::get<1>(nonscalar_data_range) + 1 -
                        std::get<0>(nonscalar_data_range) <=
                      3,
                    ExcMessage(
                      "Can't declare a vector with more than 3 components "
                      "in VTK"));
 
        // mark these components as already written:
        for (unsigned int i = std::get<0>(nonscalar_data_range);
             i <= std::get<1>(nonscalar_data_range);
             ++i)
          data_set_written[i] = true;
 
        // write the header. concatenate all the component names with double
        // underscores unless a vector name has been specified
        out << "VECTORS ";
 
        if (!std::get<2>(nonscalar_data_range).empty())
          out << std::get<2>(nonscalar_data_range);
        else
          {
            for (unsigned int i = std::get<0>(nonscalar_data_range);
                 i < std::get<1>(nonscalar_data_range);
                 ++i)
              out << data_names[i] << "__";
            out << data_names[std::get<1>(nonscalar_data_range)];
          }
 
        out << " double" << '\n';
 
        // now write data. pad all vectors to have three components
        for (unsigned int n = 0; n < n_nodes; ++n)
          {
            switch (std::get<1>(nonscalar_data_range) -
                    std::get<0>(nonscalar_data_range))
              {
                case 0:
                  out << data_vectors(std::get<0>(nonscalar_data_range), n)
                      << " 0 0" << '\n';
                  break;
 
                case 1:
                  out << data_vectors(std::get<0>(nonscalar_data_range), n)
                      << ' '
                      << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
                      << " 0" << '\n';
                  break;
                case 2:
                  out << data_vectors(std::get<0>(nonscalar_data_range), n)
                      << ' '
                      << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
                      << ' '
                      << data_vectors(std::get<0>(nonscalar_data_range) + 2, n)
                      << '\n';
                  break;
 
                default:
                  // VTK doesn't support anything else than vectors with 1, 2,
                  // or 3 components
                  Assert(false, ExcInternalError());
              }
          }
      }
 
    // now do the left over scalar data sets
    for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
      if (data_set_written[data_set] == false)
        {
          out << "SCALARS " << data_names[data_set] << " double 1" << '\n'
              << "LOOKUP_TABLE default" << '\n';
          std::copy(data_vectors[data_set].begin(),
                    data_vectors[data_set].end(),
                    std::ostream_iterator<double>(out, " "));
          out << '\n';
        }
 
    // make sure everything now gets to disk
    out.flush();
 
    // assert the stream is still ok
    AssertThrow(out.fail() == false, ExcIO());
  }
 
 
  void
  write_vtu_header(std::ostream &out, const VtkFlags &flags)
  {
    AssertThrow(out.fail() == false, ExcIO());
    out << "<?xml version=\"1.0\" ?> \n";
    out << "<!-- \n";
    out << "# vtk DataFile Version 3.0" << '\n'
        << "#This file was generated by the deal.II library";
    if (flags.print_date_and_time)
      {
        out << " on " << Utilities::System::get_time() << " at "
            << Utilities::System::get_date();
      }
    else
      out << '.';
    out << "\n-->\n";
    out << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
#ifdef DEAL_II_WITH_ZLIB
    out << " compressor=\"vtkZLibDataCompressor\"";
#endif
#ifdef DEAL_II_WORDS_BIGENDIAN
    out << " byte_order=\"BigEndian\"";
#else
    out << " byte_order=\"LittleEndian\"";
#endif
    out << ">";
    out << '\n';
    out << "<UnstructuredGrid>";
    out << '\n';
  }
 
 
 
  void
  write_vtu_footer(std::ostream &out)
  {
    AssertThrow(out.fail() == false, ExcIO());
    out << " </UnstructuredGrid>\n";
    out << "</VTKFile>\n";
  }
 
 
 
  template <int dim, int spacedim>
  void
  write_vtu(
    const std::vector<Patch<dim, spacedim>> &patches,
    const std::vector<std::string> &         data_names,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>>
      &             nonscalar_data_ranges,
    const VtkFlags &flags,
    std::ostream &  out)
  {
    write_vtu_header(out, flags);
    write_vtu_main(patches, data_names, nonscalar_data_ranges, flags, out);
    write_vtu_footer(out);
 
    out << std::flush;
  }
 
 
  template <int dim, int spacedim>
  void
  write_vtu_main(
    const std::vector<Patch<dim, spacedim>> &patches,
    const std::vector<std::string> &         data_names,
    const std::vector<
      std::tuple<unsigned int,
                 unsigned int,
                 std::string,
                 DataComponentInterpretation::DataComponentInterpretation>>
      &             nonscalar_data_ranges,
    const VtkFlags &flags,
    std::ostream &  out)
  {
    AssertThrow(out.fail() == false, ExcIO());
 
    // If the user provided physical units, make sure that they don't contain
    // quote characters as this would make the VTU file invalid XML and
    // probably lead to all sorts of difficult error messages. Other than that,
    // trust the user that whatever they provide makes sense somehow.
    for (const auto &unit : flags.physical_units)
      {
        (void)unit;
        Assert(
          unit.second.find('\"') == std::string::npos,
          ExcMessage(
            "A physical unit you provided, <" + unit.second +
            ">, contained a quotation mark character. This is not allowed."));
      }
 
#ifndef DEAL_II_WITH_MPI
    // verify that there are indeed patches to be written out. most of the
    // times, people just forget to call build_patches when there are no
    // patches, so a warning is in order. that said, the assertion is disabled
    // if we support MPI since then it can happen that on the coarsest mesh, a
    // processor simply has no cells it actually owns, and in that case it is
    // legit if there are no patches
    Assert(patches.size() > 0, ExcNoPatches());
#else
    if (patches.size() == 0)
      {
        // we still need to output a valid vtu file, because other CPUs might
        // output data. This is the minimal file that is accepted by paraview
        // and visit. if we remove the field definitions, visit is complaining.
        out << "<Piece NumberOfPoints=\"0\" NumberOfCells=\"0\" >\n"
            << "<Cells>\n"
            << "<DataArray type=\"UInt8\" Name=\"types\"></DataArray>\n"
            << "</Cells>\n"
            << "  <PointData Scalars=\"scalars\">\n";
        std::vector<bool> data_set_written(data_names.size(), false);
        for (const auto &nonscalar_data_range : nonscalar_data_ranges)
          {
            // mark these components as already written:
            for (unsigned int i = std::get<0>(nonscalar_data_range);
                 i <= std::get<1>(nonscalar_data_range);
                 ++i)
              data_set_written[i] = true;
 
            // write the header. concatenate all the component names with double
            // underscores unless a vector name has been specified
            out << "    <DataArray type=\"Float32\" Name=\"";
 
            if (!std::get<2>(nonscalar_data_range).empty())
              out << std::get<2>(nonscalar_data_range);
            else
              {
                for (unsigned int i = std::get<0>(nonscalar_data_range);
                     i < std::get<1>(nonscalar_data_range);
                     ++i)
                  out << data_names[i] << "__";
                out << data_names[std::get<1>(nonscalar_data_range)];
              }
 
            out << "\" NumberOfComponents=\"3\"></DataArray>\n";
          }
 
        for (unsigned int data_set = 0; data_set < data_names.size();
             ++data_set)
          if (data_set_written[data_set] == false)
            {
              out << "    <DataArray type=\"Float32\" Name=\""
                  << data_names[data_set] << "\"></DataArray>\n";
            }
 
        out << "  </PointData>\n";
        out << "</Piece>\n";
 
        out << std::flush;
 
        return;
      }
#endif
 
    // first up: metadata
    //
    // if desired, output time and cycle of the simulation, following the
    // instructions at
    // http://www.visitusers.org/index.php?title=Time_and_Cycle_in_VTK_files
    {
      const unsigned int n_metadata =
        ((flags.cycle != std::numeric_limits<unsigned int>::min() ? 1 : 0) +
         (flags.time != std::numeric_limits<double>::min() ? 1 : 0));
      if (n_metadata > 0)
        out << "<FieldData>\n";
 
      if (flags.cycle != std::numeric_limits<unsigned int>::min())
        {
          out
            << "<DataArray type=\"Float32\" Name=\"CYCLE\" NumberOfTuples=\"1\" format=\"ascii\">"
            << flags.cycle << "</DataArray>\n";
        }
      if (flags.time != std::numeric_limits<double>::min())
        {
          out
            << "<DataArray type=\"Float32\" Name=\"TIME\" NumberOfTuples=\"1\" format=\"ascii\">"
            << flags.time << "</DataArray>\n";
        }
 
      if (n_metadata > 0)
        out << "</FieldData>\n";
    }
 
 
    VtuStream vtu_out(out, flags);
 
    const unsigned int n_data_sets = data_names.size();
    // check against # of data sets in first patch. checks against all other
    // patches are made in write_gmv_reorder_data_vectors
    if (patches[0].points_are_available)
      {
        AssertDimension(n_data_sets + spacedim, patches[0].data.n_rows())
      }
    else
      {
        AssertDimension(n_data_sets, patches[0].data.n_rows())
      }
 
#ifdef DEAL_II_WITH_ZLIB
    const char *ascii_or_binary = "binary";
#else
    const char *              ascii_or_binary = "ascii";
#endif
 
 
    // first count the number of cells and cells for later use
    unsigned int n_nodes, n_cells, n_points_and_n_cells;
    compute_sizes(patches,
                  flags.write_higher_order_cells,
                  n_nodes,
                  n_cells,
                  n_points_and_n_cells);
 
    // in gmv format the vertex coordinates and the data have an order that is a
    // bit unpleasant (first all x coordinates, then all y coordinate, ...;
    // first all data of variable 1, then variable 2, etc), so we have to copy
    // the data vectors a bit around
    //
    // note that we copy vectors when looping over the patches since we have to
    // write them one variable at a time and don't want to use more than one
    // loop
    //
    // this copying of data vectors can be done while we already output the
    // vertices, so do this on a separate task and when wanting to write out the
    // data, we wait for that task to finish
    Table<2, float> data_vectors(n_data_sets, n_nodes);
 
    void (*fun_ptr)(const std::vector<Patch<dim, spacedim>> &,
                    Table<2, float> &) =
      &write_gmv_reorder_data_vectors<dim, spacedim, float>;
    Threads::Task<> reorder_task =
      Threads::new_task(fun_ptr, patches, data_vectors);
 
    //-----------------------------
    // first make up a list of used vertices along with their coordinates
    //
    // note that according to the standard, we have to print d=1..3 dimensions,
    // even if we are in reality in 2d, for example
    out << "<Piece NumberOfPoints=\"" << n_nodes << "\" NumberOfCells=\""
        << n_cells << "\" >\n";
    out << "  <Points>\n";
    out << "    <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\""
        << ascii_or_binary << "\">\n";
    write_nodes(patches, vtu_out);
    out << "    </DataArray>\n";
    out << "  </Points>\n\n";
    //-------------------------------
    // now for the cells
    out << "  <Cells>\n";
    out << "    <DataArray type=\"Int32\" Name=\"connectivity\" format=\""
        << ascii_or_binary << "\">\n";
    if (flags.write_higher_order_cells)
      write_high_order_cells(patches, vtu_out);
    else
      write_cells(patches, vtu_out);
    out << "    </DataArray>\n";
 
    // XML VTU format uses offsets; this is different than the VTK format, which
    // puts the number of nodes per cell in front of the connectivity list.
    out << "    <DataArray type=\"Int32\" Name=\"offsets\" format=\""
        << ascii_or_binary << "\">\n";
 
    std::vector<int32_t> offsets;
    offsets.reserve(n_cells);
 
    // std::uint8_t might be an alias to unsigned char which is then not printed
    // as ascii integers
#ifdef DEAL_II_WITH_ZLIB
    std::vector<std::uint8_t> cell_types;
#else
    std::vector<unsigned int> cell_types;
#endif
    cell_types.reserve(n_cells);
 
    unsigned int first_vertex_of_patch = 0;
 
    for (const auto &patch : patches)
      {
        const auto vtk_cell_id =
          extract_vtk_patch_info(patch, flags.write_higher_order_cells);
 
        for (unsigned int i = 0; i < vtk_cell_id[1]; ++i)
          {
            cell_types.push_back(vtk_cell_id[0]);
            first_vertex_of_patch += vtk_cell_id[2];
            offsets.push_back(first_vertex_of_patch);
          }
      }
 
    vtu_out << offsets;
    out << '\n';
    out << "    </DataArray>\n";
 
    // next output the types of the cells. since all cells are the same, this is
    // simple
    out << "    <DataArray type=\"UInt8\" Name=\"types\" format=\""
        << ascii_or_binary << "\">\n";
 
    // this should compress well :-)
    vtu_out << cell_types;
    out << '\n';
    out << "    </DataArray>\n";
    out << "  </Cells>\n";
 
 
    //-------------------------------------
    // data output.
 
    // now write the data vectors to @p{out} first make sure that all data is in
    // place
    reorder_task.join();
 
    // then write data.  the 'POINT_DATA' means: node data (as opposed to cell
    // data, which we do not support explicitly here). all following data sets
    // are point data
    out << "  <PointData Scalars=\"scalars\">\n";
 
    // when writing, first write out all vector data, then handle the scalar
    // data sets that have been left over
    std::vector<bool> data_set_written(n_data_sets, false);
    for (const auto &range : nonscalar_data_ranges)
      {
        const auto  first_component = std::get<0>(range);
        const auto  last_component  = std::get<1>(range);
        const auto &name            = std::get<2>(range);
        const bool  is_tensor =
          (std::get<3>(range) ==
           DataComponentInterpretation::component_is_part_of_tensor);
        const unsigned int n_components = (is_tensor ? 9 : 3);
        AssertThrow(last_component >= first_component,
                    ExcLowerRange(last_component, first_component));
        AssertThrow(last_component < n_data_sets,
                    ExcIndexRange(last_component, 0, n_data_sets));
        if (is_tensor)
          {
            AssertThrow((last_component + 1 - first_component <= 9),
                        ExcMessage(
                          "Can't declare a tensor with more than 9 components "
                          "in VTK/VTU format."));
          }
        else
          {
            AssertThrow((last_component + 1 - first_component <= 3),
                        ExcMessage(
                          "Can't declare a vector with more than 3 components "
                          "in VTK/VTU format."));
          }
 
        // mark these components as already written:
        for (unsigned int i = first_component; i <= last_component; ++i)
          data_set_written[i] = true;