Evaluating an integral like an expression

Question

Given a vector $\mathbf{F}$,

I want to evaluate the scalar $\Phi(\mathbf{r})$.

$\displaystyle \Phi\left(\mathbf{r}\right)=\frac{1}{4\pi}\int_{V}\frac{\boldsymbol{\nabla}'\cdot\mathbf{F}\left(\mathbf{r}'\right)}{\left|\mathbf{r}-\mathbf{r}'\right|}\mathrm{d}V'$

Over a unit square and interpolate/project it into u0 in the code below.

from dolfin import *

mesh = UnitSquareMesh(42, 42)
V = FunctionSpace(mesh, "CG", 1)
M = VectorFunctionSpace(mesh, "CG", 1, dim=2)
F = Function(M)
f = Expression(x[0], "x[0]*x[1]")     
F.interpolate(f)  # example for F

u = TrialFunction(V)
u0 = Function(V)
v_u = TestFunction(V)

r and r' are 3d vectors and dV is the volume element.
How can I do this?

yes dV is the volume element of the unit square. This is the first integral for Phi(r) in helmholtz decomposition.
http://en.wikipedia.org/wiki/Helmholtz_decomposition#Statement_of_the_theorem — bshankar, May 1, 2014

MiroK · Answer 1 · May 1, 2014

Hi, the following works but it is quite slow with more reasonable meshes

class PhiExpression(Expression):
    def __init__(self, F):
        Expression.__init__(self)
        mesh = F.function_space().mesh()
        W = FunctionSpace(mesh, 'DG', 0)
        self.divF = project(div(F), W)
        self.r_mag = \
            Expression('sqrt((x[0]-a)*(x[0]-a)+(x[1]-b)*(x[1]-b))', a=0, b=0)

    def eval(self, values, x):
        self.r_mag.a = x[0]
        self.r_mag.b = x[1]
        values[0] = assemble(self.divF/self.r_mag*dx)

mesh = UnitSquareMesh(20, 20)
V = FunctionSpace(mesh, 'CG', 1)
M = VectorFunctionSpace(mesh, 'CG', 1)
F = interpolate(Expression(('x[0]', 'x[0]*x[1]')), M)

phi = interpolate(PhiExpression(F), V)
plot(F, mesh=mesh)
plot(phi)
interactive()

Evaluating an integral like an expression

1 Answer