### (Deleted) Implementation of the reinitialization for level set function

181
views
0
13 months ago by
Dear all,

Once again thanks for all great help that you folks provide.

I am trying to implement the reinitialization level function: $\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}]$, which has the finite element discretization witting as $\frac{1}{\Delta \tau} \left \langle (\phi^{n+1} - \phi^{n}),w \right \rangle - \frac{1}{2}\left \langle \nabla w,(\phi^{n+1} + \phi^{n})\hat{n} \right \rangle + \left \langle \phi^{n+1}\phi^{n}\hat{n}, \nabla w\right \rangle + \frac{\varepsilon}{2}\left \langle( \nabla \phi^{n+1} + \nabla\phi^{n})\cdot \hat{n},\nabla w \cdot \hat{n} \right \rangle = 0$ and my code for this is:

# Modules
from fenics import*
import time
%matplotlib inline

theta = Constant(0.5)    # theta schema

#   given mesh and function l on the mesh
mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 4.0), 160, 160,"crossed")
F = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
V = VectorFunctionSpace(mesh, "DG", 0, dim=2)
FE = FunctionSpace(mesh,F)

# Initial condition for the level set function
center = Point(0.5, 0.5)
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)

phi = TrialFunction(FE)
w = TestFunction(FE)
phi0 = interpolate(dist, FE)

# Normal vector
n = project(grad(phi0)/sqrt(phi0.dx(0)**2 + phi0.dx(1)**2), V)

# Time and epsilon parameters
Delta_x = 1.0/160.0
P1 = 0.5*(pow(Delta_x,0.9))
P2 = 0.5*(pow(Delta_x,1.1))
epsilon = Constant(P1)
dTau = Constant(P2)

# Function

# Bilinear form
a, L = lhs(F), rhs(F)

bc = []
phi = Function(FE)
num_steps = 20                            # number of time steps

# Solver
for n in range(num_steps):
# Compute solution
solve(a == L, phi, bc)               # Relatinship between
x = [0.0, 0.5]
error = phi(x) - phi0(x)
print error
phi0.assign(phi)                     # Update the previous solution #
plot(phi)
interactive()    ​
The problem is I am not able to get the right solution. Thus, if someone with expertise in level set function may help me I would appreciate that.
Community: FEniCS Project