vtkwriter.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2015-2018 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
 
#ifndef VTKWRITER_HH
#define VTKWRITER_HH
 
#include <algorithm>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <memory> 
#include <sstream>
#include <type_traits>
#include <vector>
 
#include <boost/fusion/include/zip.hpp>
#include <boost/range/iterator_range.hpp>
 
#include <dune/grid/common/grid.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/referenceelements.hh>
 
#include "fem/shapefunctioncache.hh"
#include "io/iobase.hh"
#include "utilities/detailed_exception.hh"
#include "utilities/timing.hh"
#undef min
 
 
 
 
 
namespace Kaskade
{
  namespace FunctionSpace_Detail {
    // forward declaration
    template <class> class ToDomainRepresentation;
  }
 
  namespace VTKWriterDetail
  {
    // returns the textual VTK notation of this system's endianness.
    std::string vtkEndianness();
    
    // Numberings of mesh entity types in VTK files.
    // See the VTKk source code for definitions:
    // https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html
    enum VTKGeometryType
    {
      vtkLine = 3,
      vtkTriangle = 5,
      vtkQuadrilateral = 9,
      vtkTetrahedron = 10,
      vtkHexahedron = 12,
      vtkPrism = 13,
      vtkPyramid = 14,
      vtkquadLine = 21,
      vtkquadTriangle = 22,
      vtkquadQuadrilateral = 23,
      vtkquadTetrahedron = 24,
      vtkquadHexahedron = 25,
      vtkLagrangeTriangle = 69,
      vtkLagrangeQuadrilateral = 70,
      vtkLagrangeTetrahedron = 71,
      vtkLagrangeHexahedron = 72
    };
    
    VTKGeometryType vtkType(Dune::GeometryType const& t, int order);
    
    std::pair<int,int> vtkPointToDune(Dune::GeometryType const& type, int vtkPoint);
    
 
    // forward declaration
    class Base64Writer;
    
    // T  is a scalar type (float,double,int,...)
    template <class T>
    class VTKDataArrayWriter
    {
    public:
      VTKDataArrayWriter(std::ostream &s, IoOptions::OutputType outputType,
                         std::string const& name, int ncomps,
                         size_t entries, int indentCount, int precision);
      
      template <int m>
      void write (Dune::FieldVector<T,m> const&);
      
      void write(T);
      
      ~VTKDataArrayWriter();
      
    private:
      std::ostream& s;
      IoOptions::OutputType outputType;
      int ncomps;
      int counter;
      int numPerLine;
      int indentCount;
      int precision;
      std::unique_ptr<Base64Writer> base64;
    };
 
    // ---------------------------------------------------------------------------------------------
 
    void indent(std::ostream& s, int indentCount);
 
    // ---------------------------------------------------------------------------------------------
 
    // A quadrature rule with equally weighted quadrature nodes at the VTK points for
    // first order (corners only) and second order (corners and edge midpoints).
    // This is not intended for actual quadrature, but to speed up evaluation of shape functions
    // in the evaluator used when evaluating point values.
    template <class ct, int dim>
    class CornerQuadratureRule: public Dune::QuadratureRule<ct,dim>
    {
    public:
      CornerQuadratureRule(Dune::GeometryType gt, int order)
      : Dune::QuadratureRule<ct,dim>(gt,order)
      {
        assert(1<=order && order<=2);
 
        auto refCell = Dune::referenceElement<ct,dim>(gt);
 
        // Number of points. All vertices and, if order is 2, also all edge midpoints.
        int const n = refCell.size(dim) +  (order==2 ? refCell.size(dim-1) : 0);
 
        // Trivial integration weights (provided even if not needed here).
        double w = 1.0/n;
 
        // Step through all affected VTK points on the cell and record their position.
        for (int i=0; i<n; ++i)
        {
          auto [k,codim] = vtkPointToDune(gt,i);
          this->push_back(Dune::QuadraturePoint<ct,dim>(refCell.position(k,codim),w));
        }
      }
    };
    
    // ---------------------------------------------------------------------------------------------
 
    template <class GridView> 
    struct VTKGridInfo
    {
      using Cell         = typename GridView::template Codim<0>::Entity;
      using IndexSet     = typename GridView::IndexSet;
      using CellIterator = typename GridView::template Codim<0>::Iterator;
      
