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)
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.
m=1
Fk = FiniteElement("Quadrature", triangle, m+20, quad_scheme='default')
V = FunctionSpace(submesh, "Lagrange", m)
v = TestFunction(V)
a_metadata = {'quadrature_degree': Fk.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.
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()
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.
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 = (
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)
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)
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.