 #include <cusp/transpose.h>
 
 #include "functions.h"
+#include "tensor.cuh"
 #include "operator.cuh"
 #include "operator_matrix.cuh"
 #include "creation.cuh"

    cusp::multiply( middle, righty, righty);
    cusp::transpose( rightx, leftx); 
    cusp::transpose( righty, lefty); 
    Operator<T,n> weights(0.);
    for( unsigned i=0; i<n; i++)
        weights(i,i) = DLT<n>::weight[i]*h/2.; // normalisation because F is invariant
    Operator<T, n*n> weights2d = tensor( weights, weights);
    middle = tensor<T,n>( Nx*Ny, weights2d*backward2d);
    jumpx = create::jump_ot<T,n>( Nx, bcx); //without jump cg is unstable
    jumpx = dg::dgtensor( tensor<T,n>(Ny, delta), jumpx);
    jumpy = create::jump_ot<T,n>( Ny, bcy); //without jump cg is unstable
    jumpx = dg::dgtensor( jumpy, tensor<T,n>(Nx, delta));

}


{
    typedef cusp::coo_matrix<int, T, Memory> Matrix;
    typedef cusp::array1d<T, Memory> Vector;
    Polarisation2d( unsigned Nx, unsigned Ny, T hx, T hy, int bcx, int bcy);
    Matrix create( const Vector& );
  private:
    typedef cusp::array1d<int, Memory> Array;
