Subdomain Integrals in Fenics

Question

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)

    a = inner(grad(du), v)*dx
    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 01:24

score 1 · Answer 1 · answered Jun 19 '13 at 20:54

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

define the 1D FEM basis $\{\phi_i(y)\}$
for every $\{\phi_i(y)\}$ define a 2D expression $b_i(x,y) = \phi_i(y)$
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 $$
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

score 0 · Answer 2 · edited Apr 13 '17 at 12:53

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

Subdomain Integrals in Fenics

2 Answers2