Defining Dirichlet BC based on level set

8 weeks ago by

Dear all,

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
dt = t_end/50.0      
theta = Constant(0.5)

# Generate mesh
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[0]*x[0] + x[1]*x[1])-rad", rad=radius, degree=1)

# Define boundaries
def interface(x, on_boundary):
    return (abs((sqrt(x[0]*x[0] + x[1]*x[1])-radius))<1.0/32)
def outerbound(x, on_boundary):
    return (sign(l)>0.0)					# Not working
    #return (sqrt(x[0]*x[0] + x[1]*x[1]) > 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
   solve(a == L, c, bcs)

   # Update variables for the next time step
   t += dt

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'


Community: FEniCS Project

1 Answer

8 weeks ago by
Okay, two problems here. First, when using a star import, 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 x which you are really interested in.

This will solve the problem:
Import numpy
import numpy as np​

and change the definition of your subdomain to

def outerbound(x, on_boundary):
    return np.sign(l(x))>0.0
Hi Klunkean,

Thank you, it worked. One other question: is there a way to delete elements for which all the nodes are outside the circle where l(x)>0.

written 8 weeks ago by Birkan  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »