This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Newton Raphson for reinitialization on level set method

+1 vote

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.
asked Jul 4, 2017 by pimenta.pvcl FEniCS Novice (320 points)

1 Answer

+1 vote
 
Best answer

In level-set methods normal vector is computed with level-set function,

$$\hat{n} = \frac{\nabla \phi}{|\nabla \phi|}$$

(but normal vector is kept constant during reinitialization)
so it is a function in whole computational domain. FEniCS is trying to tell you that you are
interating something living on Facet over Cells - it doesn't make sense.

answered Jul 4, 2017 by mhabera FEniCS User (1,890 points)
selected Jul 4, 2017 by pimenta.pvcl
...