Integrate along axes/lines

Question

Dear Fenics-Community,

I have a (scalar) solution of my FEM which is defined on a 3D mesh. How can I obtain integrals (averages) along a spatial axis (or any line)?

I have tried to extend the plane average approach presented here by using the method described there to obtain a slice of my domain and place my line of interest on this slice, but have not be successful yet. Any ideas would be much appreciated - thanks in advance.

Best wishes,
Merlin :)

Answer 1

Hi, one possibility is as follows

from dolfin import *
import numpy as np

def line_integral(u, A, B, n): 
    '''Integrate u over segment [A, B] partitioned into n elements'''
    assert u.value_rank() == 0
    assert len(A) == len(B) > 1 and np.linalg.norm(A-B) > 0
    assert n > 0

    # Mesh line for integration
    mesh_points = [A + t*(B-A) for t in np.linspace(0, 1, n+1)]
    tdim, gdim = 1, len(A)
    mesh = Mesh()
    editor = MeshEditor()
    editor.open(mesh, tdim, gdim)
    editor.init_vertices(n+1)
    editor.init_cells(n)

    for vi, v in enumerate(mesh_points): editor.add_vertex(vi, v)

    for ci in range(n): editor.add_cell(ci, np.array([ci, ci+1], dtype='uintp'))

    editor.close()

    # Setup funcion space
    elm = u.function_space().ufl_element()
    family = elm.family()
    degree = elm.degree()
    V = FunctionSpace(mesh, family, degree)
    v = interpolate(u, V)

    return assemble(v*dx)

# ----------------------------------------------------------------------------

if __name__ == '__main__':
    mesh = UnitCubeMesh(2, 4, 5)
    V = FunctionSpace(mesh, 'CG', 1)
    f = Expression('x[0]+x[1]+x[2]', degree=1)
    v = interpolate(f, V)

    A = np.array([0, 0, 0])
    B = np.array([1, 1, 1])

    ans = line_integral(v, A, B, n=100)
    ans0 = 3*sqrt(3)/2
    print 'Error', abs(ans-ans0)