I have a 2D mesh which consists of a square with several disk-shaped holes. I want to mark the boundaries according to whether a point lies on the square’s boundaries (id=0) or the i-th hole boundary (id=i,…)
I got this to work for polygon-shaped holes but I can’t figure out why a similar approach is not working for the circular boundaries
Here’s a MWE
import numpy as np
import dolfin as FEM
import matplotlib.pyplot as plt
# MESHING A UNIT SQUARE FROM WHIC HOLES WILL BE SUBTRACTED
domain = mshr.Rectangle(FEM.Point(0., 0.,), FEM.Point(1., 1.))
# THE LIST OF (RADIUS, X_CENTER, Y_CENTER) DATA
circles = [[0.1, 0.5, 0.5]]; circles = np.array(circles)
# PERFORM THE BOOLEAN DIFFERENCES
for circle in circles:
domain -= mshr.Circle(FEM.Point(circle[1:]), circle[0], 60)
# GENERATE MESH
mesh = mshr.generate_mesh(domain, 30)
# A CLASS FOR ID'ING AN INNER CIRCULAR BOUNDARY
class InnerBoundary_circle(FEM.SubDomain):
def __init__(self, r, x0, y0):
self.r = r
self.x0 = x0
self.y0 = y0
super().__init__()
def on_circle_condition(self, x, y):
return np.power(x-self.x0, 2) + \
np.power(y-self.y0, 2) - \
np.power(self.r, 2)
def point_belongs_to_circle(self, x):
x_, y_ = x[:]
if FEM.near(self.on_circle_condition(x_, y_), 0., eps= 3.e-15):
return True
else: return False
def inside(self, x, on_boundary):
if not on_boundary:
return on_boundary
else:
return self.point_belongs_to_circle(x)
# A FUNCTION TO ID ALL BOUNDARIES
def mark_boundaries_circle(meshfn, circle_data):
for i, circle in enumerate(circle_data):
r, x, y = circle[:]
active_boundary = InnerBoundary_circle(r, x, y)
active_boundary.mark(meshfn, i+1)
print('THE SET OF INDICES IN MESHFN IS', set(meshfn.array()))
print('IT SHOULD BE', set(range(len(circle_data)+1)))
assert len(set(meshfn.array())) == len(circle_data) + 1
# ON TO MARK THE BOUNDARIES
meshfn = FEM.MeshFunction("size_t", mesh, mesh.topology().dim()-1)
meshfn.set_all(0)
try:
mark_boundaries_circle(meshfn, circles) # FAILS
print('MARKING SUCCESFUL!')
except:
print('MARKING UNSUCCESFUL...')
# HOWEVER, THIS SEEMS TO SHOW THAT MARKING SHOULD BE SUCCESFUL...
bmesh = FEM.BoundaryMesh(mesh, 'exterior')
def which_circle(x, y, circles_data):
for i, circle in enumerate(circles_data):
r, x0, y0 = circle
if FEM.near(np.power(x-x0, 2) + np.power(y-y0, 2) - np.power(r, 2), 0., eps= 3.e-15):
return i+1
return 0
with plt.style.context('seaborn-bright'):
colors = ('k', 'r')
fig, ax = plt.subplots(figsize= (10, 10))
FEM.plot(mesh)
for i, (x, y) in enumerate(bmesh.coordinates()):
plt.scatter((x,), (y,), color= colors[which_circle(x, y, circles)])
plt.xlabel('x', fontsize = 'xx-large')
plt.ylabel('y', fontsize = 'xx-large', rotation='horizontal')
ax.set_aspect('equal')
plt.show()
This produces the output:
THE SET OF INDICES IN MESHFN IS {0}
IT SHOULD BE {0, 1}
MARKING UNSUCCESFUL...
and the plot which shows that the boundary points can be identified:
I will appreciate your thoughts on why my InnerBoundary class and my mark_boundaries_circle function don’t seem to work
Using fenics 2019.1.0 installed from the ppa on a google.colab notebook
Thank you,
M.
