backup and restore of grid and dofvector

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

@@ -83,7 +83,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]);
           }
@@ -180,7 +180,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

@@ -274,9 +274,15 @@ 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;
 

src/amdis/ProblemStat.inc.hpp

@@ -487,38 +487,57 @@ writeFiles(AdaptInfo& adaptInfo, bool force)
 
 template <class Traits>
 void ProblemStat<Traits>::
-backup(AdaptInfo& adaptInfo)
+backup(std::string const& filename) const
 {
-  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.");
+    warning("Falling back to backup of gridview.");
     BackupRestoreByGridFactory<Grid>::backup(gridView(), filename);
   }
-  // TODO: backup data
-
-  msg("Problem backuped to file '{}'.", filename);
 }
 
-
 template <class Traits>
 void ProblemStat<Traits>::
-restore(Flag initFlag)
+backup(AdaptInfo& adaptInfo)
 {
-  std::string filename = Parameters::get<std::string>(name_ + "->restore filename").value();
-  test_exit(filesystem::exists(filename), "Restore file '{}' not found.", filename);
+  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::shared_ptr<Grid> grid;
+  std::unique_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));
+  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))
@@ -543,8 +562,7 @@ restore(Flag initFlag)
     createFileWriter();
 
   solution_->compress();
-
-  // TODO: read data
+  solution_->restore(solution_filename);
 }
 
 } // end namespace AMDiS

src/amdis/common/FieldMatVec.hpp

@@ -255,6 +255,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);
 
@@ -303,6 +307,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);
 
@@ -330,6 +338,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/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

@@ -2,6 +2,9 @@
 #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
@@ -40,7 +43,6 @@ void test(Factory factory)
 
     AdaptInfo adaptInfo("adapt");
     prob.backup(adaptInfo);
-    prob.solutionVector()->backup("solution.backup");
   }
 
   // restore
@@ -51,11 +53,9 @@ 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); };
+    DOFVector<typename Param::GlobalBasis> vec(*prob.globalBasis());
+    makeDOFVectorView(vec, Dune::Indices::_0) << [](auto const& x) { return x; };
+    makeDOFVectorView(vec, 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()));
@@ -71,40 +71,62 @@ template <class Grid>
 void test_cube()
 {
   using Factory = Dune::StructuredGridFactory<Grid>;
-  test<Grid>([]() { return Factory::createCubeGrid({0.0,0.0}, {1.0,1.0}, std::array<unsigned int,2>{2,2}); });
+  Dune::FieldVector<double,Grid::dimensionworld> left(0.0);
+  Dune::FieldVector<double,Grid::dimensionworld> right(1.0);
+  auto sizes = Dune::filledArray<Grid::dimension, unsigned int>(2);
+  test<Grid>([&]() { return Factory::createCubeGrid(left, right, sizes); });
 }
 
 template <class Grid>
 void test_simplex()
 {
   using Factory = Dune::StructuredGridFactory<Grid>;
-  test<Grid>([]() { return Factory::createSimplexGrid({0.0,0.0}, {1.0,1.0}, std::array<unsigned int,2>{2,2}); });
+  Dune::FieldVector<double,Grid::dimensionworld> left(0.0);
+  Dune::FieldVector<double,Grid::dimensionworld> right(1.0);
+  auto sizes = Dune::filledArray<Grid::dimension, unsigned int>(2);
+  test<Grid>([&]() { return Factory::createSimplexGrid(left, right, sizes); });
 }
 
 int main(int argc, char** argv)
 {
   AMDiS::init(argc, argv);
 
-  std::string filename = "test.backup";
-  Parameters::set("test->backup filename", filename);
-  Parameters::set("test->restore filename", filename);
+  std::string filename = "test";
+  Parameters::set("test->backup->grid", filename + ".grid");
+  Parameters::set("test->restore->grid", filename + ".grid");
+  Parameters::set("test->backup->solution", filename + ".solution");
+  Parameters::set("test->restore->solution", filename + ".solution");
 
 #if GRID_ID == 0
   test_cube<Dune::YaspGrid<2>>();
+  test_cube<Dune::YaspGrid<3>>();
 #elif GRID_ID == 1 && HAVE_DUNE_UGGRID
   test_cube<Dune::UGGrid<2>>();
+  test_cube<Dune::UGGrid<3>>();
 #elif GRID_ID == 2 && HAVE_DUNE_ALUGRID
   test_cube<Dune::ALUGrid<2,2,Dune::cube,Dune::nonconforming>>();
+  test_cube<Dune::ALUGrid<3,3,Dune::cube,Dune::nonconforming>>();
+  //test_cube<Dune::ALUGrid<2,3,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}))); });
+  test<Dune::SPGrid<double,3>>([]() {
+    return std::unique_ptr<Dune::SPGrid<double,3>>(
+      new Dune::SPGrid<double,3>(Dune::SPDomain<double, 3>({0.0,0.0,0.0}, {1.0,1.0,1.0}), Dune::SPMultiIndex<3>({2,2,2}))); });
 #elif GRID_ID == 4 && HAVE_DUNE_UGGRID
   test_simplex<Dune::UGGrid<2>>();
+  test_simplex<Dune::UGGrid<3>>();
 #elif GRID_ID == 5 && HAVE_DUNE_ALUGRID
   test_simplex<Dune::ALUGrid<2,2,Dune::simplex,Dune::nonconforming>>();
+  test_simplex<Dune::ALUGrid<3,3,Dune::simplex,Dune::nonconforming>>();
+  //test_simplex<Dune::ALUGrid<2,3,Dune::simplex,Dune::nonconforming>>();
 #elif GRID_ID == 6 && HAVE_DUNE_ALUGRID
   test_simplex<Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming>>();
+  test_simplex<Dune::ALUGrid<3,3,Dune::simplex,Dune::conforming>>();
+  //test_simplex<Dune::ALUGrid<2,3,Dune::simplex,Dune::conforming>>();
+#elif GRID_ID == 7
+  test_cube<Dune::OneDGrid>();
 #endif
   // test_simplex<Dune::AlbertaGrid<2,2>>(); // Segmentation fault
   // test_simplex<Dune::FoamGrid<2,2>>(); // Segmentation fault

test/CMakeLists.txt

@@ -2,7 +2,7 @@
 dune_add_test(SOURCES AdaptInfoTest.cpp
   LINK_LIBRARIES amdis)
 
-foreach(_GRID RANGE 6)
+foreach(_GRID RANGE 7)
   dune_add_test(NAME "BackupRestoreTest_${_GRID}"
     SOURCES BackupRestoreTest.cpp
     COMPILE_DEFINITIONS "GRID_ID=${_GRID}"