 
 namespace dg
 {
-template<class Matrix, class container>
+template<class Geometry, class Matrix, class container>
 struct Diffusion
 {
-    Diffusion( const dg::CartesianGrid2d& g, double nu, bool global):
+    Diffusion( const Geometry& g, double nu, bool global):
         nu_(nu), global(global), 
-        w2d( dg::create::weights(g)), v2d( dg::create::inv_weights(g)), 
-        temp( dg::evaluate(dg::zero, g)), expx(temp),
+        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);

     */
    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()

#endif 
    if( global)
    {
        blas1::axpby( -1., y[0], 1., gamma_n, omega); //omega = a_i\Gamma n_i - n_e
        if( global == 2)
            divide( omega, expy[1], omega);
        blas2::symv( w2d, omega, omega);
