Euclidean distance map of mesh (distance to surface)

mesh
2

Hi, one option (which works in parallel out of the box) is to get the distance by solving the Eikonal equation, see Johan Hake’s answer here. However, if the boundary is complex there are some difficulties in getting the Newton solver to converge on the nonlinear problem that is the Eikonal equation. You can try to make it more robust following the discussion in the linked post. An alternative presented below is to decompose the boundary, solve simpler problems and then get the final distance as the minimum of the distance fields from the boundary pieces.

from dolfin import *


def distance_from_bdry_piece(facet_f, tag):
    '''Solve Eikonal equation to get distance from tagged facets'''
    mesh = facet_f.mesh()
    V = FunctionSpace(mesh, 'CG', 1)
    bc = DirichletBC(V, Constant(0.0), facet_f, tag)

    u, v = TrialFunction(V), TestFunction(V)
    phi = Function(V)
    f = Constant(1.0)

    # Smooth initial guess
    a = inner(grad(u), grad(v))*dx
    L = inner(f, v)*dx
    solve(a == L, phi, bc)

    # Eikonal equation with stabilization
    eps = Constant(mesh.hmax()/100)
    F = sqrt(inner(grad(phi), grad(phi)))*v*dx - inner(f, v)*dx + eps*inner(grad(phi), grad(phi))*v*dx
    solve(F == 0, phi, bc)

    return phi


def function_min(foos):
    '''Function that is min of foos'''
    f = foos.pop()
    V = f.function_space()

    vec = lambda u: as_backend_type(u.vector()).vec()
    # The work vector
    min_f = f.copy()
    f_vec = vec(min_f)
    # Reduce
    while foos:
        g = foos.pop()
        f_vec.pointwiseMin(f_vec, vec(g))
    # Sync
    as_backend_type(min_f.vector()).update_ghost_values()
        
    return min_f


mesh = UnitSquareMesh(32, 32)
facet_f = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
CompiledSubDomain('near(x[0], 0)').mark(facet_f, 1)
CompiledSubDomain('near(x[0], 1)').mark(facet_f, 2)
CompiledSubDomain('near(x[1], 0)').mark(facet_f, 3)
CompiledSubDomain('near(x[1], 1)').mark(facet_f, 4)

phis = [distance_from_bdry_piece(facet_f, tag) for tag in (1, 2, 3, 4)]
for i, phi in enumerate(phis):
    File('d%d.pvd' % i) << phi
# Final distance is the min of these
phi = function_min(phis)
File('distance.pvd') << phi
1 Like