 
 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);

     *
     * @return cusp matrix
     */
    dg::Elliptic<dg::CartesianGrid2d, Matrix, container>& laplacianM( ) { return laplaceM;}

    /**
     * @brief Return the Gamma operator used by this object
     *
     * @return Gamma operator
     */
    dg::Helmholtz<dg::CartesianGrid2d, Matrix, container >&  gamma() {return gamma1;}

    /**
     * @brief Compute the right-hand side of the toefl equations

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


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