Add SymmetryStructure includes to MatrixBackends and remove changes from ellipt demo

examples/CMakeLists.txt

@@ -31,23 +31,4 @@ add_dependencies(examples
 if (dune-foamgrid_FOUND)
   add_amdis_executable(NAME surface.2d SOURCES surface.cc DIM 2 DOW 3 ALBERTA_GRID)
   add_dependencies(examples surface.2d)
-endif ()
-
-add_custom_target(ellipt)
-foreach(_D 2 3)
-  add_amdis_executable(NAME ellipt.yaspgrid.${_D}d SOURCES ellipt.cc DIM ${_D} DOW ${_D})
-  target_compile_definitions(ellipt.yaspgrid.${_D}d PRIVATE "GRID_ID=0")
-  add_dependencies(ellipt ellipt.yaspgrid.${_D}d)
-
-  add_amdis_executable(NAME ellipt.albertagrid.${_D}d SOURCES ellipt.cc DIM ${_D} DOW ${_D} ALBERTA_GRID)
-  target_compile_definitions(ellipt.albertagrid.${_D}d PRIVATE "GRID_ID=1")
-  add_dependencies(ellipt ellipt.albertagrid.${_D}d)
-
-  add_amdis_executable(NAME ellipt.uggrid.${_D}d SOURCES ellipt.cc DIM ${_D} DOW ${_D})
-  target_compile_definitions(ellipt.uggrid.${_D}d PRIVATE "GRID_ID=2")
-  add_dependencies(ellipt ellipt.uggrid.${_D}d)
-
-  add_amdis_executable(NAME ellipt.alugrid.${_D}d SOURCES ellipt.cc DIM ${_D} DOW ${_D})
-  target_compile_definitions(ellipt.alugrid.${_D}d PRIVATE "GRID_ID=3")
-  add_dependencies(ellipt ellipt.alugrid.${_D}d)
-endforeach()
+endif ()
\ No newline at end of file

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>
 
-#if GRID_ID == 1
-#include <dune/grid/albertagrid.hh>
-using Grid = Dune::AlbertaGrid<GRIDDIM,WORLDDIM>;
-#elif GRID_ID == 2
-#include <dune/grid/uggrid.hh>
-using Grid = Dune::UGGrid<GRIDDIM>;
-#elif GRID_ID == 3
-#include <dune/alugrid/grid.hh>
-using Grid = Dune::ALUGrid<GRIDDIM,WORLDDIM,Dune::simplex,Dune::conforming>;
-#else
-#include <dune/grid/yaspgrid.hh>
-using Grid = Dune::YaspGrid<GRIDDIM>;
-#endif
-
 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 ? 6 : 4;
+  int numLevels = GRIDDIM == 2 ? 8 : 5;
   if (argc > 2)
     numLevels = std::atoi(argv[2]);
 
@@ -40,8 +34,7 @@ int main(int argc, char** argv)
     return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
   };
 
-  using Param = LagrangeBasis<Grid, 2>;
-  ProblemStat<Param> prob("ellipt");
+  ElliptProblem prob("ellipt");
   prob.initialize(INIT_ALL);
 
   auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
@@ -51,49 +44,40 @@ int main(int argc, char** argv)
   prob.addVectorOperator(opForce, 0);
 
   // set boundary condition
-#if GRIDDIM == 2
-  prob.boundaryManager()->setBoxBoundary({1,1,1,1});
-#elif GRIDDIM == 3
-  prob.boundaryManager()->setBoxBoundary({1,1,1,1,1,1});
-#endif
-  prob.addDirichletBC(1, 0, 0, g);
+  prob.addDirichletBC(BoundaryType{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 l = 0; l < numLevels; ++l) {
-    prob.globalRefine(1);
+  for (int i = 0; i < numLevels; ++i) {
+    prob.grid()->globalRefine(1);
+    auto gridView = prob.gridView();
 
     double h = 0;
-    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());
-      }
-    }
+    for (auto const& e : edges(gridView))
+      h = std::max(h, e.geometry().volume());
     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)), prob.gridView(), 6);
