backup and restore for DOFVectors

src/amdis/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:
+    /// \brief Write a hierarchic grid to disk
+    template <class GridView>
+    static void backup(GridView const& gv, std::string const& filename)
+    {
+      std::ofstream out(filename, std::ios::binary);
+      backup(gv,out);
+    }
+
+    /// \brief write a hierarchic grid into a stream
+    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 writer(gv);
+      writer.writeVertices(out);
+      writer.writeElements(out);
+    }
+
+    /// \brief read a hierarchic grid from disk
+    static Grid* restore(std::string const& filename)
+    {
+      std::ifstream in(filename, std::ios::binary);
+      return restore(in);
+    }
+
+    /// \brief read a hierarchic grid from a stream
+    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

src/amdis/ProblemStat.inc.hpp

@@ -9,6 +9,7 @@
 #include <dune/typetree/childextraction.hh>
 
 #include <amdis/AdaptInfo.hpp>
+#include <amdis/BackupRestore.hpp>
 #include <amdis/FileWriter.hpp>
 #include <amdis/Assembler.hpp>
 #include <amdis/GridFunctionOperator.hpp>
@@ -491,7 +492,12 @@ backup(AdaptInfo& adaptInfo)
   auto param = Parameters::get<std::string>(name_ + "->backup filename");
   std::string filename = param ? *param : name_ + "_" + std::to_string(adaptInfo.timestepNumber()) + ".backup";
 
+  try {
     Dune::BackupRestoreFacility<Grid>::backup(*grid_, filename);
+  } catch (Dune::NotImplemented const&) {
+    msg("Falling back to backup of gridview.");
+    BackupRestoreByGridFactory<Grid>::backup(gridView(), filename);
+  }
   // TODO: backup data
 
   msg("Problem backuped to file '{}'.", filename);
@@ -506,8 +512,13 @@ restore(Flag initFlag)
   test_exit(filesystem::exists(filename), "Restore file '{}' not found.", filename);
 
   // restore grid from file
-  std::shared_ptr<Grid> grid(Dune::BackupRestoreFacility<Grid>::restore(filename));
-  adoptGrid(grid);
+  std::shared_ptr<Grid> grid;
+  try {
+    grid.reset(Dune::BackupRestoreFacility<Grid>::restore(filename));
+  } catch (Dune::NotImplemented const&) {
+    grid.reset(BackupRestoreByGridFactory<Grid>::restore(filename));
+  }
+  adoptGrid(std::move(grid));
 
   // create fespace
   if (initFlag.isSet(INIT_FE_SPACE) || initFlag.isSet(INIT_SYSTEM))

src/amdis/linearalgebra/DOFVectorBase.hpp

@@ -189,6 +189,12 @@ namespace AMDiS
     /// Assemble all vector operators
     void assemble();
 
+    /// Write DOFVector to file
+    void backup(std::string const& filename);
+
+    /// Read backup data from file
+    void restore(std::string const& filename);
+
     /// Return the associated DataTransfer object
     std::shared_ptr<DataTransfer const> dataTransfer() const
     {

src/amdis/linearalgebra/DOFVectorBase.inc.hpp

 #pragma once
 
+#include <cstdint>
+#include <fstream>
+#include <functional>
+
 #include <amdis/Assembler.hpp>
 #include <amdis/LocalOperator.hpp>
 #include <amdis/typetree/Traversal.hpp>
@@ -89,9 +93,79 @@ assemble()
   for (auto const& element : elements(basis_->gridView())) {
     localView.bind(element);
     assemble(localView);
-    localView.unbind(element);
+    localView.unbind();
   }
   finish(true);
 }
 
+
+template <class GB, class B>
+void DOFVectorBase<GB,B>::
+backup(std::string const& filename)
+{
+  std::ofstream out(filename, std::ios::binary);
+
+  std::int64_t numElements = basis_->gridView().size(0);
+  out.write((char*)&numElements, sizeof(std::int64_t));
+
+  auto localView = basis_->localView();
+  std::vector<value_type> data;
+  for (auto const& element : elements(basis_->gridView()))
+  {
+    localView.bind(element);
+    data.clear();
+    for_each_leaf_node(localView.tree(), [&](auto const& node, auto) -> void {
+      auto const& fe = node.finiteElement();
+      std::size_t size = fe.size();
+
+      for (std::size_t i = 0; i < size; ++i)
+        data.push_back((*this)[localView.index(node.localIndex(i))]);
+    });
+
+    std::uint64_t len = data.size();
+    out.write((char*)&len, sizeof(std::uint64_t));
+    out.write((char*)data.data(), len*sizeof(value_type));
+
+    localView.unbind();
+  }
+}
+
+template <class GB, class B>
+void DOFVectorBase<GB,B>::
+restore(std::string const& filename)
+{
+  std::ifstream in(filename, std::ios::binary);
+
+  std::int64_t numElements = 0;
+  in.read((char*)&numElements, sizeof(std::int64_t));
+  assert(numElements == basis_->gridView().size(0));
+
+  // assume the order of element traversal is not changed
+  auto localView = basis_->localView();
+  std::vector<value_type> data;
+  for (auto const& element : elements(basis_->gridView()))
+  {
+    std::uint64_t len = 0;
+    in.read((char*)&len, sizeof(std::uint64_t));
+    data.resize(len);
+
+    in.read((char*)data.data(), len*sizeof(value_type));
+
+    localView.bind(element);
+    std::size_t shift = 0;
+    for_each_leaf_node(localView.tree(), [&](auto const& node, auto) -> void {
+      auto const& fe = node.finiteElement();
+      std::size_t size = fe.size();
+
+      assert(data.size() >= shift+size);
+      for (std::size_t i = 0; i < size; ++i)
+        (*this)[localView.index(node.localIndex(i))] = data[shift + i];
+
+      shift += size;
+    });
+
+    localView.unbind();
+  }
+}
+
 } // end namespace AMDiS

