(Deleted) Implementation of the reinitialization for level set function


155
views
0
11 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)
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)

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 
F = phi*w*dx - phi0*w*dx - dTau*theta*inner(grad(w),n)*phi*dx - dTau*theta*inner(grad(w),n)*phi0*dx + \
    dTau*inner(grad(w),n)*phi*phi0*dx + dTau*epsilon*theta*inner(grad(w),n)*inner(grad(phi),n)*dx + \
    dTau*epsilon*theta*inner(grad(w),n)*inner(grad(phi0),n)*dx     

# 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
Please login to add an answer/comment or follow this question.
The thread is closed. No new answer/comment may be added.

Similar posts:
Search »
  • Nothing matches yet.