Correctly pass grid function types into the local operator

python/dune/perftool/pdelab/driver/solve.py

@@ -9,6 +9,7 @@ from dune.perftool.pdelab.driver import (get_form_ident,
                                          name_initree,
                                          )
 from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
+                                                           type_domainfield,
                                                            type_trial_gfs,
                                                            )
 from dune.perftool.pdelab.driver.constraints import (type_constraintscontainer,
@@ -82,8 +83,9 @@ def name_vector(form_ident):
 
 @preamble
 def typedef_vector(name, form_ident):
-    go_type = type_gridoperator(form_ident)
-    return "using {} = {}::Traits::Domain;".format(name, go_type)
+    gfs = type_trial_gfs()
+    df = type_domainfield()
+    return "using {} = Dune::PDELab::Backend::Vector<{},{}>;".format(name, gfs, df)
 
 
 def type_vector(form_ident):

python/dune/perftool/pdelab/localoperator.py

@@ -157,6 +157,26 @@ def localoperator_basename(form_ident):
     return get_form_option("classname", form_ident)
 
 
+def name_gridfunction_member(coeff):
+    name = "gridfunction_{}".format(coeff.count())
+    define_gridfunction_member(name)
+    return name
+
+
+@class_member(classtag="operator")
+def define_gridfunction_member(name):
+    _type = type_gridfunction_template_parameter(name)
+    param = "{}_".format(name)
+    constructor_parameter("const {}&".format(_type), param, classtag="operator")
+    initializer_list(name, [param], classtag="operator")
+    return "const {}& {};".format(_type, name)
+
+
+@template_parameter(classtag="operator")
+def type_gridfunction_template_parameter(name):
+    return name.upper()
+
+
 def class_type_from_cache(classtag):
     from dune.perftool.generation import retrieve_cache_items
 
@@ -702,10 +722,14 @@ def generate_localoperator_kernels(operator):
     lop_template_test_gfs()
     lop_template_range_field()
 
-    # Make sure there is always the same constructor arguments (even if parameter class is empty)
-    from dune.perftool.pdelab.localoperator import name_initree_member
+    # Make sure there is always the same constructor arguments, even if some of them are
+    # not strictly needed. Also ensure the order.
     name_initree_member()
 
+    # Iterate over the needed grid functions in correct order
+    for c in sorted(filter(lambda c: c.count() > 2, form.coefficients()), key=lambda c: c.count()):
+        name_gridfunction_member(c)
+
     # Set some options!
     from dune.perftool.pdelab.driver import isQuadrilateral
     if isQuadrilateral(form.coefficients()[0].ufl_element().cell()):

test/coeffeval/coeffeval_poisson.cc

@@ -18,64 +18,95 @@
 #include "poisson_localoperator.hh"
 #include "dune/pdelab/stationary/linearproblem.hh"
 #include "dune/grid/io/file/vtk/subsamplingvtkwriter.hh"
+#include "dune/pdelab/function/discretegridviewfunction.hh"
 
 
 int main(int argc, char** argv)
 {
+  // MPI helper stuff
   Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
-  using VectorBackendP1 = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::none>;
+
+  // Parse the ini file
+  Dune::ParameterTree initree;
+  Dune::ParameterTreeParser::readINITree(argv[1], initree);
+
+  // Build a grid
   using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
   using GV = Grid::LeafGridView;
+  IniGridFactory<Grid> factory(initree);
+  std::shared_ptr<Grid> grid = factory.getGrid();
+  GV gv = grid->leafGridView();
+
+  // General types and stuff
   using DF = Grid::ctype;
   using RangeType = double;
+
+  // Finite Element Maps
   using P1_FEM = Dune::PDELab::PkLocalFiniteElementMap<GV, DF, RangeType, 1>;
+  P1_FEM p1_fem(gv);
+
+  // 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>;
-  using LOP_R = PoissonLocalOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, RangeType>;
-  using P1_dirichlet_GFS_CC = P1_dirichlet_GFS::ConstraintsContainer<RangeType>::Type;
-  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>;
-  using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_SuperLU;
-  using V_R = GO_r::Traits::Domain;
-  using SLP = Dune::PDELab::StationaryLinearProblemSolver<GO_r, LinearSolver, V_R>;
-  Dune::ParameterTree initree;
-  Dune::ParameterTreeParser::readINITree(argv[1], initree);
-  IniGridFactory<Grid> factory(initree);
-  std::shared_ptr<Grid> grid = factory.getGrid();
-  GV gv = grid->leafGridView();
-  P1_FEM p1_fem(gv);
   P1_dirichlet_GFS p1_dirichlet_gfs_(gv, p1_fem);
   p1_dirichlet_gfs_.name("p1_dirichlet_gfs_");
+  p1_dirichlet_gfs_.update();
+  std::cout << "gfs with " << p1_dirichlet_gfs_.size() << " dofs generated  "<< std::endl;
+
+  // Solution vectors / Grid Functions
+  using V_R = Dune::PDELab::Backend::Vector<P1_dirichlet_GFS,DF>;
+  V_R x_r(p1_dirichlet_gfs_);
+  V_R sol(p1_dirichlet_gfs_);
+
+  using GF = Dune::PDELab::DiscreteGridViewFunction<P1_dirichlet_GFS, V_R>;
+  GF sol_gf(p1_dirichlet_gfs_, sol);
+
+  // Local Operator
+  using LOP_R = PoissonLocalOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, RangeType, GF>;
+  LOP_R lop_r(p1_dirichlet_gfs_, p1_dirichlet_gfs_, initree, sol_gf);
+
+  // Constraints stuff
+  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);
-  LOP_R lop_r(p1_dirichlet_gfs_, p1_dirichlet_gfs_, initree);
-  p1_dirichlet_gfs_.update();
+  std::cout << "cc with " << p1_dirichlet_gfs__cc.size() << " dofs generated  "<< std::endl;
+
+  // Matrix Backend
+  using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
   int generic_dof_estimate =  6 * p1_dirichlet_gfs_.maxLocalSize();
   int dofestimate = initree.get<int>("istl.number_of_nnz", generic_dof_estimate);
   MatrixBackend mb(dofestimate);
+
+  // Grid Operator
+  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>;
   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;
+
+  // Solver
+  using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_SuperLU;
   LinearSolver ls(false);
-  V_R x_r(p1_dirichlet_gfs_);
-  V_R sol(p1_dirichlet_gfs_);
+  using SLP = Dune::PDELab::StationaryLinearProblemSolver<GO_r, LinearSolver, V_R>;
+
+  // Interpolation
   x_r = 0.0;
   std::size_t seed = 0;
   auto rng = std::mt19937_64(seed);
   auto dist = std::uniform_real_distribution<>(-1., 1.);
   for (auto& v : sol)
     v = dist(rng);
-
   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_, sol);
 
+  // Solving
   double reduction = initree.get<double>("reduction", 1e-12);
   SLP slp(go_r, ls, x_r, reduction);
   slp.apply();
+
+  // VTK visualization
   using VTKWriter = Dune::SubsamplingVTKWriter<GV>;
   Dune::RefinementIntervals subint(initree.get<int>("vtk.subsamplinglevel", 1));
   VTKWriter vtkwriter(gv, subint);
@@ -83,6 +114,8 @@ int main(int argc, char** argv)
   CuttingPredicate predicate;
   Dune::PDELab::addSolutionToVTKWriter(vtkwriter, p1_dirichlet_gfs_, x_r, Dune::PDELab::vtk::defaultNameScheme(), predicate);
   vtkwriter.write(vtkfile, Dune::VTK::ascii);
+
+  // Error calculation
   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);
   RangeType l2error(0.0);