test/BackupRestoreTest.cpp

 #include <amdis/AMDiS.hpp>
+#include <amdis/Integrate.hpp>
 #include <amdis/ProblemStat.hpp>
 
+#if HAVE_DUNE_SPGRID
+#include <dune/grid/spgrid.hh>
+#endif
+#if HAVE_DUNE_ALUGRID
+#include <dune/alugrid/grid.hh>
+#endif
+#if HAVE_DUNE_FOAMGRID
+#include <dune/foamgrid/foamgrid.hh>
+#endif
+
 #include "Tests.hpp"
 
 using namespace AMDiS;
@@ -21,11 +32,15 @@ void test(Factory factory)
     prob.initialize(INIT_ALL);
     prob.globalRefine(2);
 
+    prob.solution(Dune::Indices::_0) << [](auto const& x) { return x; };
+    prob.solution(Dune::Indices::_1) << [](auto const& x) { return two_norm(x); };
+
     num_elements = prob.grid()->size(0);
     num_vertices = prob.grid()->size(Grid::dimension);
 
     AdaptInfo adaptInfo("adapt");
     prob.backup(adaptInfo);
+    prob.solutionVector()->backup("solution.backup");
   }
 
   // restore
@@ -35,6 +50,20 @@ void test(Factory factory)
 
     AMDIS_TEST_EQ(num_elements, prob.grid()->size(0));
     AMDIS_TEST_EQ(num_vertices, prob.grid()->size(Grid::dimension));
+
+    auto vec = *prob.solutionVector();
+    vec.restore("solution.backup");
+
+    prob.solution(Dune::Indices::_0) << [](auto const& x) { return x; };
+    prob.solution(Dune::Indices::_1) << [](auto const& x) { return two_norm(x); };
+
+    double error0 = std::sqrt(integrate(
+      unary_dot(prob.solution(Dune::Indices::_0) - makeDiscreteFunction(vec, Dune::Indices::_0)), prob.gridView()));
+    AMDIS_TEST(error0 < 1.e-10);
+
+    double error1 = std::sqrt(integrate(
+      pow<2>(prob.solution(Dune::Indices::_1) - makeDiscreteFunction(vec, Dune::Indices::_1)), prob.gridView()));
+    AMDIS_TEST(error1 < 1.e-10);
   }
 }
 
@@ -60,12 +89,26 @@ int main(int argc, char** argv)
   Parameters::set("test->backup filename", filename);
   Parameters::set("test->restore filename", filename);
 
+#if GRID_ID == 0
   test_cube<Dune::YaspGrid<2>>();
-  // test_cube<Dune::UGGrid<2>>();
-
-  // test_simplex<Dune::AlbertaGrid<2,2>>();
-  // test_simplex<Dune::UGGrid<2>>();
+#elif GRID_ID == 1 && HAVE_DUNE_UGGRID
+  test_cube<Dune::UGGrid<2>>();
+#elif GRID_ID == 2 && HAVE_DUNE_ALUGRID
+  test_cube<Dune::ALUGrid<2,2,Dune::cube,Dune::nonconforming>>();
+#elif GRID_ID == 3 && HAVE_DUNE_SPGRID
+  test<Dune::SPGrid<double,2>>([]() {
+    return std::unique_ptr<Dune::SPGrid<double,2>>(
+      new Dune::SPGrid<double,2>(Dune::SPDomain<double, 2>({0.0,0.0}, {1.0,1.0}), Dune::SPMultiIndex<2>({2,2}))); });
+#elif GRID_ID == 4 && HAVE_DUNE_UGGRID
+  test_simplex<Dune::UGGrid<2>>();
+#elif GRID_ID == 5 && HAVE_DUNE_ALUGRID
+  test_simplex<Dune::ALUGrid<2,2,Dune::simplex,Dune::nonconforming>>();
+#elif GRID_ID == 6 && HAVE_DUNE_ALUGRID
+  test_simplex<Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming>>();
+#endif
+  // test_simplex<Dune::AlbertaGrid<2,2>>(); // Segmentation fault
+  // test_simplex<Dune::FoamGrid<2,2>>(); // Segmentation fault
 
   AMDiS::finalize();
-  return 0;
+  return report_errors();
 }

test/CMakeLists.txt

@@ -2,9 +2,15 @@
 dune_add_test(SOURCES AdaptInfoTest.cpp
   LINK_LIBRARIES amdis)
 
-dune_add_test(SOURCES BackupRestoreTest.cpp
-  LINK_LIBRARIES amdis)
-add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 BackupRestoreTest)
+foreach(_GRID RANGE 6)
+  dune_add_test(NAME "BackupRestoreTest_${_GRID}"
+    SOURCES BackupRestoreTest.cpp
+    COMPILE_DEFINITIONS "GRID_ID=${_GRID}"
+    LABELS "BackupRestore"
+    LINK_LIBRARIES amdis)
+  add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "BackupRestoreTest_${_GRID}")
+endforeach()
+unset(_GRID)
 
 dune_add_test(SOURCES ConceptsTest.cpp
   LINK_LIBRARIES amdis)