Unable to set low tolerance for picking points on boundary, [Warning: Found no facets matching domain for boundary condition]

mathematics
1

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)
2

For your current problem, where you have two external boundaries, namely the circle with radius Ri and the circle with radius Re, I would suggest the following approach:


r_mid = Ri + (Re - Ri) / 2

class Inside(SubDomain):
    def inside(self, x, on_boundary):
        return sqrt(x[0]**2 + x[1]**2) < r_mid  and on_boundary

class Outside(SubDomain):
    def inside(self, x, on_boundary):
        return sqrt(x[0]**2 + x[1]**2) > r_mid and on_boundary

sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(0)
Inside().mark(sub_domains, 1)
Outside().mark(sub_domains, 2)
File("facets.pvd") << sub_domains


inner_circ = DirichletBC(V, u_inner, sub_domains, 1)
outer_circ = DirichletBC(V, u_outer, sub_domains, 2)

as you can see by opening the file facets.pvd in Paraview, you observe that the boundaries are correctly marked, independent of the mesh size.

3

Thank you very much, Dokken, for the prompt response. Your suggestion works perfectly, and I was able to follow the logic of the script.

Warm Regards