Add sections in generated drivers for much better reading

The MR is best described with an example of generated code:

bool driver(int argc, char** argv){  
  // Initialize basic stuff...  
  using RangeType = double;
  Dune::ParameterTree initree;
  Dune::ParameterTreeParser::readINITree(argv[1], initree);
  
  // Setup grid (view)...  
  using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
  using GV = Grid::LeafGridView;
  using DF = Grid::ctype;
  IniGridFactory<Grid> factory(initree);
  std::shared_ptr<Grid> grid = factory.getGrid();
  GV gv = grid->leafGridView();
  
  // Set up finite element maps...  
  using P1_FEM = Dune::PDELab::PkLocalFiniteElementMap<GV, DF, RangeType, 1>;
  P1_FEM p1_fem(gv);
  
  // Set up grid function spaces...  
  using VectorBackendP1 = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::none>;
  using DirichletConstraintsAssember = Dune::PDELab::ConformingDirichletConstraints;
  using P1_dirichlet_GFS = Dune::PDELab::GridFunctionSpace<GV, P1_FEM, DirichletConstraintsAssember, VectorBackendP1>;
  P1_dirichlet_GFS p1_dirichlet_gfs_(gv, p1_fem);
  p1_dirichlet_gfs_.name("p1_dirichlet_gfs_");
  
  // Set up constraints container...  
  using P1_dirichlet_GFS_CC = P1_dirichlet_GFS::ConstraintsContainer<RangeType>::Type;
  P1_dirichlet_GFS_CC p1_dirichlet_gfs__cc;
  p1_dirichlet_gfs__cc.clear();
  auto p1_bctype_lambda = [&](const auto& x){ return 1.0; };
  auto p1_bctype = Dune::PDELab::makeBoundaryConditionFromCallable(gv, p1_bctype_lambda);
  Dune::PDELab::constraints(p1_bctype, p1_dirichlet_gfs_, p1_dirichlet_gfs__cc);
  
  // Set up grid grid operators...  
  using LOP_R = rOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, RangeType>;
  using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
  using GO_r = Dune::PDELab::GridOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, LOP_R, MatrixBackend, DF, RangeType, RangeType, P1_dirichlet_GFS_CC, P1_dirichlet_GFS_CC>;
  LOP_R lop_r(p1_dirichlet_gfs_, p1_dirichlet_gfs_, initree);
  p1_dirichlet_gfs_.update();
  int generic_dof_estimate =  6 * p1_dirichlet_gfs_.maxLocalSize();
  int dofestimate = initree.get<int>("istl.number_of_nnz", generic_dof_estimate);
  MatrixBackend mb(dofestimate);
  GO_r go_r(p1_dirichlet_gfs_, p1_dirichlet_gfs__cc, p1_dirichlet_gfs_, p1_dirichlet_gfs__cc, lop_r, mb);
  std::cout << "gfs with " << p1_dirichlet_gfs_.size() << " dofs generated  "<< std::endl;
  std::cout << "cc with " << p1_dirichlet_gfs__cc.size() << " dofs generated  "<< std::endl;
  
  // Set up solution vectors...  
  using V_R = Dune::PDELab::Backend::Vector<P1_dirichlet_GFS,DF>;
  V_R x_r(p1_dirichlet_gfs_);
  x_r = 0.0;
  auto lambda_0000 = [&](const auto& x){ return (double)exp((-1.0) * ((0.5 - x[1]) * (0.5 - x[1]) + (0.5 - x[0]) * (0.5 - x[0]))); };
  auto func_0000 = Dune::PDELab::makeGridFunctionFromCallable(gv, lambda_0000);
  Dune::PDELab::interpolate(func_0000, p1_dirichlet_gfs_, x_r);
  
  // Set up (non)linear solvers...  
  using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_SuperLU;
  using SLP = Dune::PDELab::StationaryLinearProblemSolver<GO_r, LinearSolver, V_R>;
  LinearSolver ls(false);
  double reduction = initree.get<double>("reduction", 1e-12);
  SLP slp(go_r, ls, x_r, reduction);
  slp.apply();
  
  // Do visualization...  
  using VTKWriter = Dune::SubsamplingVTKWriter<GV>;
  Dune::RefinementIntervals subint(initree.get<int>("vtk.subsamplinglevel", 1));
  VTKWriter vtkwriter(gv, subint);
  std::string vtkfile = initree.get<std::string>("wrapper.vtkcompare.name", "output");
  CuttingPredicate predicate;
  Dune::PDELab::addSolutionToVTKWriter(vtkwriter, p1_dirichlet_gfs_, x_r, Dune::PDELab::vtk::defaultNameScheme(), predicate);
  vtkwriter.write(vtkfile, Dune::VTK::ascii);
  
  // Maybe print residuals and matrices to stdout...  
  if (initree.get<bool>("printresidual", false)) {
    using Dune::PDELab::Backend::native;
    V_R r(x_r);
    // Setup random input
    std::size_t seed = 0;
    auto rng = std::mt19937_64(seed);
    auto dist = std::uniform_real_distribution<>(-1., 1.);
    for (auto& v : x_r)
      v = dist(rng);
    r=0.0;
    go_r.residual(x_r, r);
    Dune::printvector(std::cout, native(r), "residual vector", "row");
  }
  if (initree.get<bool>("printmatrix", false)) {
    // Setup random input
    std::size_t seed = 0;
    auto rng = std::mt19937_64(seed);
    auto dist = std::uniform_real_distribution<>(1., 10.);
    for (auto& v : x_r)
      v = dist(rng);
    using M = typename GO_r::Traits::Jacobian;
    M m(go_r);
    go_r.jacobian(x_r,m);
    using Dune::PDELab::Backend::native;
    Dune::printmatrix(std::cout, native(m),"global stiffness matrix","row",9,1);
  }
  
  // Maybe calculate errors for test results...  
  using P1_DIRICHLET_GFS__DGF = Dune::PDELab::DiscreteGridFunction<decltype(p1_dirichlet_gfs_),decltype(x_r)>;
  P1_DIRICHLET_GFS__DGF p1_dirichlet_gfs__dgf(p1_dirichlet_gfs_,x_r);
  using DifferenceSquaredAdapter_ = Dune::PDELab::DifferenceSquaredAdapter<decltype(func_0000), decltype(p1_dirichlet_gfs__dgf)>;
  DifferenceSquaredAdapter_ dsa_(func_0000, p1_dirichlet_gfs__dgf);
  RangeType l2error(0.0);
  {
    // L2 error squared of difference between numerical
    // solution and the interpolation of exact solution
    // for treepath ()
    typename decltype(dsa_)::Traits::RangeType err(0.0);
    Dune::PDELab::integrateGridFunction(dsa_, err, 10);
  
    l2error += err;
    if (gv.comm().rank() == 0){
      std::cout << "L2 Error for treepath : " << err << std::endl;
    }}
  bool testfail(false);
  using std::abs;
  using std::isnan;
  if (gv.comm().rank() == 0){
    std::cout << "\nl2errorsquared: " << l2error << std::endl << std::endl;
  }
  if (isnan(l2error) or abs(l2error)>1e-7)
    testfail = true;
  return testfail;
  
}

This fixes #109 (closed)