template <class VectorDiscreteFunctionSpace, class EntityType, class QuadratureType >
MatrixType massMatrix( const EntityType& e, VectorDiscreteFunctionSpace& space, QuadratureType& quadrature )
{
    typedef typename VectorDiscreteFunctionSpace::RangeType
        RangeType;
    typedef typename VectorDiscreteFunctionSpace::DomainType
        ElementCoordinateType;
    const typename VectorDiscreteFunctionSpace::BaseFunctionSetType& baseset = space.baseFunctionSet( e );
    int numDofs = baseset.numBaseFunctions();
    MatrixType a( 0. );
    RangeType psi;
    for (int qp =0; qp < quadrature.nop() ; ++qp ) {
        ElementCoordinateType x = quadrature.point(qp);
    	double weight = quadrature.weight(qp);
    	double intel = e.geometry().integrationElement( x );

    	for ( int i = 0; i < numDofs; ++i) {
    	    baseset.evaluate(i,x,psi);
    	    for ( int j = 0; j < numDofs; ++j) {
                double val = baseset.evaluateSingle(j,x,psi);
                val *= intel * weight;
                a[i][j] += val;
            }
    	}
    }
    return a;
}

template <class VectorDiscreteFunctionSpace, class ScalarDiscreteFunctionSpace, class EntityType, class QuadratureType >
MatrixType gradMatrix( const EntityType& e, VectorDiscreteFunctionSpace& vector_space, ScalarDiscreteFunctionSpace& scalar_space, QuadratureType& quadrature )
{
    typedef typename ScalarDiscreteFunctionSpace::RangeType
        RangeType;
    typedef typename VectorDiscreteFunctionSpace::JacobianRangeType
        JacobianRangeType;
    typedef typename VectorDiscreteFunctionSpace::DomainType
        ElementCoordinateType;
    const typename VectorDiscreteFunctionSpace::BaseFunctionSetType& vector_baseset = vector_space.baseFunctionSet( e );
    const typename ScalarDiscreteFunctionSpace::BaseFunctionSetType& scalar_baseset = scalar_space.baseFunctionSet( e );
    const int dim = ScalarDiscreteFunctionSpace::GridPartType::GridType::dimension;
    int vector_numDofs = vector_baseset.numBaseFunctions();
    int scalar_numDofs = scalar_baseset.numBaseFunctions();
    MatrixType a( 0. );
    JacobianRangeType helpmat(0.0);
    RangeType phi;
    for (int qp =0; qp < quadrature.nop() ; ++qp ) {
        ElementCoordinateType x = quadrature.point(qp);
    	double weight = quadrature.weight(qp);
    	double intel = e.geometry().integrationElement( x );
    	for ( int j = 0; j < scalar_numDofs; ++j) {
    	    scalar_baseset.evaluate(j,x,phi);

    	    for ( int k = 0; k < dim; ++k )
    	    	helpmat[k][k] = phi[0];

    	    for ( int i = 0; i < vector_numDofs; ++i) {
                double val = vector_baseset.evaluateGradientSingle(i,e,x,helpmat);
                val *= intel * weight;
                a[i][j] += val;
            }
    	}
    }
    return a;
}

template <class VectorDiscreteFunctionSpace, class ScalarDiscreteFunctionSpace, class EntityType, class QuadratureType >
MatrixType divMatrix( const EntityType& e, VectorDiscreteFunctionSpace& vector_space, ScalarDiscreteFunctionSpace& scalar_space, QuadratureType& quadrature )
{
    typedef typename VectorDiscreteFunctionSpace::RangeType
        RangeType;
    typedef typename VectorDiscreteFunctionSpace::JacobianRangeType
        JacobianRangeType;
    typedef typename VectorDiscreteFunctionSpace::DomainType
        ElementCoordinateType;
    const typename VectorDiscreteFunctionSpace::BaseFunctionSetType& vector_baseset = vector_space.baseFunctionSet( e );
    const typename ScalarDiscreteFunctionSpace::BaseFunctionSetType& scalar_baseset = scalar_space.baseFunctionSet( e );
    int vector_numDofs = vector_baseset.numBaseFunctions();
    int scalar_numDofs = scalar_baseset.numBaseFunctions();
    MatrixType a( 0. );
    RangeType phi;
    for (int qp =0; qp < quadrature.nop() ; ++qp ) {
        ElementCoordinateType x = quadrature.point(qp);
    	double weight = quadrature.weight(qp);
    	double intel = e.geometry().integrationElement( x );
    	for ( int j = 0; j < vector_numDofs; ++j) {
    	    vector_baseset.evaluate(j,x,phi);

    	    for ( int i = 0; i < scalar_numDofs; ++i) {
                double val = scalar_baseset.evaluateGradientSingle(i,e,x,phi[0]);
                val *= intel * weight;
                a[i][j] += val;
            }
    	}
    }
    return a;
}