Integrating over a sphere within the mesh

variational formulation
1

Hi there,

I’m just wondering; I have a function defined over the whole of my mesh, which is large. I want to compute the simple integral:

\int_S f dS

over a sphere of radius R centred in the middle of my domain, where R is much smaller than the system. Is there any way to do this? I should say that I do not have this defined as a subdomain.

Best,
Ryan

2

Hi, consider

from dolfin import *

# The idea here is to exploit the divergence theorem
# 
#   \int_{surface of bounded volume} vec.n dS = \int_{bounded volume} div(v) dx
#
# For circle centered at (0, 0) we have n = (x, y)/radius. Now suppose that
# vec = f*(x, y). Then
#
#   \int_{surface of bounded volume} vec.n dS = radius*\int_{surface of bounded volume} f dS


radius = 0.2
mesh = Mesh('domain.xml')

# Conforming mesh
# Here we load the cell function with marked circle (tag 1) of radius 0.2
volumes = MeshFunction('size_t', mesh, 'domain_physical_region.xml')
# If not given, the cell function would be done manually
for cell in cells(mesh):
    if cell.midpoint().norm() < radius:
        assert volumes[cell] == 1  # volumes[cell] = 1

        
def circle_integral(f, vols=volumes, rad=radius):
    # The circle is marked as 1
    dV = Measure('dx', domain=vols.mesh(), subdomain_data=vols, subdomain_id=1)

    x = SpatialCoordinate(mesh)
    value = assemble(div(f*x)*dV)
    # Take the radius into account
    value /= rad
    
    return value

# Sanity 
f = Constant(1)
value = circle_integral(f)
print(value, 2*pi*radius*f(0), abs(value - 2*pi*radius*f(0)))

# Something more advanced
x = SpatialCoordinate(mesh)
f = inner(x, x)
value = circle_integral(f)
print(value, 2*pi*radius**3, abs(value - 2*pi*radius**3))

# Finally chech on nonconforming mesh
in_circle = CompiledSubDomain('x[0]*x[0] + x[1]*x[1] < rad', rad=radius**2)

true = 2*pi*radius**3
e0, h0 = None, None

for i in range(3, 12):
    ncells = 2**i
    mesh = RectangleMesh(Point(-1, -1), Point(1, 1), ncells, ncells, 'crossed')
    volumes = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)

    in_circle.mark(volumes, 1)

    x = SpatialCoordinate(mesh)
    f = inner(x, x)
    value = circle_integral(f, volumes, radius)

    e, h = abs(value - true), mesh.hmin()
    if e0 is not None:
        rate = ln(e/e0)/ln(h/h0)
    else:
        rate = -1.

    e0, h0 = e, h
    print('error %g(%.2f)' % (e, rate))

The Gmsh file used is as follows

// to refine `gmsh -2 domain.geo -clscale 0.5`
// dolfin-convert domain.msh domain.xml

size = 0.1;
len = 1;
rad = 0.2;

// Quarter domain
Point(1) = {0, 0, 0, size};
Point(2) = {len, 0, 0, size};
Point(3) = {len, len, 0, size};
Point(4) = {0, len, 0, size};

// Circle 
Point(5) = {rad, 0, 0, size};
Point(6) = {0, rad, 0, size};

Line(1) = {1, 5};
Line(2) = {5, 2};
Line(3) = {2, 3};
Line(4) = {3, 4};
Line(5) = {4, 6};
Line(6) = {6, 1};

Circle(7) = {5, 1, 6};

Line Loop(1) = {7, 6, 1};
Plane Surface(1) = {1};
Physical Surface(1) = {1}; // Inside circle


Line Loop(2) = {5, -7, 2, 3, 4};
Plane Surface(2) = {2};
Physical Surface(2) = {2}; // Outside

// The rest by symmetry
vols[] = Symmetry {1, 0, 0, 0} {Duplicata { Surface{2}; Surface{1}; }};
Physical Surface(1) += {vols[1]};
Physical Surface(2) += {vols[0]};

vols[] = Symmetry {0, -1, 0, 0} {Duplicata { Surface{vols[0]}; Surface{vols[1]}; }};
Physical Surface(1) += {vols[1]};
Physical Surface(2) += {vols[0]};

vols[] = Symmetry {-1, 0, 0, 0} {Duplicata { Surface{vols[0]}; Surface{vols[1]}; }};
Physical Surface(1) += {vols[1]};
Physical Surface(2) += {vols[0]};