             dg::blas1::axpby( -nu_, y[i], 0., y[i]);
         }
     }
-    dg::Elliptic<dg::CartesianGrid2d, Matrix, container>& laplacianM() {return LaplacianM_perp;}
-    const container& weights(){return w2d;}
-    const container& precond(){return v2d;}
+    dg::Elliptic<Geometry, Matrix, container>& laplacianM() {return LaplacianM_perp;}
+    const container& weights(){return LaplacianM_perp.weights();}
+    const container& precond(){return LaplacianM_perp.precond();}
 
   private:
     double nu_;
     bool global;
     const container w2d, v2d;
-    container temp, expx;
-    dg::Elliptic<dg::CartesianGrid2d, Matrix, container> LaplacianM_perp;
+    container temp;
+    dg::Elliptic<Geometry, Matrix, container> LaplacianM_perp;
 };
 
 template< class Matrix, class container >

     */
    ToeflR( const CartesianGrid2d& g, double kappa, double nu, double tau, double eps_pol, double eps_gamma, int global);

    /**
     * @brief Exponentiate pointwise every Vector in src 
     *
     * @param src source
     * @param dst destination may equal source
     */
    void exp( const std::vector<container>& src, std::vector<container>& dst);

    /**
     * @brief Take the natural logarithm pointwise of every Vector in src 
     *
     * @param src source
     * @param dst destination may equal source
     */
    void log( const std::vector<container>& src, std::vector<container>& dst);

    void divide( const container& zaehler, const container& nenner, container& result);

    /**
     * @brief Returns phi and psi that belong to the last y in operator()


namespace dg
{
template<class Geometry, class Matrix, class container>
struct Diffusion
{
    Diffusion( const Geometry& g, double nu, bool global):
        nu_(nu), global(global), 
        temp( dg::evaluate(dg::zero, g)), 
        LaplacianM_perp( g, dg::normed, dg::centered){
    }
    void operator()( std::vector<container>& x, std::vector<container>& y)
    {
        /* x[0] := N_e - 1
         * x[2] := N_i -1 
         */
        for( unsigned i=0; i<x.size(); i++)
        {
            dg::blas2::gemv( LaplacianM_perp, x[i], temp);