      static int constexpr dim = GridView::dimension;
      
      
      VTKGridInfo(GridView const& gridView_, IoOptions const& options_)
      : gridView(gridView_),
        cells(gridView.template begin<0>(),gridView.template end<0>()),
        ncells(gridView.size(0)),
        options(options_)
      {        
        // For counting the corners (i.e. all corners of all cells) we step through the 
        // cells and count their respective corners (and edges for order two). This is 
        // general enough to treat mixed grids containing different cell types.
        ncorners = 0;
        for (auto const& cell: cells) 
        {
          ncorners += cell.subEntities(dim);
          if (options.order==2)
            ncorners += cell.subEntities(dim-1);
        }
        
        // In conforming mode, the number of points is just the number of vertices 
        // (plus edges for oder two). Otherwise, points and corners coincide.
        if (options.dataMode==IoOptions::conforming)
        {
          npoints = gridView.size(dim);
          if (options.order==2)
            npoints += gridView.size(dim-1);
        }
        else
          npoints = ncorners;
      }
      
 
      // returns entity center coordinates
      template <class Cell>
      auto center(Cell const& cell, int subindex, int codim) const
      {
        return Dune::ReferenceElements<typename GridView::ctype,dim>::general(cell.type()).position(subindex,codim);
      }
      
      // returns a globally unique index for vertices and edges in conforming mode. 
      // Vertices come first, then edges.
      template <class Cell>
      int conformingIndex(Cell const& cell, int subindex, int codim) const
      {
        // due to the sorting vertices first then edges there is no need to switch between order 1 and 2
        if (codim==dim)
          return gridView.indexSet().subIndex(cell,subindex,codim);
        else
          return gridView.size(dim) + gridView.indexSet().subIndex(cell,subindex,codim);
      }
      
      // Writes the point coordinates to the VTK XML stream
      template <class Scalar=typename GridView::ctype>
      void writeGridPoints (int indentCount, std::ostream& s) const
      {
        VTKDataArrayWriter<Scalar> p(s, options.outputType, "Coordinates", 3,
                                     npoints, indentCount, options.precision);
        using Vector = Dune::FieldVector<Scalar,GridView::dimensionworld>;
        
        if (options.dataMode==IoOptions::conforming)
        {
          // We need to iterate over the cells to visit each point, but visit points multiple times this way.
          // We just collect the coordinates in a buffer sorted according to the point's index as defined
          // by conformingIndex(), and simply overwrite the previous data -- in conforming grids, it is
          // anyways the same.
          std::vector<Vector> xs(npoints);
          for (auto const& cell: cells)
          {
            // Simplex vertices come first, independent of the order.
            for (int corner=0; corner<cell.subEntities(dim); ++corner)   // global coordinates of all corners
              xs[conformingIndex(cell,corner,dim)] = cell.geometry().corner(corner); 
            
            // For second order, the edge midpoints are appended.
            if (options.order == 2)
              for (int edge=0; edge<cell.subEntities(dim-1); ++edge)     // global coordinates of all edges
                xs[conformingIndex(cell,edge,dim-1)] = cell.geometry().global(center(cell,edge,dim-1));
 
            // For higher order, we need to follow VTK's Lagrange cell definition.
            // See https://blog.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/
            // for some explanation.
            if (options.order > 2)
            {
              assert("VTK output of order > 2 not yet implemented\n"==0);
            }
          }
          
          for (auto const& x: xs)
            p.write(x);
        }
        else // nonconforming
        {
          // In nonconforming mode, the cells are just concatenated.
          for (auto const& cell: cells)
          {
            // First come vertices.
            for (int corner=0; corner<cell.subEntities(dim); ++corner)
              p.write(static_cast<Vector>(cell.geometry().corner(corner)));   // perform scalar type conversion on the fly
            
            // For order 2, edge midpoints come next if requested.
            if (options.order == 2)
              for (int edge=0; edge<cell.subEntities(dim-1); ++edge)
                p.write(static_cast<Vector>(cell.geometry().global(center(cell,edge,dim-1))));   // perform scalar type conversion on the fly
 
            // For higher order, we need to follow VTK's Lagrange cell definition. See
            // https://blog.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/
            // for some explanation.
            if (options.order == 3)
            {
              assert("VTK output of order > 2 not yet implemented\n"==0);
            }
          }
        }
      }
      
