Distance to boundary

Question

I am missing a functionality to calculate the distance to the closest boundary.
Google helped me to find this: http://bazaar.launchpad.net/~johan-hake/dolfin/distance-to-boundary

It seems to implement just what I am looking as:

distances = Function(V)

# Compute distance to Boundary for each vertex
distance_computation = DistanceToBoundaryComputation(mesh)
distance_computation.vertex_distances(distances.vector())

Unfortunately, I did not find the functionality in the latest DOLFIN git revision. So I'm wondering if this was ever merged to the main project? Or is this functionality available by some other way now?

Answer 1

Use BoundingBoxTree and BoundaryMesh:

from dolfin import *

mesh = UnitCubeMesh(4,4,4)
bmesh = BoundaryMesh(mesh, "exterior")

bbtree = BoundingBoxTree()
bbtree.build(bmesh)

vertex_distance_to_boundary = MeshFunction("double", mesh, 0)

for v_idx in xrange(mesh.num_vertices()):
    v = Vertex(mesh, v_idx)
    _, distance = bbtree.compute_closest_entity(v.point())
    print distance
    vertex_distance_to_boundary[v_idx] = distance

plot(vertex_distance_to_boundary)
interactive()

Comments

This does not seem to work correctly in 2D (i.e. for UnitSquareMesh).

---

You're right, this seems to be a bug...

---

Ok I'll file an issue to Dolfin in Bibucket. Thanks!
Issue:
https://bitbucket.org/fenics-project/dolfin/issue/139/boundingboxtreecompute_closest_entity

Answer 2

I implemented the distance-to-boundary functionality but found that solving the eikonal equation was much more convenient. For example if you want to get the distance from just a small set of the boundary and the geometry is really complicated and you do not want the direct distance but rather the distance within the domain:

from dolfin import *
import math

installation = "/PATH/TO/DOLFIN/share/dolfin"
demo_path = "/demo/documented/stokes-mini/"

mesh_path = installation + demo_path
mesh = Mesh(mesh_path + "dolfin_fine.xml.gz")

V = FunctionSpace(mesh, 'CG', 1)
v = TestFunction(V)
u = TrialFunction(V)
f = Constant(1.0)
y = Function(V)

class DolphinBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((x[0]-0.5)**2+(x[1]-0.5)**2) < 0.25

boundary = DolphinBoundary()
bc = DirichletBC(V, Constant(0.0), boundary)

# Initialization problem to get good initial guess for nonlinear problem:
F1 = inner(grad(u), grad(v))*dx - f*v*dx
solve(lhs(F1)==rhs(F1), y, bc)

# Stabilized Eikonal equation 
print "max cell size:", mesh.hmax()
eps = Constant(mesh.hmax()/25)
F = sqrt(inner(grad(y), grad(y)))*v*dx - f*v*dx + eps*inner(grad(y), grad(v))*dx
solve(F==0, y, bc)

print "Max distance:", y.vector().max()

plot(y, rescale=True, title="Eikonal", interactive=True)

Comments

Nice, but this does not seem to be robust approach. For example, replace boundary with "on_boundary" such that all boundaries are considered, and then this does not converge with linear elements (but quadratic elements seem to work).

---

I just wanted to comment that this answer helped me a lot in a particular problem, thank you! It is correct as Tuomas wrote that the approach is not very robust, but it can be made very robust by two fairly simple modifications:

1.) Write down the resulting Newton equations yourself (at the PDE level) and implement the iteration yourself, but start with larger stabilization parameter epsilon, i.e. epsilon=1. Decrease epsilon by a factor 0.5 every Newton step. The initial guess u=0 works fine.

This leads to Newton update equations that are of singularly perturbed form, i.e. standard Galerkin discretization is unstable as epsilon becomes small. Thus,

2.) Use an SUPG stabilization in the update equation. By using the fact that |grad(u)| should be about 1 the standard SUPG terms can be greatly simplified, even more when you modify the epsilon to be h*epsilon where h is the mesh size.

If anybody wants me to I can provide a working code example.