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.
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.
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**2 + x**2) < r_mid and on_boundary and x > tol u_outer = Constant(0.0) class Outside(SubDomain): def inside(self, x, on_boundary): return sqrt(x**2 + x**2) > r_mid and on_boundary and x > 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