Hi @rkailash,
The issue here is the way mshr treats circles: not as circles, but as N-sided polygons (N-gons). You have specified num_segments = 500, which means that your mesh is of a 500-gon. mshr places vertices at the endpoints of each boundary segment, and at the midpoint. The vertex at the midpoint is not located on the circle. (In fact, it is displaced inwards in the radial direction by a distance of R(1-\cos{(\pi/N)}). This is why, in your original post, tol = 1e-3 worked but tol = 1e-4 did not. (Consider that for Re = 5 and num_segments = 500, R(1-\cos{(\pi/N)}) = 9.8\cdot 10^{-5}, very close to your tol.)
One way to fix this issue is to use Gmsh rather than mshr to generate your mesh. Gmsh treats circles as circles, not as N-gons, so all vertices on the boundary will be located on the circle, not displaced inwardly as in the mshr-generated mesh. If you take this appoach, your original functions
def inside(x, on_boundary):
return near(x[0]**2 + x[1]**2, Ri**2, TOL) and on_boundary
def outside(x, on_boundary):
return near(x[0]**2 + x[1]**2, Re**2, TOL) and on_boundary
should correctly identify the boundary nodes, even for very small values of TOL.
If you wish to continue using mshr, first consider that this statement from your original post is not strictly true:
The accuracy of picking out the boundaries is binary: the boundaries are correctly picked out, or they are not. The optimal value of tol is the largest value that is sufficiently small to exclude any internal facets. Choosing tol to be any smaller than this value will not have any effect on whether internal facets are picked, but it will make it more likely that external facets will be erroneously excluded. Using the expression that I have provided above, you could start with the following, and see if tol is sufficiently small to exclude all internal facets:
def inside(x, on_boundary):
tol = 1.5*Ri * (1 - np.cos(np.pi / num_segments))
return near(x[0]**2 + x[1]**2, Ri**2, tol) and on_boundary
def outside(x, on_boundary):
tol = 1.5*Re * (1 - np.cos(np.pi / num_segments))
return near(x[0]**2 + x[1]**2, Re**2, tol) and on_boundary