Integration with bivariate integrand

I/O
1

Hi all,

I’m solving cardiac electrophysiology problems using cbcbeat. I want to post-process a potential, for which I need to evaluate the following integral:

f(\vec x) = \int_\varOmega \frac{I(\vec \xi)}{\| \vec x - \vec \xi\|} ~ \text{d}\vec \xi

where I is a dolfin function calculated from my solution. The function f is basically an inverse distance weighting of I.

I can evaluate f at a single point \vec x using the assemble function in dolfin. However, I want to evaluate f on all points on my mesh. Do you know how I can do this while avoiding a costly for-loop?

Thank you very much for your help!

2

You might be able to re-formulate it as the solution to a PDE, then solve the PDE approximately using FEniCS. In 3D, f would be the solution to the Poisson problem with a source term of I (up to a constant factor involving \pi), but on an infinite domain. Look up the term “fundamental solution”.

You can get surprisingly-okay solutions for infinite domain problems by introducing geometrical mappings, as in this code example I wrote for a question on the old forum:

from dolfin import *
mesh = UnitSquareMesh(50,50)

# (The mapping tends to pump up the estimated quadrature degree to an
# impractical level.)
dx = dx(metadata={"quadrature_degree":4})

# Convergence appears to hold for messy unstructured meshes too
#from mshr import *
#mesh = generate_mesh(Rectangle(Point(0,0),Point(1,1)),128)

# Use to define an exact solution decaying to zero at infinity
def gaussian(x,sigma):
    r = x[0]**2 + x[1]**2
    return exp(-0.5*((r/sigma)**2))/(sigma*sqrt(2.0*pi))

# Map the unit interval to the entire real line
def I_to_R(x):
    return tan(pi*(x-0.5))

# Map the unit square to R^2
def I2_to_R2(xhat):
    return as_vector([I_to_R(xhat[0]),I_to_R(xhat[1])])

# Map mesh coordinates to infinite domain
xhat = SpatialCoordinate(mesh)
x = I2_to_R2(xhat)
Dx = grad(x)

# Define grad and div w.r.t. x
def gradx(f):
    return dot(grad(f),inv(Dx))
def divx(f):
    n = rank(f)
    ii = indices(n)
    return as_tensor(gradx(f)[ii+(ii[n-1],)],ii[0:n-1])

# Integral change-of-variables
J = det(Dx)

# Manufacture source term for the Gaussian exact solution
u_exact = gaussian(x,0.2)
f = -divx(gradx(u_exact))

# Define weak Poisson problem
V = FunctionSpace(mesh,"Lagrange",2)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(gradx(u),gradx(v))*J*dx
L = inner(f,v)*J*dx

# Homogeneous Dirichlet BCs on parametric domain (so, at infinity in
# physical domain); this only really makes sense for a 2D electric potential
# if the net charge is zero, since the free-space Green's function does
# not decay to zero until at least three dimensions.  
bc = DirichletBC(V,0.0,"on_boundary")

# Solve
u = Function(V)
solve(a==L,u,bcs=[bc,])

# Check L^2 error on bi-unit square in physical domain; some quick
# mesh-doubling tests indicate that this converges optimally.
import math
chi = conditional(gt(Max(abs(x[0]),abs(x[1])),1.0),Constant(0.0),Constant(1.0))

# Error in full domain also appears to converge despite J blowing up at edges
#chi = Constant(1.0)

print("L^2 error: "+str(math.sqrt(assemble(chi*((u-u_exact)**2)*J*dx))))

print("Net charge: "+str(assemble(f*J*dx)))

####### Visualization #######

# Output for visualizing the solution
u.rename("u","u")
File("u.pvd") << u

# Using the standard project() function is oscillatory for singular functions
def lumpedProject(f):
    V = FunctionSpace(mesh,"Lagrange",1)
    v = TestFunction(V)
    lhs = assemble(inner(Constant(1.0),v)*dx)
    rhs = assemble(inner(f,v)*dx)
    u = Function(V)
    as_backend_type(u.vector())\
        .vec().pointwiseDivide(as_backend_type(rhs).vec(),
                               as_backend_type(lhs).vec())
    return u

# Regularize mapping used for visualization slightly, to avoid too much
# craziness near the edges:
cx = Constant((0.5,0.5))
x_reg = I2_to_R2(0.9999*(xhat-cx)+cx)
dispx = lumpedProject(x_reg[0]-xhat[0])
dispy = lumpedProject(x_reg[1]-xhat[1])
dispx.rename("dispx","dispx")
dispy.rename("dispy","dispy")
File("dispx.pvd") << dispx
File("dispy.pvd") << dispy

# To plot, load all three files, apply an Append Attributes filter to them, use
# a Calculator filter to define the vector field
#
#   dispx*iHat + dispy*jHat
#
# and use this field in a Warp by Vector filter.  Set the warping scale to 1,
# color by u, and zoom in to the domain of interest.  Using the
# "Surface with Edges" mode to visualize shows the effect of the
# geometrical mapping clearly.

On the other hand, if the formula comes out of some textbook, you might also consider whether the assumption of an infinite domain was introduced there “for simplicity”, and whether applying specific BCs on a finite domain might make more sense in your application.

EDIT: Also, evaluating the integral for fixed x using assemble() is probably not very accurate, because the FE quadrature is not designed for singular integrands.

3

Thank you for this smart solution!