How to obtain indices of the mesh points on boundary? (to plot boundary in matplotlib)

Question

I have a mesh with a (roughly) circular boundary.
I want to plot the boundary in matplotlib, so I need an array of points (x,y) corresponding to the points on the boundary.

To obtain the mesh point positions for a vector space V, we can do
n = V.dim() d = u.geometric_dimension() dof_coordinates = Vf.dofmap().tabulate_all_coordinates(M) dof_coordinates.resize((n, d)) x = dof_coordinates[:, 0] y = dof_coordinates[:, 1]

How do I get the indices of these arrays corresponding to points on my boundary?
(My meshes are defined in xml files and vary in lattice topology.)
Thank you in advance!

Answer 1

One question: Assuming that you somehow have the indices you want, how are you planning on making sure that the "boundarynodes-vector" has its components sorted in the correct "order" to get a nice plot?

anyway, if you are only interested in the boundary coordinates, then the following works for me:

In [10]: from dolfin import *
In [11]: mesh = UnitSquareMesh(3,3) 
In [12]: bmesh = BoundaryMesh(mesh, "exterior", True)
In [13]: print bmesh.coordinates()
[[ 0.          0.        ]
 [ 0.33333333  0.        ]
 [ 0.          0.33333333]
 [ 0.66666667  0.        ]
 [ 1.          0.        ]
 [ 1.          0.33333333]
 [ 0.          0.66666667]
 [ 1.          0.66666667]
 [ 0.          1.        ]
 [ 0.33333333  1.        ]
 [ 0.66666667  1.        ]
 [ 1.          1.        ]]

If you are interested in plotting boundary data, then you may want to look at this.

Here is a minimal code. You'd have to view the output with paraview, however.

from dolfin import *

mesh = UnitSquareMesh(3,3)
bmesh = BoundaryMesh(mesh, "exterior", True)

V = FunctionSpace(mesh, "Lagrange", 1)
Vb = FunctionSpace(bmesh, "Lagrange", 1)

u = Function(V)
ub = interpolate(u,Vb)

File("./test_u.pvd") << u
File("./test_ub.pvd") << ub

Comments

Thanks multigrid202. Since my bounaries are all single-valued in polar coordinates (ie, there is a single-valued radial value for the boundary location for each theta), I draw the boundary like this:

import matplotlib.pyplot as plt 

bmesh = BoundaryMesh(mesh, "exterior", True)
Bxytmp = bmesh.coordinates()
Btheta = np.arctan2(Bxytmp[:,1],Bxytmp[:,0])
Binds = np.argsort(Btheta)
Bxy = np.vstack((Bxytmp[Binds,:],Bxytmp[Binds[0],:]))
fig, ax = plt.subplots(1, 1)
ax.plot(Bxy[:,0],Bxy[:,1],'k-')