      void writeGridCells (int indentCount, std::ostream& s) const
      {
        // connectivity
        auto p1 = std::make_unique<VTKDataArrayWriter<int>>(s, options.outputType, "connectivity", 1,
                                                            ncorners, indentCount, options.precision);
        
        int offset = 0;
        std::vector<int> offsets; offsets.reserve(ncells);
        for (auto const& cell: cells)                                             // step through all cells
        {
          int const numCorners = cell.subEntities(dim);
          int const numEdges = cell.subEntities(dim-1);
          int const numPoints = numCorners + (options.order==2? numEdges: 0);
          
          for (int vtkPoint=0; vtkPoint<numPoints; ++vtkPoint)                   // visit each point in that cell
          {
            auto subentity = vtkPointToDune(cell.type(),vtkPoint);               // and find its Dune coordinates
            p1->write( options.dataMode==IoOptions::conforming?
                          conformingIndex(cell,subentity.first,subentity.second):
                          offset + (subentity.second==dim? 0: numCorners) + subentity.first );
          }
          offset += numPoints;                                                   // in nonconforming mode, the
          offsets.push_back(offset);                                             // points are stored blockwise,
        }                                                                        // cell by cell, just advance the offset
        p1.reset(); // write end tag on destruction
        
        // write offsets
        auto p2 = std::make_unique<VTKDataArrayWriter<int>>(s, options.outputType, "offsets", 1,
                                                            ncells, indentCount, options.precision);
        for (int off: offsets)
          p2->write(off);
        p2.reset(); // write end tag on destruction
        
        // write cell types only if dimension greater than 1 - there is only one 1D cell type: line
        if (dim>1)
        {
          VTKDataArrayWriter<unsigned char> p(s, options.outputType, "types", 1,
                                              ncells, indentCount, options.precision);
          for (auto const& cell: cells)
            p.write(static_cast<unsigned char>(vtkType(cell.type(),options.order)));
        }
      }
    
    
      GridView const& gridView;
      boost::iterator_range<CellIterator> cells;
      size_t ncorners;
      size_t npoints;
      size_t ncells;
      IoOptions const& options;
    };
    
    
    // ----------------------------------------------------------------------------------------
 
    
    template <class GridView, class Function, class Names>
    void writeCellData (VTKGridInfo<GridView> const& gridInfo, Function const& pair, Names const& names, int indentCount, std::ostream& s)
    {
      using Scalar = typename std::remove_reference_t<typename boost::fusion::result_of::value_at_c<Function,0>::type>::Scalar;
 
      auto const& f = boost::fusion::at_c<0>(pair);
      
      // Check whether this is best represented as cell data: only if it is (piecewise) constant.
      if (f.space().mapper().maxOrder() > 0)
        return;
      
      auto varDesc = boost::fusion::at_c<1>(pair);
      VTKDataArrayWriter<Scalar> p(s, gridInfo.options.outputType, names[varDesc.id], varDesc.m, gridInfo.ncells, indentCount, gridInfo.options.precision);
      for (auto const& cell: gridInfo.cells)
        p.write(f.value(cell,gridInfo.center(cell,0,0)));
    }
    
    
    
