Merge branch 'feature/backup_restore' into 'master'

Backup and restore facilities

See merge request !38

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<GridView> 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/Initfile.hpp

@@ -47,6 +47,12 @@ namespace AMDiS
       value = pt().get(key, value);
     }
 
+    template <class T>
+    static void set(std::string const& key, T const& value)
+    {
+      pt()[key] = value;
+    }
+
     /// Returns specified info level
     static int getMsgInfo()
     {

src/amdis/Mesh.hpp

@@ -119,14 +119,19 @@ namespace AMDiS
         Dune::GmshReader<Grid> reader;
         return std::unique_ptr<Grid>{reader.read(filename)};
       }
-#if HAVE_ALBERTA
-      else if (ext == ".1d" || ext == ".2d" || ext == ".3d") {
-        Dune::GridFactory<Grid> factory;
-        Dune::AlbertaReader<Grid> reader;
-        reader.readGrid(filename, factory);
-        return std::unique_ptr<Grid>{factory.createGrid()};
-      }
-#endif
+// #if HAVE_ALBERTA
+//       else if (ext == ".1d" || ext == ".2d" || ext == ".3d") {
+//         Dune::Hybrid::ifElse(bool_t<ALBERTA_DIM == Grid::dimensionworld>{},
+//         [&](auto id) {
+//           Dune::GridFactory<Grid> factory;
+//           Dune::AlbertaReader<Grid> reader;
+//           id(reader).readGrid(filename, id(factory));
+//           return std::unique_ptr<Grid>{factory.createGrid()};
+//         });
+
+//         warning("WORLDDIM must be given in Alberta flags.");
+//       }
+// #endif
 
       error_exit("No way to construct UG-Grid found");
       return {};

src/amdis/PeriodicBC.inc.hpp

@@ -84,7 +84,7 @@ initAssociations(Basis const& basis)
         auto x = faceTrafo_.evaluate(insideCoords[i]);
         for (std::size_t j = 0; j < outsideCoords.size(); ++j) {
           auto const& y = outsideCoords[j];
-          if ((x-y).two_norm() < tol) {
+          if (distance(x,y) < tol) {
             periodicNodes_[insideGlobalDOFs[i]] = true;
             associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]);
           }
@@ -183,7 +183,7 @@ initAssociations2(Basis const& basis)
       auto x = faceTrafo_.evaluate(insideCoords[i]);
       for (std::size_t j = 0; j < outsideCoords.size(); ++j) {
         auto const& y = outsideCoords[j];
-        if ((x-y).two_norm() < tol) {
+        if (distance(x,y) < tol) {
           periodicNodes_[insideGlobalDOFs[i]] = true;
           associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]);
         }

src/amdis/ProblemStat.hpp

@@ -295,6 +295,18 @@ namespace AMDiS
     /// Writes output files.
     void writeFiles(AdaptInfo& adaptInfo, bool force = false);
 
+    /// Backup the grid
+    void backup(std::string const& filename) const;
+
+    /// Implements \ref ProblemStatBase::backup
+    void backup(AdaptInfo& adaptInfo) override;
+
+    /// Retore the grid
+    auto restore(std::string const& filename);
+
+    /// Implements \ref ProblemStatBase::restore
+    void restore(Flag initFlag) override;
+
 
   public: // get-methods
 

src/amdis/ProblemStat.inc.hpp

@@ -6,9 +6,11 @@
 
 #include <dune/common/hybridutilities.hh>
 #include <dune/common/timer.hh>
+#include <dune/grid/common/capabilities.hh>
 #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>
@@ -477,4 +479,87 @@ writeFiles(AdaptInfo& adaptInfo, bool force)
 }
 
 
+template <class Traits>
+void ProblemStat<Traits>::
+backup(std::string const& filename) const
+{
+  if (Dune::Capabilities::hasBackupRestoreFacilities<Grid>::v)
+    Dune::BackupRestoreFacility<Grid>::backup(*grid_, filename);
+  else {
+    warning("Falling back to backup of gridview.");
+    BackupRestoreByGridFactory<Grid>::backup(gridView(), filename);
+  }
+}
+
+
+template <class Traits>
+void ProblemStat<Traits>::
+backup(AdaptInfo& adaptInfo)
+{
+  auto param = Parameters::get<std::string>(name_ + "->backup->grid");
+  std::string grid_filename = param ? *param :
+    name_ + "_" + std::to_string(adaptInfo.timestepNumber()) + ".grid";
+
+  auto param2 = Parameters::get<std::string>(name_ + "->backup->solution");
+  std::string solution_filename = param2 ? *param2 :
+    name_ + "_" + std::to_string(adaptInfo.timestepNumber()) + ".solution";
+
+  backup(grid_filename);
+  solution_->backup(solution_filename);
+  msg("Problem backuped to files '{}' and '{}'.", grid_filename, solution_filename);
+}
+
+
+template <class Traits>
+auto ProblemStat<Traits>::
+restore(std::string const& filename)
+{
+  // restore grid from file
+  std::unique_ptr<Grid> grid;
+  if (Dune::Capabilities::hasBackupRestoreFacilities<Grid>::v)
+    grid.reset(Dune::BackupRestoreFacility<Grid>::restore(filename));
+  else
+    grid.reset(BackupRestoreByGridFactory<Grid>::restore(filename));
+  return std::move(grid);
+}
+
+
+template <class Traits>
+void ProblemStat<Traits>::
+restore(Flag initFlag)
+{
+  std::string grid_filename = Parameters::get<std::string>(name_ + "->restore->grid").value();
+  std::string solution_filename = Parameters::get<std::string>(name_ + "->restore->solution").value();
+  test_exit(filesystem::exists(grid_filename), "Restore file '{}' not found.", grid_filename);
+  test_exit(filesystem::exists(solution_filename), "Restore file '{}' not found.", solution_filename);
+
+  // restore grid from file
+  adoptGrid(restore(grid_filename));
+
+  // create fespace
+  if (initFlag.isSet(INIT_FE_SPACE) || initFlag.isSet(INIT_SYSTEM))
+    createGlobalBasis();
+
+  // create system
+  if (initFlag.isSet(INIT_SYSTEM))
+    createMatricesAndVectors();
+
+  // create solver
+  if (linearSolver_)
+    warning("solver already created\n");
+  else if (initFlag.isSet(INIT_SOLVER))
+    createSolver();
+
+  // create marker
+  if (initFlag.isSet(INIT_MARKER))
+    createMarker();
+
+  // create file writer
+  if (initFlag.isSet(INIT_FILEWRITER))
+    createFileWriter();
+
+  solution_->compress();
+  solution_->restore(solution_filename);
+}
+
 } // end namespace AMDiS

