Integrating a 2D function along a single dimension

Question

I was wondering if there is a better way of computing integrals of functions in one dimension only. Basically, I have a function $u(x,y)$ defined on an (irregular) triangulation of the unit square, and I wish to compute the integral
$$
U(x) := \int_0^1 u(x,y) dy.
$$
My current solution is to evaluate u along a uniform grid of points, and do trapezoidal quadrature along the y-axis. I'd really like to know if there was a neater/more efficient solution!

Answer 1

Hi,

My solution and an exemple with u(x,y)=x+y

from dolfin import *
import numpy

#2D domain to define u(x,y)
m_2D=RectangleMesh(0.,0.,1.,1.,10,10)

V_2D=FunctionSpace(m_2D, "CG",1)

U_2D_e=Expression('x[0]+x[1]')

U_2D_f=interpolate(U_2D_e,V_2D)

#1D integration domain
m_1D=IntervalMesh(20,0.,1.)
V_1D=FunctionSpace(m_1D, "CG",1)

#an Expression class to evaluate u(x,y) on the 1D domain (range y, x fixed)

class my1DExpression(Expression):
    def __init__(self,u_2d,x):
        self.u_2d=u_2d
        self.x=x
        Expression.__init__(self)
        self._vx= numpy.array([0.])
        self._pt_x= numpy.array([0.,0.])
    def eval(self, values, x):
        self._pt_x[0]= self.x
        self._pt_x[1]= x[0]
        self.u_2d.eval(self._vx,self._pt_x)
        values[0] = self._vx[0]



#An exemple int_0^1 x+y dy=x+1/2
for x in numpy.linspace(0.,1.,10):
    Uy=interpolate(my1DExpression(U_2D_f,x=x),V_1D)
    intUy=assemble(Uy*dx)
    print intUy,x+0.5,intUy-(x+0.5)

#it is possible to define the u(x) function
class UxExpression(Expression):
    def __init__(self,u_2d):
        self.u_2d=u_2d
        self.x=x
        Expression.__init__(self)
    def eval(self, values, x):
        _Uy=interpolate(my1DExpression(U_2D_f,x=x[0]),V_1D)
        values[0] = assemble(_Uy*dx)

#1D for ux
m_1Dx=IntervalMesh(20,0.,1.)
V_1Dx=FunctionSpace(m_1Dx, "CG",1)

Ux_f=interpolate(UxExpression(U_2D_f),V_1Dx)
plot(Ux_f,interactive=True)