Dear Community

Please find below the code segment which solves the 2D axisymmetric conduction problem in a spherical annulus.

There is a parameter called “TOL” in the program which is used to pick out points on the boundary. I have seen examples where TOL values as low as 1e-14 have been used. However, for the current set of parameters shown in the code, the lowest I can go is TOL=1e-3, and the computed solution agrees well with the analytical result.

However, if I use TOL=1e-4, I get the message

*** Warning: Found no facets matching domain for boundary condition" ***

and the computed solution is not correct. Using a higher number of segments for the polygonal approximation of a circle (the “num_segments” parameter) results in the same behavior. Using a higher value of the “mesh_resolution” parameter also does not help, and at mesh_resolution=512, I get the following error,

UMFPACK V5.7.8 (Nov 9, 2018): ERROR: out of memory

along with the previous warning about no facets found.

Could you please guide me on this matter ? Specifically, what should I do in order to be able to set a lower TOL value? I am assuming that smaller the value of TOL, more accurately would the boundaries be picked out.

Thank You

Warm Regards

```
from dolfin import *
from mshr import *
import math
import numpy as np
Re = 5.
Ri = 1.
num_segments=500
domain = Circle(Point(0., 0.), Re, num_segments) - Circle(Point(0., 0.), Ri, num_segments)
mesh_resolution = 64
mesh = generate_mesh(domain, mesh_resolution)
V = FunctionSpace(mesh, "CG", 2)
x = SpatialCoordinate(mesh)
r = Expression("sqrt(x[0]*x[0]+x[1]*x[1])", degree=2)
theta = Expression("atan2(x[1],x[0])", degree=2)
# Define boundary subdomains
TOL = 1e-3
# TOL = 1e-4 or lower results in
# "Warning: Found no facets matching domain for boundary condition"
u_inner = Expression('cos(theta)', degree=2, theta=theta)
def inside(x, on_boundary):
return near(x[0]**2 + x[1]**2, Ri**2, TOL) and on_boundary
u_outer = Constant(0.0)
def outside(x, on_boundary):
return near(x[0]**2 + x[1]**2, Re**2, TOL) and on_boundary
inner_circ = DirichletBC(V,u_inner,inside)
outer_circ = DirichletBC(V,u_outer,outside)
bcs = [inner_circ, outer_circ]
u = TrialFunction(V)
v = TestFunction(V)
a1 = (u.dx(0)*v.dx(0))+(u.dx(1)*v.dx(1))
a2 = (Constant(1.0)/(r*sin(theta)))*v*(u.dx(1))
a = (a1-a2)*dx
f = Constant(0.0)
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bcs)
```