Restructuring of Matrix and Vector and linear-algebra backends

Summary

Restructuring of Matrix and Vector and linear-algebra backends

Details

• Restructuring the vector and matrix backend wrappers
• Added init-finalize also for vectors
• Require gather/scatter operations instead of direct vector/matrix access
• Rewritten the interpolation
• Added ParallelIndexSet and DOFMapping to be used also in other backends

examples/ellipt.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-#include <iostream>
-
 #include <amdis/AMDiS.hpp>
 #include <amdis/Integrate.hpp>
 #include <amdis/LocalOperators.hpp>
 #include <amdis/ProblemStat.hpp>
 #include <amdis/common/Literals.hpp>
 
+#include <dune/grid/yaspgrid.hh>
+using Grid = Dune::YaspGrid<GRIDDIM>;
+
 using namespace AMDiS;
 using namespace Dune::Indices;
 
-// 1 component with polynomial degree 1
-using Param   = YaspGridBasis<GRIDDIM, 2>;
-using ElliptProblem = ProblemStat<Param>;
-
 int main(int argc, char** argv)
 {
   Environment env(argc, argv);
+  Dune::Timer t;
 
-  int numLevels = GRIDDIM == 2 ? 8 : 5;
+  int numLevels = GRIDDIM == 2 ? 6 : 4;
   if (argc > 2)
     numLevels = std::atoi(argv[2]);
 
@@ -34,7 +29,8 @@ int main(int argc, char** argv)
     return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
   };
 
-  ElliptProblem prob("ellipt");
+  using Param = LagrangeBasis<Grid, 2>;
+  ProblemStat<Param> prob("ellipt");
   prob.initialize(INIT_ALL);
 
   auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
@@ -44,40 +40,44 @@ int main(int argc, char** argv)
   prob.addVectorOperator(opForce, 0);
 
   // set boundary condition
-  prob.addDirichletBC(BoundaryType{1}, 0, 0, g);
+  prob.addDirichletBC(1, 0, 0, g);
 
   AdaptInfo adaptInfo("adapt");
 
   std::vector<double> errL2; errL2.reserve(numLevels);
   std::vector<double> errH1; errH1.reserve(numLevels);
   std::vector<double> widths; widths.reserve(numLevels);
-  for (int i = 0; i < numLevels; ++i) {
-    prob.grid()->globalRefine(1);
-    auto gridView = prob.gridView();
+  for (int l = 0; l < numLevels; ++l) {
+    prob.globalRefine(1);
 
     double h = 0;
-    for (auto const& e : edges(gridView))
-      h = std::max(h, e.geometry().volume());
+    for (auto const& e : elements(prob.gridView(), Dune::Partitions::interior)) {
+      auto refElem = Dune::referenceElement<double,GRIDDIM>(e.type());
+      auto geo = e.geometry();
+      for (int i = 0; i < refElem.size(GRIDDIM-1); ++i) { // edges
+        auto v0 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,0,GRIDDIM), GRIDDIM));
+        auto v1 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,1,GRIDDIM), GRIDDIM));
+        h = std::max(h, (v0 - v1).two_norm());
+      }
+    }
     widths.push_back(h);
 
-    prob.globalBasis()->update(gridView);
-    prob.solutionVector()->compress();
     prob.assemble(adaptInfo);
     prob.solve(adaptInfo);
 
-    double errorL2 = integrate(sqr(g - prob.solution(0)), gridView, 6);
+    double errorL2 = integrate(sqr(g - prob.solution(0)), prob.gridView(), 6);
     errL2.push_back(std::sqrt(errorL2));
-    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), gridView, 6);
+    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), prob.gridView(), 6);
     errH1.push_back(std::sqrt(errorH1));
 
 #if WRITE_FILES
     Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(gridView);
     vtkWriter.addVertexData(prob.solution(0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
-    vtkWriter.write("u_" + std::to_string(i));
+    vtkWriter.write("u_" + std::to_string(l));
 #endif
   }
 
-  std::cout << "\n";
+  msg("");
   msg("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}",
       "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1");
   msg_("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n");
@@ -92,5 +92,6 @@ int main(int argc, char** argv)
     msg("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}",
         i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]);
 
+  msg("elapsed time: {} seconds", t.elapsed());
   return 0;
 }