    template <class GridView, class Function, class Names>
    void writeVertexData (VTKGridInfo<GridView> const& gridInfo, Function const& pair, 
                          Names const& names, int indentCount, std::ostream& s, Timings& timer)
    {
      using ValueType = typename std::remove_reference_t<typename boost::fusion::result_of::value_at_c<Function,0>::type>::ValueType;
      using Grid = typename GridView::Grid;
      using Scalar = typename ValueType::field_type;
      constexpr int dim = GridView::dimension;
     
      auto const& f = boost::fusion::at_c<0>(pair);
      auto varDesc = boost::fusion::at_c<1>(pair);
 
      // Check whether this is best represented as vertex data: only if it is not (piecewise) constant.
      if (f.space().mapper().maxOrder() == 0)
        return;
 
      // If no name is given (empty string), we just call the function by it's variable id. Empty names are 
      // infeasible, since multiple point data shall not have the same name.
      std::string const& name = names[varDesc.id].empty() ? "var"+paddedString(varDesc.id,2)
                                                          : names[varDesc.id];
      VTKDataArrayWriter<Scalar> p(s, gridInfo.options.outputType, name,
                                   varDesc.m, gridInfo.npoints, indentCount,
                                   gridInfo.options.precision);
 
 
      if (gridInfo.options.dataMode==Kaskade::IoOptions::conforming)
      {
        std::vector<ValueType> fs(gridInfo.npoints,ValueType(0.0));   // function values, in conforming mode average over multiple points
        std::vector<short> count(gridInfo.npoints,0);                 // count how many points contribute to an entity
 
        // TODO: Currently boundary FE functions are zero on cells which only touch the boundary by one vertex (or edge in 3D). This will lead to unwanted results when values are averaged.
        // Therefore, we use the following two lines to transform boundary FE functions to continuous FE functions over the whole domain, then averaging works as expected.
        // Remove this when not necessary any longer. (And also delete ToDomainRepresetation from functionspace.hh)
        FunctionSpace_Detail::ToDomainRepresentation<std::remove_const_t<std::remove_reference_t<typename boost::fusion::result_of::value_at_c<Function,0>::type>>> tdr(f);
        auto const& df = tdr.get();
 
        // As shape function evaluation is expensive, we make use of a shape function cache
        // and the evaluator storing shape function values in it. This requires formulating
        // the evaluation in terms of a quadrature rule.
        ShapeFunctionCache<Grid,Scalar> sfCache;
        auto evaluator = f.space().evaluator(&sfCache,0);    // no derivatives needed
        CornerQuadratureRule<typename Grid::ctype,dim> qr(gridInfo.cells.begin()->type(),
                                                          gridInfo.options.order);
 
        timer.start("eval at vertices");
        for (auto const& cell: gridInfo.cells)
        {
          evaluator.moveTo(cell);
          assert(cell.type()==qr.type());
          for (int point=0; point<qr.size(); ++point)
          {
            auto [k,codim] = vtkPointToDune(cell.type(),point);
            int idx = gridInfo.conformingIndex(cell,k,codim); // global index of point
            ++count[idx];
 
            auto qp = qr[point];
            evaluator.evaluateAt(qp.position(),qr,point,0);
            fs[idx] += df.value(evaluator);
          }          
        }
        timer.stop("eval at vertices");
 
        timer.start("write out");
        for (int i=0; i<gridInfo.npoints; ++i)          // now write out the values in vertex order
          p.write(fs[i]/static_cast<typename ValueType::field_type>(count[i])); // and average them on writing
        timer.stop("write out");
      }
      else // nonconforming
      {
        for (auto const& cell: gridInfo.cells)
        {
          // First write corner values
          for (int corner=0; corner<cell.subEntities(dim); ++corner)
            p.write(f.value(cell,gridInfo.center(cell,corner,dim)));
 
          // If required, write edge midpoint values as well.
          if (gridInfo.options.order == 2)
            for (int edge=0; edge<cell.subEntities(dim-1); ++edge)
              p.write(f.value(cell,gridInfo.center(cell,edge,dim-1)));
 
          // If arbitrary higher order, write edges, faces, and interior nodes.
          // This is done recursively from outside to inside, like onion shells, and is
          // therefore done on an equidistant grid.
          // TODO: implement this. See https://blog.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/
          // TODO: If correctly implemented, this will include the order 2 case.
          if (gridInfo.options.order > 2)
          {
          }
        }
      }      
    }
    
    // The following piece of code defines the "active" scalar and vectorial variables (those used to immediately color the output in Paraview). 
    // According to VTK file formats, there need not be any active scalar or vectorial variable. Taking simply the first one as is done here makes little sense. 
    // TODO: make this configurable via IoOptions.
    template <class Functions, class Names>
    std::string guessActiveFields(Functions const& functions, Names const& names, bool cellData) 
    {
      std::string scalar = "";
      std::string vector = "";
      
      boost::fusion::for_each(functions,[&](auto const& pair) 
      {
        auto const& f = boost::fusion::at_c<0>(pair);
        auto varDesc = boost::fusion::at_c<1>(pair);
      
        // Check whether this is best represented as cell or point data.
        if ( (f.space().mapper().maxOrder()==0 && cellData)
          || (f.space().mapper().maxOrder()>0 && !cellData) )
        {
          if (f.components>1 && vector.empty())
            vector = " Vectors=\"" + names[varDesc.id] + "\"";
          if (f.components==1 && scalar.empty())
            scalar = " Scalars=\"" + names[varDesc.id] + "\"" ;
        }
      });
      return scalar+vector;
    }
  
  }
  // ---------------------------------------------------------------------------------------------
  
