         detail::legendre( px, xn[row], n, forward);
         
         int offset = row;
-        unsigned col_begin = celln[row]*n;
         for(int k=0; k<n; k++)
         {
-            Aj[offset] = col_begin + k*n;
+            Aj[offset] = cellx[row]*n + k
             Av[offset] = px[k];
             offset +=pitch;
         }

    double * Av = thrust::raw_pointer_cast(&A.values(0,0));

    //compute normalized x and y values and cell numbers
    thrust::device_vector<double> xn(x.size()), yn(y.size()), cellnh(x.size()), cellmh(x.size());
    dg::blas1::transform( x, xn, dg::PLUS<double>( -g.x0()));
    dg::blas1::transform( y, yn, dg::PLUS<double>( -g.y0()));
    dg::blas1::scal( xn, 1./g.hx());

        detail::legendre( py, yn[row], n, forward);
        
        int offset = row;
        int col_begin = celly[row]*Nx*n*n + cellx[row]*n;
        //evaluate 2d Legendre polynomials at (xn)...
        for( int k=0; k<n; k++)
            for( int l=0; l<n; l++)
