Dear all,
I am trying to implement a reinitialization procedure for level set method using the following equation:
$$\phi_{\tau}(\vec{x},\tau) + \nabla\cdot[\phi(\vec{x},\tau)(1-\phi(\vec{x},\tau))\hat{n}] =\varepsilon\nabla\cdot[\hat{n}\phi(\vec{x},\tau)\cdot\hat{n}]$$
# Modules
from fenics import*
import time
%matplotlib inline
# Parameters
theta = Constant(0.5) # theta schema
t_end = 5.0
num_steps = 10
dt = t_end/num_steps
inv_dt = Constant(1.0/dt)
# Mesh
mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 4.0), 20, 80,"crossed")
center = Point(0.5, 0.5)
radius = 0.2
dist = Expression('sqrt((x[0]-A)*(x[0]-A) + (x[1]-B)*(x[1]-B))-r',degree=2, A=center[0], B=center[1],r=radius)
# Difine function space
#L = FiniteElement("Lagrange", mesh.ufl_cell(), 2) # levelset field
V = FunctionSpace(mesh, "Lagrange", 1)
#V = FunctionSpace(mesh,L)
phi = TrialFunction(V)
phi0 = TrialFunction(V)
w = TestFunction(V)
# Initial condition
phi0 = interpolate(dist,V)
n = FacetNormal(mesh)
h = CellSize(mesh)
epsilon = 2*h
psi = -phi*(1-phi)*inner(grad(w),n) + epsilon*(inner(grad(phi),n)*inner(grad(w),n))
psi0 = -phi0*(1-phi0)*inner(grad(w),n) + epsilon*(inner(grad(phi0),n)*inner(grad(w),n))
F = inv_dt*(phi-phi0)*w*dx + theta*psi*dx + (1-theta)*psi0*dx
a, L = lhs(F), rhs(F)
A = assemble(a)
b = assemble(L)
solver = KrylovSolver("cg", "ilu")
solver.set_operator(A)
sol = Function(V)
But I stack since the following errors came out
Integral of type cell cannot contain a FacetNormal.