  template <class VariableSet>
  void writeVTK(VariableSet const& vars, IoOptions options, std::ostream& s)
  {
    using namespace VTKWriterDetail;
 
    // Measure time spent.
    auto& timer =  Timings::instance();
    timer.start("VTK output");
 
    // Standard VTK files support at most order 2. (There's a new higher order
    // type, but that's not yet supported by Kaskade 7.
    options.order = std::max(1,std::min(options.order,2));
    
    // If data mode "inferred" is requested, try to select either conforming or nonconforming
    // based on rules.
    if (options.dataMode==IoOptions::inferred)
      options.dataMode = VariableSet::Descriptions::continuity>=0 ? IoOptions::conforming
                                                                  : IoOptions::nonconforming;
    int indentCount = 0;
    constexpr int dim = VariableSet::Descriptions::Grid::dimension;
    
    timer.start("grid info");
    VTKGridInfo<typename VariableSet::Descriptions::GridView> gridInfo(vars.descriptions.gridView,options);
    timer.stop("grid info");
 
    // xml header
    s << "<?xml version=\"1.0\"?>" << std::endl;
    
    // VTKFile
    std::string const gridTag = dim>1? "UnstructuredGrid": "PolyData";
    s << "<VTKFile type=\"" << gridTag << "\" version=\"0.1\" byte_order=\"" << vtkEndianness() << "\">\n";
    ++indentCount;
    
    // UnstructuredGrid 
    indent(s,indentCount);
    s << "<" << gridTag << ">" << std::endl;
    ++indentCount;
    
    // Piece
    indent(s,indentCount);
    if (dim>1)
      s << "<Piece NumberOfPoints=\"" << gridInfo.npoints << "\" NumberOfCells=\"" << gridInfo.ncells << "\">\n";
    else
      s << "<Piece NumberOfPoints=\"" << gridInfo.npoints << "\" NumberOfVerts=\"0\" NumberOfLines=\"" << gridInfo.ncells << "\" NumberOfPolys=\"0\">\n";
    ++indentCount;
    
    // Write grid points
    timer.start("write vertices");
    indent(s,indentCount); s << "<Points>" << std::endl;
    if (options.precision <= 6)
      gridInfo.template writeGridPoints<float>(indentCount+1,s);
    else
      gridInfo.template writeGridPoints<double>(indentCount+1,s);
    indent(s,indentCount); s << "</Points>" << std::endl;
    timer.stop("write vertices");
 
    // Write grid cells
    timer.start("write cells");
    std::string const cellTagName = dim>1? "Cells": "Lines";        
    indent(s,indentCount);  s << '<' << cellTagName << ">\n";
    gridInfo.writeGridCells(indentCount+1,s);
    indent(s,indentCount); s << "</" << cellTagName << ">\n";
    timer.stop("write cells");
 
    // associate the FE functions with their variable descriptions
    auto varDesc = typename VariableSet::Descriptions::Variables();
    auto functions = boost::fusion::zip(vars.data,varDesc);
    
    // Write all functions' cell data. If a function is best represented as vertex data,
    // it will be skipped here.
    timer.start("write cell data");
    indent(s,indentCount); s << "<CellData"
                             << guessActiveFields(functions,vars.descriptions.names,true) << ">\n";
    boost::fusion::for_each(functions,[&](auto const& f)
      {
        writeCellData(gridInfo,f,vars.descriptions.names,indentCount+1,s);
      });
    indent(s,indentCount); s << "</CellData>" << std::endl;
    timer.stop("write cell data");
 
    // Write all functions' vertex data. If a function is best represented as cell data,
    // it will be skipped here.
    timer.start("write point data");
    indent(s,indentCount); s << "<PointData"
                             << guessActiveFields(functions,vars.descriptions.names,false) << ">\n";
    boost::fusion::for_each(functions,[&](auto const& f)
      {
        writeVertexData(gridInfo,f,vars.descriptions.names,indentCount+1,s,timer);
      });
    indent(s,indentCount); s << "</PointData>\n";
    timer.stop("write point data");
 
    // /Piece
    --indentCount;
    indent(s,indentCount); s << "</Piece>\n";
    
    // /UnstructuredGrid
    --indentCount;
    indent(s,indentCount); s << "</" << gridTag << ">\n";
    
    // /VTKFile
    s << "</VTKFile>" << std::endl;
 
    timer.stop("VTK output");
  }
 
  // -----------------------------------------------------------------------------------------------
  
  template <class VariableSet>
  void writeVTK(VariableSet const& vars, std::string fname, IoOptions const& options)
  {
    // generate filename for process data
    fname += (VariableSet::Descriptions::GridView::dimension>1? ".vtu" : ".vtp");
    
    // write process data
    std::ofstream file(fname);
    if (!file)
      throw Kaskade::FileIOException("Opening  failed.\n",fname,__FILE__,__LINE__);
    writeVTK(vars,options,file);
    if (!file)
      throw Kaskade::FileIOException("Writing  failed.\n",fname,__FILE__,__LINE__);
    file.close();
  }
 
}
 
// --------------------------------------------------------------------------------------------
 
 
 
#endif