Picking boundaries for a semi-circular domain

mesh
1

Dear Community

I followed Dokken’s recommendation for correctly picking the points on the boundary of a 2D spherical annulus. The method works perfectly for a circular domain.

I am now working with a semi-circular domain, and the only two boundaries that I want selected are the outer and the inner circumferences. The horizontal bottom portion must not be picked as a boundary. I came up with the following workaround (code shown below) and a visualization of the semi-circular domain is provided. While it appears to be correct, a zoomed in view near the inner circle, shows that a small portion of the circumference, near the bottom, is left unmarked.

entire_domain
entire_domain1764×938 287 KB

zoomed_in
zoomed_in1764×938 169 KB

A consequence of this small portion being left out is that when I try to print the circumference of the inner semi-circle (of unit radius), I am getting 3.01543 as the answer, which is not quite the same as the expected answer, which is \pi.

A similar test applied to the full inner circle, using the approach in, gives the answer 6.282151, which is much closer to the correct answer of 2\pi.

Could you please suggest an alternative for correctly picking the points on the boundary for a semi-circle? I will soon be calculating quantities based on the integral over the inner circumference, and would need to do it accurately.

Thank You
Warm Regards

from dolfin import *
from mshr import *
import math
import numpy as np

Re = 10.
Ri = 1.
num_segments=100
domain = Circle(Point(0., 0.), Re, num_segments) - Circle(Point(0., 0.), Ri, num_segments)

#creating semi-circular domain
domain = (domain - Rectangle(Point(0., 0.), Point(Re, -Re))
                - Rectangle(Point(0., 0.), Point(-Re, -Re)))


mesh_resolution = 128
mesh = generate_mesh(domain, mesh_resolution)
V = FunctionSpace(mesh, "CG", 2)
x = SpatialCoordinate(mesh)

### Defining boundary conditions 

r_mid = Ri + (Re - Ri) / 2
tol = 1e-16

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

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

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)

### writing facets to file
File("facets.pvd") << sub_domains

#### Test of work-around
ds = Measure('ds', domain=mesh, subdomain_data=sub_domains)
print(assemble(Constant(1)*ds(1))) # Expect this to return 
                                   # circumference of semi-circle
2

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
3

Thank you for the suggestion, @conpierce8. I am trying to get this done through mshr as much as possible. If it doesn’t work, I will shift to the Gmsh option.

I tried the suggestion given above, and it seems that the subdomains are not getting marked properly (as seen from the ParaView visualization), and also because the circumference is printed to be “0”.

Using a factor of 1 (instead of 1.5) in the expression for the tolerance does not change the result.

Wanted to update the thread with this observation.

u_inner = Constant(1.0)
class Inside(SubDomain):
    def inside(self, 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

u_outer = Constant(0.0)
class Outside(SubDomain):
    def inside(self, 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


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)

ds = Measure('ds', domain=mesh, subdomain_data=sub_domains)
print(assemble(Constant(1)*ds(1))) # returns circumference of internal circle
4

Hi @rkailash,

You are right. I can see that I made an error in my code. The tol I suggested is applicable to the radius, not the square of the radius. The following code should give the result you desire:

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)

class Inside(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1.5*Ri * (1 - np.cos(np.pi / num_segments))
        return near(sqrt(x[0]**2 + x[1]**2), Ri, tol) and on_boundary

class Outside(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1.5*Re * (1 - np.cos(np.pi / num_segments))
        return near(sqrt(x[0]**2 + x[1]**2), Re, tol) 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)

ds = Measure('ds', domain=mesh, subdomain_data=sub_domains)
print(assemble(Constant(1)*ds(1))) # returns circumference of internal circle
print(assemble(Constant(1)*ds(2))) # returns circumference of internal circle

### writing facets to file
File("facets.pvd") << sub_domains

domains

2 Likes
5

Thank you @conpierce8 , the method suggested above works perfectly!

1 Like