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

Hello,

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)
subdomains.set_all(0)

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

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

x = SpatialCoordinate(submesh)

center = Point(0.5, 0.5)

#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))

#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.
m=1
V = FunctionSpace(submesh, "Lagrange", m)

v = TestFunction(V)

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

#evaluating only u_e on boundary.
)

#printing linear form vector.
b=assemble(L)
print(b[:])

``````

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)
dofxy=V.tabulate_dof_coordinates()
for i,x in enumerate(dofxy):
r = sqrt(pow(x[0]-0.5,2)+pow(x[1]-0.5,2))
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)
)
``````

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)
...
``````

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.

TLDR:
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 = (

)

+ 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)
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)

dofxy=V.tabulate_dof_coordinates()

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

d.vector()[2*i]=p1.x - p2.x
d.vector()[2*i+1]=p1.y - p2.y
u_e.vector()[i]=sin(15*pi*p1.x)*sin(15*pi*p1.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)