# Subdomain Integrals in Fenics

I'd like to compute $u(y) = \int f(x,y) dx$, but don't see a really efficient way to do this in FEniCS. My first thought was to just assemble(f*v*dx) for two different meshes. (try1, below), but it seems I can't combine 1D and 2D meshes, and that the DG meshes have to have at least 2 basis functions in each direction, and I'm not sure how the mesh topology -> basis coefficient layout is controlled.

My next thought was to solve du/dx = f with a zero-BC on one side, to accumulate the integral on x = 1. That one gets all the way through assembling the problem, then fails with a singular matrix exception.

Is there a better way to do this in FEniCS? To solve this problem in general, I would suggest adding general convolutions to FEniCS, so that $u = \int h(\vec x-\vec y) v(\vec y) dy$ could be done by an appropriate matrix-multiply...

from dolfin import *

m = UnitSquareMesh(12, 12)
V = FunctionSpace(m, "Lagrange", 1)

f = Function(V)
f.interpolate(Expression("x[1]-x[0]")) # int_0^1 f dx = x[1] - 0.5

def try1():
m2 = UnitSquareMesh(1, 12)
V2 = FunctionSpace(m2, "DG", 0)
u = Function(V2)
print u.vector().array().shape
v = TestFunction(V2)
A = assemble(f*v*dx)
print A.array()*12.0

def try2():
du = TrialFunction(V)
v = TestFunction(VectorFunctionSpace(m, "Lagrange", 1))

bc = DirichletBC(V, Constant(0.0), lambda x,b: b and x[0] < 1e-8)
#bc2 = DirichletBC(V, Constant(0.0), lambda x,b: b and x[1] < 1e-8)

L = f*v[0]*dx
#L2 = f*v[1]*dx

u = Function(V)

solve(a == L, u, bcs=bc)

print u.vector().array()

• Could you be more precise with respect to what you're trying to solve? I don't think that your code is doing what you think it is. – Peter Brune Jun 19 '13 at 20:55
• Adding these lines to the end of try1 actually solves the problem. But it only works if the unitsquare mesh is laid out in a (12 x 2) shape. print np.sum(np.reshape(A.array(), (-1,2)), 1)*12.0 print (np.arange(12)+0.5)/12.0 - 0.5 – David M. Rogers Jun 19 '13 at 23:56
• I see, so "np.sum(np.reshape(A.array(), (-1,2)), 1)*12.0" only accidentally solves the problem, since "f" is is linear, and only evaluated at maybe the centroid of triangles in m2 during assemble... – David M. Rogers Jun 20 '13 at 1:24

I don't think there is a generic way to do this, since the (discrete) representation of $u(y)$ has to be defined beforehand. If you want it expressed in a 1D finite element basis then you can do the following

1. define the 1D FEM basis $\{\phi_i(y)\}$
2. for every $\{\phi_i(y)\}$ define a 2D expression $b_i(x,y) = \phi_i(y)$
3. test every $b_i$ against your $f$ to get $$\int \int b_i(x,y)f(x,y)dxdy=\int\phi_i(y)u(y)dy =: \tilde u_i$$
4. compute the coefficients of the expansion of $u$ in the FEM basis - $u(y) \approx \sum_i u_i \phi_i(y)$ - via $$\begin{bmatrix} \vdots \\ u_i \\ \vdots \end{bmatrix} = \mathcal M_{\phi}^{-1} \begin{bmatrix} \vdots \\ \tilde u_i \\ \vdots \end{bmatrix},$$ where $M_\phi$ is the FEM mass matrix.

If you have the expression b_i in fenics and function f, then you find

tu_i = assemble(b_i*f*dx)

• Of course, I could just compute the inner product of f with a bunch of basis functions, but this becomes impractical because I'd be casting a 1D function to 2 (or 3 is where I'm going) dimensions and repeating the function call N times. Python is supposed to be high-level and vectorized. – David M. Rogers Jun 19 '13 at 22:26
• The problem seems to broad to provide a general low-level functionality for this. What if you are not on a rectangular domain aligned with the coordinate axes? However, I might be wrong and the fenics Q&A forum has the solution for your needs. – Jan Jun 20 '13 at 12:29

Another way would be to solve the Poisson equation using $f$ as a source, then use Gauss' law to compute the flux across a line segment for each dividing plane. That method still wouldn't be vectorized, but at least only requires 1 mesh function (for marking domains). If I make that compromise, I may as well just iterate in python:

def slice(mesh, n, L):
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
subdomains.set_all(0)
dx = L/float(n)
for i in range(1, n):
b = AutoSubDomain(lambda x: x[1] > i*dx-1e-8)
b.mark(subdomains, i)
return subdomains

def try3():
markers = slice(m, 12, 1.0)
dX = dx[markers]
v = [assemble(f*dX(i))*12.0 for i in range(12)]
print v