src/amdis/ProblemStatBase.hpp

@@ -111,6 +111,12 @@ namespace AMDiS
      */
     virtual void estimate(AdaptInfo& adaptInfo) = 0;
 
+    /// \brief backup the grid and the solution to a file
+    virtual void backup(AdaptInfo& adaptInfo) = 0;
+
+    /// \brief restore the grid and the solution from a file
+    virtual void restore(Flag initFlag) = 0;
+
     /// Returns the name of the problem.
     virtual std::string const& name() const = 0;
   };

src/amdis/common/FieldMatVec.hpp

@@ -271,6 +271,10 @@ namespace Dune
 
 
   /// Dot-product with the vector itself
+  template <class T,
+    std::enable_if_t<Dune::IsNumber<T>::value, int> = 0>
+  auto unary_dot(T const& x);
+
   template <class T, int N>
   auto unary_dot(FieldVector<T, N> const& x);
 
@@ -319,6 +323,10 @@ namespace Dune
   /** \ingroup vector_norms
    *  \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 )
    **/
+  template <class T,
+    std::enable_if_t<Dune::IsNumber<T>::value, int> = 0>
+  auto two_norm(T const& x);
+
   template <class T, int N>
   auto two_norm(FieldVector<T, N> const& x);
 
@@ -346,6 +354,10 @@ namespace Dune
   // ----------------------------------------------------------------------------
 
   /// The euklidean distance between two vectors = |lhs-rhs|_2
+  template <class T,
+    std::enable_if_t<Dune::IsNumber<T>::value, int> = 0>
+  T distance(T const& lhs, T const& rhs);
+
   template <class T, int N>
   T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs);
 

src/amdis/common/FieldMatVec.inc.hpp

@@ -239,6 +239,14 @@ T sum(FieldMatrix<T, 1, N> const& x)
 
 
 /// Dot-product with the vector itself
+template <class T,
+  std::enable_if_t<Dune::IsNumber<T>::value, int> >
+auto unary_dot(T const& x)
+{
+  using std::abs;
+  return AMDiS::Math::sqr(abs(x));
+}
+
 template <class T, int N>
 auto unary_dot(FieldVector<T, N> const& x)
 {
@@ -327,6 +335,14 @@ auto one_norm(FieldMatrix<T, 1, N> const& x)
 /** \ingroup vector_norms
   *  \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 )
   **/
+template <class T,
+  std::enable_if_t<Dune::IsNumber<T>::value, int> >
+auto two_norm(T const& x)
+{
+  using std::abs;
+  return abs(x);
+}
+
 template <class T, int N>
 auto two_norm(FieldVector<T, N> const& x)
 {
@@ -374,6 +390,14 @@ auto infty_norm(FieldMatrix<T, 1, N> const& x)
 // ----------------------------------------------------------------------------
 
 /// The euklidean distance between two vectors = |lhs-rhs|_2
+template <class T,
+  std::enable_if_t<Dune::IsNumber<T>::value, int> >
+T distance(T const& lhs, T const& rhs)
+{
+  using std::abs;
+  return abs(lhs - rhs);
+}
+
 template <class T, int N>
 T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs)
 {

src/amdis/linearalgebra/DOFVectorBase.hpp

@@ -226,6 +226,12 @@ namespace AMDiS
     /// Assemble all vector operators added by \ref addOperator().
     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>
@@ -91,9 +95,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

src/amdis/operations/FieldMatVec.hpp

@@ -47,6 +47,7 @@ namespace AMDiS
       template <class V>
       constexpr auto operator()(V const& vec) const
       {
+        using Dune::unary_dot;
         return unary_dot(vec);
       }
 
@@ -74,6 +75,7 @@ namespace AMDiS
       template <class V>
       constexpr auto operator()(V const& vec) const
       {
+        using Dune::two_norm;
         return two_norm(vec);
       }
     };
@@ -85,6 +87,7 @@ namespace AMDiS
       template <class M>
       constexpr auto operator()(M const& mat) const
       {
+        using Dune::trans;
         return trans(mat);
       }
 

test/BackupRestoreTest.cpp

+#include <amdis/AMDiS.hpp>
+#include <amdis/Integrate.hpp>
+#include <amdis/ProblemStat.hpp>
+
+#include <dune/grid/onedgrid.hh>
+#include <dune/grid/yaspgrid.hh>
+
+#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;
+
+template <class Grid, class Factory>
+void test(Factory factory)
+{
+  using Param  = TaylorHoodBasis<Grid>;
+