[dune-pdelab] Finite Volume Method with slope reconstruction

[Gregor Corbin]

Hello dear Dune developers,

I have a dune-pdelab code that implements a first-order finite volume 
method for a hyperbolic system. Currently I am trying to implement a 
slope reconstruction to achieve second order accuracy in space. 
Therefore, before each stage of the time-stepping scheme, I need to 
iterate over the grid and compute the slopes on each cell.
My current approach is to write a limiter class and give it to the 
explicit one step method like this:

oneStepMethod.apply(t,dt,x_old,x_new,LinearReconstruction);

In the prestage method of this class, which gets only the solution 
vector, i iterate over all the cells on the grid. For each cell I want 
to extract the solution on that cell and on it's neighbors.
But I have no idea how to do this.

I appreciate any help on this problem. Is there maybe even a better way 
to implement slope limiting in pdelab?

Cheers and thanks in advance,
Gregor Corbin

Below is a sketch for the code of the limiter class:

template <typename GridView, typename GFS, int n_components>
class LinearReconstruction
{
   enum {dim = GridView::dimension};
   enum {n_comp = n_components};

public:
   LinearReconstruction(const GridView &grid_view, const GFS &gfs):
   grid_view_(grid_view),
   gfs_(gfs),
   mapper_(grid_view_),
   slopes_(mapper_.size(),0.0)
   {}

   template<typename X>
   void prestage(X& x)
   {
     // reconstruct slopes
     for (const auto& cell : elements(grid_view_) )
     {
        // I would like to do something like in the local operator methods,
        // but how do I get the lfsu_s ???
        // is there a method to get it from the cell and the grid 
function space?

        // lfsu_s = gfs_.localFunctionSpace(cell); <- does something 
similar exist?
       Dune::FieldVector<double,n_comp> u_s(0.0);
       for (int k = 0; k < n_comp; ++k)
         u_s[k] = x(lfsu_s.child(k),0);

       // get the solutions from all the neighbors
       std::vector<Dune::FieldVector<double,n_comp> > 
u_n(pow(2,dim),0.0); // only valid for structured cubic grids
       int i = 0;
       for (const auto& is : intersections(grid_view_,cell))
       {
         const auto &neighbor = is.outside();
          // lfsu_n = gfs_.localFunctionSpace(neighbor); <- does 
something similar exist?

          for (int k = 0; k < n_comp; ++k)
               u_n.at(i)[k] = x(lfsu_n.child(k),0);
          ++i;
       }

       // use the info to reconstruct slopes
       slopes_.at(mapper_.index(cell)) =  // some function(u_s, u_n)
     }
   }

   template<typename V>
   void poststage(V& v)
   {}
private:
   const GridView &grid_view_;
    const GFS &gfs_;
   using Mapper =
Dune::MultipleCodimMultipleGeomTypeMapper<GridView,Dune::MCMGElementLayout>;
   const Mapper mapper_;
   std::vector<Dune::FieldMatrix<double,n_comp,dim> > slopes_;
};