Finding nodes/edges that are on the boundary of a 3D domain

On Fri, Jun 24, 2011 at 05:35:57PM -0000, Neilen Marais wrote:
> New question #162628 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/162628
>
> Hi,
>
> I'm trying to find all the edges that lie on the boundary of a 3D domain. finding the faces works the normal way:
>
> import dolfin
>
> mesh = dolfin.UnitCube(1,1,1)
> class Boundary(dolfin.SubDomain):
> def inside(self, x, on_boundary):
> print on_boundary
> return on_boundary
> boundary_subdomain = Boundary()

You can use the class DomainBoundary for this:

b = DomainBoundary()
f = FacetFunction('uint', mesh)
f.set_all(0)
b.mark(f, 1)

> face_meshfunction = dolfin.MeshFunction('uint', mesh, 2)
> face_meshfunction.set_all(0)
> boundary_subdomain.mark(face_meshfunction, 1)
>
>
> but if I try
>
> edge_meshfunction = dolfin.MeshFunction('uint', mesh, 1)
> edge_meshfunction.set_all(0)
> boundary_subdomain.mark(edge_meshfunction, 1)
>
> no edges are marked. It seems as if on_boundary is never true for edges or nodes.
>
> In the mean time I'll work aroud it by first finding the faces and then taking their edges, but it would be nice if I could do it in one step.

The reason this doesn't work is that on_boundary is set to true only
for facets since there's no "easy" way to check whether an edge is on
the boundary. For a face, one can just check whether it has 1 or 2
cell neighbors. For an edge, you need to manually check whether it is
part of a facet that is on the boundary.

--
Anders

Hi Anders,

So do I understand correctly:

In order to get the boundary edges I need to

1) Mark all the boundary faces using DomainBoundary()
2) Manually mark the edges by checking that they are connected to the boundary faces (or possibly more efficiently by looping over all the boundary faces and marking their edges)

Thanks
Neilen

On Mon, Jun 27, 2011 at 03:01:21PM -0000, Neilen Marais wrote:
> Question #162628 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/162628
>
> Neilen Marais posted a new comment:
> Hi Anders,
>
> So do I understand correctly:
>
> In order to get the boundary edges I need to
>
> 1) Mark all the boundary faces using DomainBoundary()
> 2) Manually mark the edges by checking that they are connected to the boundary faces (or possibly more efficiently by looping over all the boundary faces and marking their edges)

Yes. And the latter option is more efficient.

--
Anders

Thanks, seems to be working by looping over the boundary faces.