integrating over variable domain

Question

Given a function $f$, I'm interested in defining another function that gives the average value of $f$ in a ball of a fixed radius. Namely, I want to define

$$g(x) = \int_{B_r(x)}f(y)dy,$$

where $B_r(x)$ is a ball of radius $r$ centered at $x$. I currently don't have any ideas about how to do this since it involves integrating over a subdomain that varies---that is, depends on $x$. If anyone has an idea about how to implement this, please let me know.

Answer 1

Hi, the following snippet uses an instance of SubDomain class to mark cells of the mesh that are in inside the circle with flag 1. Assembler then uses only these cells to get the approximate area.

from dolfin import *

class Circle(SubDomain):
    'Interior of circle with radius r centered at (x_c, y_c).'
    def __init__(self, x_c, y_c, r):
        self.x_c = x_c
        self.y_c = y_c
        self.r = r
        SubDomain.__init__(self)

    def inside(self, x, on_boundary):
        return (x[0]-self.x_c)**2 + (x[1]-self.y_c)**2 < self.r**2

mesh = UnitSquareMesh(100, 100)
cell_f = CellFunction('size_t', mesh, 0)
dx = Measure('dx')[cell_f]

x0, y0, r = 0., 0., 0.125
circle = Circle(x0, y0, r)

def displace(circle, step):
    if step < 5 or step >= 11:
        circle.x_c += 0.1
        circle.y_c += 0.1
    elif step < 8:
        circle.r += 0.1
    elif step < 11:
        circle.r -= 0.1

t = 0
for step in range(17):
    circle.mark(cell_f, 1)
    plot(cell_f)
    print 'area', assemble(1*dx(1, domain=mesh))
    displace(circle, step)
    # Set all tags to zero for fresh start in next iteration
    cell_f.set_all(0)

interactive() 

Comments

Thank you for your answer, but it isn't clear to me how what you wrote can be used to directly answer my question. You clearly showed how to specify a subdomain that can be moved around, but how can this be used to define the function $g$ mentioned in my question? Perhaps I'm missing something obvious.

---

I was able to use what you wrote to come up with the following code which solves the problem. However, this solution is very time consuming as it involves defining the function at each mesh point. If you know how to do this without cycling through all of the mesh points, please let me know.

from dolfin import *

class Circle(SubDomain):
    'Interior of circle with radius r centered at (x_c, y_c).'
    def __init__(self, x_c, y_c, r):
        self.x_c = x_c
        self.y_c = y_c
        self.r = r
        SubDomain.__init__(self)

    def inside(self, x, on_boundary):
        return (x[0]-self.x_c)**2 + (x[1]-self.y_c)**2 < self.r**2

mesh = UnitSquareMesh(20, 20)
cell_f = CellFunction('size_t', mesh, 0)
dx = Measure('dx')[cell_f]

V = FunctionSpace(mesh, "CG", 1)

# function to be integrated
f = Expression('sin(x[0]+x[1])')

# function to be defined
g = Function(V)
g_array = g.vector().array()

# define g at each mesh point
coor = mesh.coordinates()
for ii in range(mesh.num_vertices()):
    circle = Circle(coor[ii][0], coor[ii][1], 0.125)
    circle.mark(cell_f, 1)
    g_array[ii] += assemble(f*dx(1, domain=mesh))
    cell_f.set_all(0)

g.vector()[:] = g_array
plot(g)