Reinitialization — Prediction-Correction method

class levelSet_tool.Reinitialization(level_set, V_ls, l)

Bases: LevelSet

This is the Reinitialization class.

In the subsection we give some details about the Reinitialization method used. We use the Prediction Correction scheme proposed in ….. Initialization of predictor problem.

Definition of the Predictor variational problem

prediction problem is given by: Find \(\phi_{p}:D \rightarrow \mathbb{R}\)

\[\begin{split}\begin{cases} -\Delta\phi_{p} & \!\!\!\! = \psi(x)\text{ in D}\\ \phi_{p} & \!\!\!\! = 0\text{ on }\Gamma \\ \nabla \phi_{p} \cdot n & \!\!\!\! = \psi(x) \text{ on } \partial D . \end{cases}\end{split}\]

With Nitsche method, this yields to the following weak formulation: Find \(\phi^{0} \in V\), such that for all \(v \in V\)

\[a\left(\phi_{p},v\right)=l\left(v\right),\]

where:

\[\begin{split}\begin{align} a\left(\phi_{p},v\right) &= \int_{D}\nabla \phi_{p}\cdot\nabla v\,\text{ }dx{\color{red}{ -\int_{\Gamma}\nabla \phi_{p}\cdot n_{\Gamma} \, v\text{ }ds-\int_{\Gamma}\nabla v\cdot n_{\Gamma} \, \phi_{p}\text{ }ds}} \\ &{\color{red}{ +\gamma_{D}\int_{\Gamma} \phi_{p} \, v\text{ }ds}}\\ l\left(v\right) &= \int_{D}\psi(x) v\text{ }dx + \int_{\partial D}\psi(x) v\text{ }ds \end{align}\end{split}\]

with \(\gamma_{D}>0\) is the Nistche parameter.

The text in red is the Nitsche terms.

Intialization of normal field to isocontour:

V_DG = fem.functionspace(self.mesh, ("DG", 0, (self.dim,)))
self.n_K = fem.Function(V_DG)
self.norm_euclidienne = ufl.sqrt(inner(ufl.grad(self.level_set),ufl.grad(self.level_set)))
self.n_K = ufl.grad(self.level_set)/ self.norm_euclidienne

Bilinear form (Predictor):

import ufl

u_r = ufl.TrialFunction(self.V_ls) # Trial function
v_r = ufl.TestFunction(self.V_ls) # Test function

self.gamma_r = 1e+4 # Value of Nitsche parameter

self.h_r = CellDiameter(self.mesh) # mesh size

self.a_predict  = ufl.inner(grad(u_r), grad(v_r))*self.dx
self.a_predict += - dot(grad(u_r), self.n_K)*v_r*self.dsq
self.a_predict += - dot(grad(v_r), self.n_K)*u_r*self.dsq
self.a_predict += self.gamma_r*1.0/self.h_r*u_r*v_r*self.dsq

Linear form (Predictor):

Approximation of the signed indicator function \(\psi\):

self.eps = 1e-6
self.sign = self.level_set / (ufl.sqrt(self.level_set**2+self.eps**2))

self.L_predict = inner(self.l**2*self.sign,v_r)*self.dx
nuemann_bc = ufl.conditional(ufl.le(self.sign,0),-1,0)
nuemann_bc = ufl.conditional(ufl.ge(self.sign,0),1,nuemann_bc)
self.L_predict += ufl.dot(nuemann_bc,v_r)*ufl.ds

Be aware to impose correct Nuemann condition to guaranty the order of convergence.

Definition of the Corrector variational problem

prediction problem is given by: Find \(\phi:D \rightarrow \mathbb{R}\)

\[\begin{split}\begin{cases} \nabla\cdot\left(\nabla\phi-\frac{\nabla\phi}{\left|\nabla\phi\right|}\right)& \!\!\!\! =0\text{ in }D\\ \phi&\!\!\!\! =0\text{ on }\Gamma\\ \left(\nabla\phi-\frac{\nabla\phi}{\left|\nabla\phi\right|}\right)\cdot n& \!\!\!\! =0\text{ on }\partial D. \end{cases}\end{split}\]

With Nitsche method, this yields to the following weak formulation: Find \(\phi^{n+1} \in V\), such that for all \(v \in V\):

\[a\left(\phi^{n+1},v\right)=l\left(v,\phi^{n}\right),\]

where :

\[\begin{split}\begin{align} a\left(\phi^{n+1},v\right)&=\int_{D}\nabla\phi^{n+1}\cdot\nabla v\text{ }dx -\int_{\Gamma}\nabla\phi^{n+1}\cdot n_{\Gamma}v\text{ }ds-\int_{\Gamma}\nabla v\cdot n_{\Gamma}\phi^{n+1}\text{ }ds \notag \\ &+\gamma_{D}\int_{\Gamma} \phi^{n+1}v\text{ }ds\label{corrector_bilin}\\ l\left(v,\phi^{n}\right)&=\int_{D}\frac{\nabla\phi^{n}}{\max\left(\left|\nabla\phi^{n}\right|,\epsilon\right)}\cdot\nabla v\text{ }dx \quad\text{ with } \epsilon >0 \text{, very small} \label{corrector_lin} \end{align}\end{split}\]

with \(\phi^0=\phi_{p}\) given by the Predictor problem and \(\gamma_{D}>0\) is the Nistche parameter.

Bilinear form (Corrector):

import ufl

u_r = ufl.TrialFunction(self.V_ls) # Trial function
v_r = ufl.TestFunction(self.V_ls) # Test function

self.gamma_r = 1e+4 # Value of Nitsche parameter

self.h_r = CellDiameter(self.mesh) # mesh size

self.a_correct  = ufl.inner(grad(u_r), grad(v_r))*self.dx
self.a_correct += - dot(grad(u_r), self.n_K)*v_r*self.dsq
self.a_correct += - dot(grad(v_r), self.n_K)*u_r*self.dsq
self.a_correct += self.gamma_r*1.0/self.h_r*u_r*v_r*self.dsq

Linear form (Corrector):

self.L_correct = inner(self.n_K, grad(v_r))*self.dx

The approximation of the normal to isocontour is automatically updated for each iteration.

corrector(level_set)

Returns the solution of the correction problem, denoted \(\phi_{i}\).

Parameters:

level_set (fem.Expression) – The level_set function \(\phi\).

Returns:

The solution to correction problem.

Return type:

fem.Expression

predictor(level_set)

Returns the solution of the prediction problem, denoted \(\phi_{p}\).

Parameters:

level_set (fem.Expression) – The level_set function \(\phi\).

Returns:

The solution to prediction problem.

Return type:

fem.Expression

reinitializationPC(level_set, step_reinit)

Returns the solution of the PC reinitialization method for step_reinit iteration of the correction problem.

Parameters:

• level_set (fem.Expression) – The level_set function \(\phi\).

• step_reinit (int) – The number of iteration for correction problem.

Returns:

The solution to P.C. reinitialization problem.

Return type:

fem.Expression