BackupRestore.hpp

#pragma once
 
#include <cstdint>
#include <fstream>
#include <memory>
#include <vector>
 
#include <dune/grid/common/backuprestore.hh>
#include <dune/grid/common/gridfactory.hh>
 
namespace AMDiS
{
  // Fallback Implementation for Grids not supporting direct BackupRestore
  template <class Grid>
  class BackupRestoreByGridFactory
  {
    template <class GridView>
    class Writer
    {
      enum { dim = GridView::dimension };
      enum { dow = GridView::dimensionworld };
 
      using ct = typename GridView::ctype;
 
    public:
      Writer(GridView const& gv)
        : gridView_(gv)
      {
        std::int64_t num = 0;
        indexMap_.resize(gridView_.size(dim));
        auto const& indexSet = gridView_.indexSet();
        for (auto const& vertex : vertices(gridView_))
          indexMap_[indexSet.index(vertex)] = std::int64_t(num++);
      }
 
      void writeVertices(std::ostream& out) const
      {
        for (auto const& vertex : vertices(gridView_)) {
          auto v = vertex.geometry().center();
          out.write((char*)&v, dow*sizeof(ct));
        }
      }
 
      void writeElements(std::ostream& out) const
      {
        auto const& indexSet = gridView_.indexSet();
 
        std::vector<std::int64_t> connectivity;
        connectivity.reserve(8);
        for (auto const& e : elements(gridView_)) {
          unsigned int id = e.type().id();
          out.write((char*)&id, sizeof(unsigned int));
 
          connectivity.clear();
          for (unsigned int j = 0; j < e.subEntities(dim); ++j)
            connectivity.emplace_back(indexMap_[indexSet.subIndex(e,j,dim)]);
 
          out.write((char*)connectivity.data(), connectivity.size()*sizeof(std::int64_t));
        }
      }
 
    private:
      GridView gridView_;
      std::vector<std::int64_t> indexMap_;
    };
 
    class Reader
    {
      enum { dim = Grid::dimension };
      enum { dow = Grid::dimensionworld };
 
      using ct = typename Grid::ctype;
      using Factory = Dune::GridFactory<Grid>;
      using GlobalCoordinates = Dune::FieldVector<ct,dow>;
 
    public:
      Reader(Factory& factory, std::istream& in)
        : factory_(factory)
      {
        in.read((char*)&numElements_, sizeof(std::int64_t));
        in.read((char*)&numVertices_, sizeof(std::int64_t));
      }
 
      void readVertices(std::istream& in) const
      {
        GlobalCoordinates p;
        for (std::int64_t i = 0; i < numVertices_; ++i) {
          in.read((char*)&p[0], dow*sizeof(ct));
          factory_.insertVertex(p);
        }
      }
 
      void readElements(std::istream& in) const
      {
        std::vector<std::int64_t> connectivity(8);
        std::vector<unsigned int> vertices;
        vertices.reserve(8);
 
        for (std::int64_t i = 0; i < numElements_; ++i) {
          unsigned int id = 0;
          in.read((char*)&id, sizeof(unsigned int));
 
          Dune::GeometryType type(id,dim);
          auto refElem = Dune::referenceElement<ct,dim>(type);
 
          in.read((char*)connectivity.data(), refElem.size(dim)*sizeof(std::int64_t));
          vertices.clear();
          std::copy_n(connectivity.begin(),refElem.size(dim),std::back_inserter(vertices));
 
          factory_.insertElement(type, vertices);
        }
      }
 
    private:
      Factory& factory_;
      std::int64_t numElements_ = 0;
      std::int64_t numVertices_ = 0;
    };
 
  public:
    template <class GridView>
    static void backup(GridView const& gv, std::string const& filename)
    {
      std::ofstream out(filename, std::ios::binary);
      backup(gv,out);
    }
 
    template <class GridView>
    static void backup(GridView const& gv, std::ostream& out)
    {
      std::int32_t dim = GridView::dimension;
      std::int32_t dow = GridView::dimensionworld;
      std::int64_t num_elements = gv.size(0);
      std::int64_t num_vertices = gv.size(dim);
 
      out.write((char*)&dim, sizeof(std::int32_t));
      out.write((char*)&dow, sizeof(std::int32_t));
      out.write((char*)&num_elements, sizeof(std::int64_t));
      out.write((char*)&num_vertices, sizeof(std::int64_t));
 
      Writer<GridView> writer(gv);
      writer.writeVertices(out);
      writer.writeElements(out);
    }
 
    static Grid* restore(std::string const& filename)
    {
      std::ifstream in(filename, std::ios::binary);
      return restore(in);
    }
 
    static Grid* restore(std::istream& in)
    {
      Dune::GridFactory<Grid> factory;
 
      std::int32_t dim = 0, dow = 0;
 
      in.read((char*)&dim, sizeof(std::int32_t));
      in.read((char*)&dow, sizeof(std::int32_t));
 
      assert(Grid::dimension == dim);
      assert(Grid::dimensionworld == dow);
 
      Reader reader(factory, in);
      reader.readVertices(in);
      reader.readElements(in);
 
      std::unique_ptr<Grid> ptr(factory.createGrid());
      return ptr.release();
    }
  };
 
} // end namespace AMDiS