3
$\begingroup$

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()
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Peter Brune Jun 19 '13 at 20:55
  • $\begingroup$ 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 $\endgroup$ – David M. Rogers Jun 19 '13 at 23:56
  • $\begingroup$ 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... $\endgroup$ – David M. Rogers Jun 20 '13 at 1:24
1
$\begingroup$

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)
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – David M. Rogers Jun 19 '13 at 22:26
  • $\begingroup$ 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. $\endgroup$ – Jan Jun 20 '13 at 12:29
0
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.