Issue with marking an inner boundary

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.

Setting the eps kwarg in FEM.near() to 1e-3 yields the desired output for me. Maybe consider a mesh size dependent tolerance.

And by the way: Nicely written and very clean code. You don’t see that too often :slight_smile:

…let me advertise vtkplotter.dolfin :slight_smile:

from vtkplotter.dolfin import plot
plot(bmesh, style=1, lc='r', lw=4)
plot(mesh,  wireframe=True, c='k', add=True)

#import matplotlib.pyplot as plt
#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= 1e-3):
#            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()

image

Thank you for your reply. Indeed, drastically relaxing the tolerance does the trick. It’s odd though, I’m still wondering why the “manual” marking (plot-thing) does work with the strict tolerance? Would love to get to the bottom of this

That’s great, thank you for letting us know about it.
I have two questions, though:

  1. Can I choose different colors for different boundaries?
  2. I’m trying to run vtkplotter on a google.colab jupyter notebook, the installation seems to run smoothly but actually trying to plot seems to crash the whole thing, have you tried this?

Can I choose different colors for different boundaries?

As your bmesh is a single object so you cannot do it with the plot() shorcut. Still you can do it with

from vtkplotter.dolfin import plot
from vtkplotter import Points

colors = ('k', 'r')
cols, coords = [], []
for x,y in bmesh.coordinates():
    cols.append(colors[which_circle(x, y, circles)])
    coords.append([x,y,0])
pts = Points(coords, c=cols)
plot(mesh, pts, style=1, wireframe=True, c='k')

image

I’m trying to run vtkplotter on a google.colab jupyter notebook, the installation seems to run smoothly but actually trying to plot seems to crash the whole thing, have you tried this?

I have no experience with that! But it looks from the debug dump , as it is running on a server, you cannot pop up a rendering window and the webGL backend is not working as expected… I will investigate if there is possible fix…

1 Like

Okay, it seems, that by default the optional parameter check_midpoint to the SubDomain.mark() method is set to true. If the midpoint between two vertices does not meet your criterion (which it won’t since the midpoints aren’t exactly on the circle) then also the vertices won’t be recognized as being inside your subdomain. Of course, if you set the tolerance large, at some point also the midpoints will fit the criterion.

Calling the mark method as
active_boundary.mark(meshfn, i+1,False)
in your function mark_boundaries_circle() works even with the small eps you first used.

1 Like

Thanks, @klunkean
that works and it’s nice to know why it does! marking your answer as solution

(for future ref. a good tolerance is around 1e-12)

Cheers!
M

Thanks for your reply. Very good, it will be very interesting to see if/when you get that working

Cheers,
M.

1 Like