Distance to Boundary

Question

Hi,

I am trying to calculate the distance of interior points w.r.t. boundary. I found a very useful solution DISTANCE TO BOUNDARY, the solution I am using is,

from dolfin import *

mesh = UnitSquareMesh(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()

I am however, having the problem, in case of when not all edges defining the mesh are not included in distance calculation.

For example, in the code above, I am calculating the distance w.r.t. all edges on UnitSquareMesh (Left, Right, Top, Bottom),

                    TOP
                -------------
      LEFT   |                | RIGHT
             |                |
              -------------
                 BOTTOM

but how do I calculate the distance, for instance, only w.r.t. Left, Right EDGES

                   < NOT INCLUDED >
                -------------
      LEFT   |                | RIGHT
             |                |
              -------------
                 < NOT INCLUDED >

Answer 1

Hi, consider the following modification

from dolfin import *

mesh = UnitSquareMesh(40, 40, 'crossed')
bmesh = BoundaryMesh(mesh, 'exterior')
# Suppose now bottom boundary is not of interest
cell_f = CellFunction('size_t', bmesh, 0)
CompiledSubDomain('near(x[1], 0)').mark(cell_f, 1)
bmesh_sub = SubMesh(bmesh, cell_f, 0)
tree = bmesh_sub.bounding_box_tree()

V = FunctionSpace(mesh, 'CG', 1)
v_2_d = vertex_to_dof_map(V)
bdry_distance = Function(V)
values = bdry_distance.vector().array()
for index, vertex in enumerate(vertices(mesh)):
    _, d = tree.compute_closest_entity(vertex.point())
    values[v_2_d[index]] = d
bdry_distance.vector().set_local(values)
bdry_distance.vector().apply('insert')

plot(bdry_distance)
interactive()

Note that due to SubMesh the code does not work in parallel.

Comments

Thanks for the answer.
I am using fencis to kind of check my implementations. Solutions agree in cases where, the geometry is simple, for instance in channel (rectangular, unit-square), But on some not so complicated meshes, like L-shape or Backward-facing-step, I realized that the fenics solution computed Euclidean_distance between interior points and the boundary (points), but how can I find normal distance from the boundary for interior points. Any help would be great.