How to replicate an expression based vector field explicitly. i.e. distances to complicated boundary


I am working on an immersed method for the simple poisson problem, which requires that I compute the distances from my computational domain to the true boundary. A simplified, low resolution example is below.

I am successfully able to implement the circle problem using analytical expressions within my bilinear and linear forms. I am trying to do the same thing but with arbitrary geometry (list of segments).

attempt minimum working example

from dolfin import *

#define subdomain with level-set
class Omega1(SubDomain): 
    def inside(self, x, on_boundary): 
        return True if ((x[0]-0.5)*(x[0]-0.5) + (x[1]-0.5)*(x[1]-0.5)) <= 0.09 else False 

#create mesh
mesh = UnitSquareMesh(nx, nx) 

#instantiate submesh  
subdomains = MeshFunction("size_t", mesh, 2) 

subdomain1 = Omega1() 
subdomain1.mark(subdomains, 1) 

#assign marker to boundary   
submesh = SubMesh(mesh, subdomains, 1) 
boundary_parts = MeshFunction("size_t", submesh, 1) 
boundMesh = BoundaryMesh(submesh,'exterior') 
n = FacetNormal(submesh) 
x = SpatialCoordinate(submesh) 

center = Point(0.5, 0.5) 
radius = 0.3  

#create distance function as a list vector, x_mod and y_mod are the locations on the true circle as a function of grid locations x[0] and x[1].
r = sqrt(pow(x[0]-center[0],2)+pow(x[1]-center[1],2)) 

x_mod = x[0]+(x[0]-center[0])*(radius-r)/r 
y_mod = x[1]+(x[1]-center[1])*(radius-r)/r 

#create u_e expression which is the solution at the closest point on the true circle projected onto some coordinate within the domain.
u_e = sin(15*pi*x_mod)*sin(15*pi*y_mod) 

#setting up only linear form for simplicity, m is polynomial order.
Fk = FiniteElement("Quadrature", triangle, m+20, quad_scheme='default')             
V = FunctionSpace(submesh, "Lagrange", m) 

v = TestFunction(V) 
a_metadata = {'quadrature_degree':, 'quad_scheme': Fk.quadrature_scheme()} 

ds = Measure('ds', domain=submesh,  subdomain_data=boundary_parts, metadata=a_metadata) 

f = 450*pi*pi*sin(15*pi*x[0])*sin(15*pi*x[1])

#evaluating only u_e on boundary.
L = (f*v*dx(metadata=a_metadata) 
     + dot(grad(v),n)*u_e*ds(1) 

#printing linear form vector.

This approach is good for geometry I can analytically express. However I will need to extend this to 2d geometry defined by a list of an arbitrary amount of segments, and I don’t clearly see how I can use expressions.

I attempted to explicitly assign the nodal values of u_e into a vector. Like:

u_e = Function(V)
radius = 0.3
for i,x in enumerate(dofxy):
    r = sqrt(pow(x[0]-0.5,2)+pow(x[1]-0.5,2)) 
    x_mod = x[0]+(x[0]-0.5)*(radius-r)/r 
    y_mod = x[1]+(x[1]-0.5)*(radius-r)/r 
    u_e.vector()[i] = sin(15*pi*x_mod)*sin(15*pi*y_mod)

Which then went into the same linear form:

L = (f*v*dx(metadata=a_metadata) 
     + dot(grad(v),n)*u_e*ds(1) 

assembling L and printing showed that these approaches are not the same. My naive guess was that this had to do with the interpolation to the quadrature points. leading me to try:

Vt = FunctionSpace(submesh, "Lagrange", m)
u_ev = interpolate(u_e, Vt)
    L = (f*v*dx(metadata=a_metadata) 
#         + dot(grad(v),n)*u_e*ds(1) 
         + dot(grad(v),n)*u_ev*ds(1) 

This got me closer to the L assembled with the expression based u_e but I began losing optimal convergence at high order polynomials ie m>=4 relative to the manufactured solution. I achieve optimal convergence for the expression based u_e and d.

The larger picture here beyond the minimal example is that I am hoping someone can enlighten me on how to implement a distance vector field d (closest point projection) and projected solution u_e that is equivalent to the expression based approach. e.g. a square created by a list of 4 segments.

with the full system being:

        a = (
            inner(nabla_grad(u), nabla_grad(v))*dx(metadata=a_metadata) 
             - dot(grad(u),n)*v*ds(1) 
             + dot(grad(v),n)*(u+dot(grad(u),d))*ds(1)  #d ListTensor located here.


        L = (f*v*dx(metadata=a_metadata) 
             + dot(grad(v),n)*u_e*ds(1) #u_e located here

Where d is a ListTensor formed by

distx = ((x[0]-center[0])/r)*(radius-r) 
disty = ((x[1]-center[1])/r)*(radius-r) 
d = as_vector([distx,disty]) 

which is understood by fenics/ufl as:

[(0.3 + -1 * sqrt((-0.5 + x[0]) ** 2 + (-0.5 + x[1]) ** 2)) * (-0.5 + x[0]) / sqrt((-0.5 + x[0]) ** 2 + (-0.5 + x[1]) ** 2), (0.3 + -1 * sqrt((-0.5 + x[0]) ** 2 + (-0.5 + x[1]) ** 2)) * (-0.5 + x[1]) / sqrt((-0.5 + x[0]) ** 2 + (-0.5 + x[1]) ** 2)]

Right now I am building a vector function space like:

    V = FunctionSpace(submesh, "Lagrange", m) 
    Vd = VectorFunctionSpace(submesh, "Lagrange", m) 

    d = Function(Vd)
    u_e = Function(V)

    for i in range(len(dofxy)):
        p1, p2 = ops.nearest_points(square.boundary,

        d.vector()[2*i]=p1.x - p2.x
        d.vector()[2*i+1]=p1.y - p2.y

which appears to be fine up until order 4 polynomials for a square.

Even more TLDR:

How could I replicate this:

distx = ((x[0]-center[0])/r)*(radius-r) 
disty = ((x[1]-center[1])/r)*(radius-r) 
d = as_vector([distx,disty]) 

for a square (or poly), and not a circle.

If there is a solution that also uses expressions, that would be fine as well.