+    double errorL2 = integrate(sqr(g - prob.solution(0)), gridView, 6);
     errL2.push_back(std::sqrt(errorL2));
-    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), prob.gridView(), 6);
+    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), 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(l));
+    vtkWriter.write("u_" + std::to_string(i));
 #endif
   }
 
-  msg("");
+  std::cout << "\n";
   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");
@@ -108,6 +92,5 @@ 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;
 }

examples/init/boundary.dat.2d

-elliptMesh->global refinements: 3
+mesh->global refinements: 3
 
-ellipt->mesh:                   elliptMesh
+ellipt->mesh:                   mesh
 
 ellipt->solver:                 bicgstab
 ellipt->solver->max iteration:  10000

examples/init/cahn_hilliard.dat.2d

@@ -14,7 +14,7 @@ refinement->bulk: 2
 ch->mesh:                       chMesh
 ch->solver:                     direct
 ch->solver->max iteration:      1000
-ch->solver->absolute tolerance: 1e-6
+ch->solver->relative tolerance: 1e-6
 ch->solver->break if tolerance not reached: 1
 ch->solver->info: 2
 

examples/init/ellipt.dat.2d

 elliptMesh->global refinements: 0
 elliptMesh->num cells: 8 8
 elliptMesh->overlap: 1
-% elliptMesh->structured: simplex
 
 ellipt->mesh: elliptMesh
 ellipt->symmetry: spd

examples/init/ellipt.dat.3d

 elliptMesh->global refinements: 0
 elliptMesh->num cells: 8 8 8
 elliptMesh->overlap: 1
-% elliptMesh->structured: simplex
 
 ellipt->mesh: elliptMesh
 ellipt->symmetry: spd

examples/init/heat.dat.2d

@@ -5,6 +5,7 @@ heatMesh->num cells: 8 8
 
 heat->mesh: heatMesh
 heat->names: u
+
 heat->solver: pcg
 heat->solver->max iteration: 1000
 heat->solver->absolute tolerance: 1e-6

src/amdis/BiLinearForm.hpp

@@ -24,7 +24,8 @@ namespace AMDiS
       : public MatrixBase<RB,CB,MatrixBackend<BackendTraits<RB,T>>>
   {
     using Self = BiLinearForm;
-    using Backend = MatrixBackend<BackendTraits<RB,T>>;
+    using Traits = BackendTraits<RB,T>;
+    using Backend = MatrixBackend<Traits>;
     using Super = MatrixBase<RB,CB,Backend>;
 
   public:

src/amdis/linearalgebra/MatrixBase.hpp

@@ -8,6 +8,7 @@
 #include <amdis/common/Concepts.hpp>
 #include <amdis/common/TypeTraits.hpp>
 #include <amdis/functions/NodeIndices.hpp>
+#include <amdis/linearalgebra/SymmetryStructure.hpp>
 #include <amdis/typetree/MultiIndex.hpp>
 
 namespace AMDiS

src/amdis/linearalgebra/eigen/MatrixBackend.hpp

@@ -8,6 +8,7 @@
 #include <dune/common/timer.hh>
 
 #include <amdis/Output.hpp>
+#include <amdis/linearalgebra/SymmetryStructure.hpp>
 
 namespace AMDiS
 {

src/amdis/linearalgebra/istl/MatrixBackend.hpp

@@ -7,6 +7,7 @@
 #include <dune/istl/matrixindexset.hh>
 
 #include <amdis/Output.hpp>
+#include <amdis/linearalgebra/SymmetryStructure.hpp>
 
 namespace AMDiS
 {

src/amdis/linearalgebra/mtl/MatrixBackend.hpp

@@ -10,6 +10,7 @@
 #include <boost/numeric/mtl/utility/range_wrapper.hpp>
 
 #include <amdis/Output.hpp>
+#include <amdis/linearalgebra/SymmetryStructure.hpp>
 
 namespace AMDiS
 {