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.
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*x+x*x)", degree=2) theta = Expression("atan2(x,x)", 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**2 + x**2, Ri**2, TOL) and on_boundary u_outer = Constant(0.0) def outside(x, on_boundary): return near(x**2 + x**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)