Defining Dirichlet BC based on level set
I want to solve a transient heat problem on a circular subregion of a rectangular domain. I define level set function to identify the circular interface and apply Dirichlet BC on the nodes which are closer to this interface. I want to assign zero to all nodes outside the circle. I could do this by writing the definition of level set function itself into the outerbound function. But it does not work by return (sign(l)>0.0). Here sign(l) gives a positive value outside of the circle and negative inside. Do you know how can this be done? For this problem there is a solution but I want to use this assignment for a complex problem where level set moves so the boundary changes.
Here is the code:
from fenics import * import time from mshr import * # Time step t_end=1E-1; t=0.0 dt = t_end/50.0 theta = Constant(0.5) TOL=1.E-12 # Generate mesh point=1.1 mesh = RectangleMesh(Point(-point,-point), Point(point,point), 32, 32,'crossed') # Geometry center = Point(0.0,0.0) radius = 1.0 # Elements and space V = FunctionSpace(mesh, "Lagrange", 1) c = TrialFunction(V); c0 = Function(V); l=Function(V); c_ = TestFunction(V) # Level set lexp=Expression("sqrt(x*x + x*x)-rad", rad=radius, degree=1) l.assign(interpolate(lexp,V)) # Define boundaries def interface(x, on_boundary): return (abs((sqrt(x*x + x*x)-radius))<1.0/32) def outerbound(x, on_boundary): return (sign(l)>0.0) # Not working #return (sqrt(x*x + x*x) > radius+1.0/32) # Working bcs = list() bcs.append( DirichletBC(V,Constant(1.0),interface) ) bcs.append( DirichletBC(V,Constant(0.0),outerbound) ) # Problem definition F = (c-c0)/dt*c_*dx+theta*dot(grad(c),grad(c_))*dx+(1.0-theta)*dot(grad(c0),grad(c_))*dx a, L = lhs(F), rhs(F) # Solve problem t = dt while t < t_end: # Solve c=Function(V) solve(a == L, c, bcs) # Update variables for the next time step t += dt c0.assign(c)
Here is what I want to get:
Here is the error message:
TypeError: SWIG director type mismatch expected 'bool' as the output argument of 'inside'
sign()is a UFL function and thus is not numerically evaluated. Second, your levelset function is not in any way connected to the input argument to the
inside()function of your subdomain. The way you wrote it, it is always the
sign()of some dolfin function object rather than the function object's value at position
xwhich you are really interested in.
This will solve the problem:
import numpy as np
and change the definition of your subdomain to
def outerbound(x, on_boundary): return np.sign(l(x))>0.0