Reinit class

The Reinit class solve the Eikonal equation

\[ \begin{aligned} \partial_t \varphi - h^2 \Delta \varphi &= S(\phi)(1 - |\nabla \varphi|) && x\in \mathcal{D},\, t>0, \\ \partial_n \varphi &= 0 && x\in \partial\mathcal{D},\, t>0, \\ \varphi(0, x) &= \phi(x) && x\in \overline{\mathcal{D}}, \end{aligned} \]

where \(\phi\) is the current level set function, \(h\) is the mesh diameter, and

\[ S(\phi) = \frac{\phi}{\sqrt{\phi^2 + h^2}}. \]

---

This class implements the Petrov Galerkin and two-step Adams-Bashforth methods to solve the diffusive Eikonal equation with fictitious time derivative.

Attributes:

Name	Type	Description
dt	Constant	Time step.
phi_ini	Function	Store the initial level set function to be reinitialized.
phi_prev	Function	Store the previous iteration.
w	Function	Auxiliary function.
uh	Function	Solution function.
problem0	LinearProblem	Solver used in the first iteration.
problem	LinearProblem	Solver used in subsequent iterations.

Methods:

Name	Description
run	Run the iterative method over a fixed number of time steps.

Source code in code/formopt.py

class Reinit:
    """
    This class implements the Petrov Galerkin and
    two-step Adams-Bashforth methods to solve
    the diffusive Eikonal equation with
    fictitious time derivative.

    Attributes
    ----------
    dt : Constant
        Time step.
    phi_ini : Function
        Store the initial level set function to be reinitialized.
    phi_prev : Function
        Store the previous iteration.
    w : Function
        Auxiliary function.
    uh : Function
        Solution function.
    problem0 : LinearProblem
        Solver used in the first iteration.
    problem : LinearProblem
        Solver used in subsequent iterations.

    Methods
    -------
    run(phi: Function, steps: int, tend: float) -> None
        Run the iterative method over a fixed number of time steps.
    """

    def __init__(
        self, domain: Mesh, space: FunctionSpace, phi: Function, diam2: float
    ) -> None:
        """
        Set up the reinitialization method.
        Since the bilinear parts change across iterations
        (due to the Petrov-Galerkin test function),
        two solvers are created, which compile
        the equations at each call.

        Parameters
        ----------
        domain : Mesh
            Problem domain.
        space : FunctionSpace
            Space of functions.
        phi : Function
            Level set function.
        diam2 : float
            Square of the mesh diameter.
        """

        self.dt = const(domain, 0.0)
        dx = Measure("dx", domain=domain)

        self.phi_ini = Function(space)
        self.phi_prev = Function(space)

        self.uh = Function(space)
        self.w = Function(space)

        # Sign function
        sign_phi_ini = self.phi_ini / sqrt(self.phi_ini**2 + diam2)
        # Hamiltonian
        H = lambda p: sign_phi_ini * sqrt(dot(p, p))
        # Gradient of H
        GradH = lambda p: sign_phi_ini * p / sqrt(dot(p, p))

        u = TrialFunction(space)
        v = TestFunction(space)

        # Petrov-Galerkin test function
        tau = 2.0 * sqrt(1.0 / self.dt**2 + sign_phi_ini**2 / diam2)
        new_v = v + dot(grad(v), GradH(grad(phi))) / tau

        # Adams-Bashforth weak formulation
        a = u * new_v * dx
        L = phi * new_v * dx
        L += self.dt * sign_phi_ini * new_v * dx
        L += (
            (self.dt / 2.0) * (H(grad(self.phi_prev)) - 3.0 * H(grad(phi))) * new_v * dx
        )
        # Add diffusion
        L += (
            (self.dt / 2.0) * diam2 * dot(grad(self.phi_prev - 3.0 * phi), grad(v)) * dx
        )

        self.problem = create_solver(form(a), form(L), [], self.uh)

        # Explicit Euler weak formulation
        a0 = u * new_v * dx
        L0 = self.dt * sign_phi_ini * new_v * dx
        L0 += (self.phi_ini - self.dt * H(grad(self.phi_ini))) * new_v * dx
        # Add diffusion
        L0 -= self.dt * diam2 * dot(grad(self.phi_ini), grad(v)) * dx

        self.problem0 = create_solver(form(a0), form(L0), [], self.uh)

    def run(self, phi: Function, steps: int, tend: float) -> None:
        """
        Run the iterative method over a fixed number of time steps.
        The solver for the first iteration is called only once,
        and then the main solver is called at each time step.

        Parameters
        ----------
        phi : Function
            Level set function.
        steps : int
            Number of time steps.
        tend : float
            Final time.
        """

        self.dt.value = tend / steps
        self.phi_ini.interpolate(phi)

        self.problem0.solve()
        self.phi_prev.x.array[:] = self.uh.x.array[:]

        for _ in range(steps):
            self.w.x.array[:] = phi.x.array[:]
            self.problem.solve()
            phi.x.array[:] = self.uh.x.array[:]
            self.phi_prev.x.array[:] = self.w.x.array[:]

__init__(domain, space, phi, diam2)

Set up the reinitialization method. Since the bilinear parts change across iterations (due to the Petrov-Galerkin test function), two solvers are created, which compile the equations at each call.

Parameters:

Name	Type	Description	Default
domain	Mesh	Problem domain.	required
space	FunctionSpace	Space of functions.	required
phi	Function	Level set function.	required
diam2	float	Square of the mesh diameter.	required

Source code in code/formopt.py

def __init__(
    self, domain: Mesh, space: FunctionSpace, phi: Function, diam2: float
) -> None:
    """
    Set up the reinitialization method.
    Since the bilinear parts change across iterations
    (due to the Petrov-Galerkin test function),
    two solvers are created, which compile
    the equations at each call.

    Parameters
    ----------
    domain : Mesh
        Problem domain.
    space : FunctionSpace
        Space of functions.
    phi : Function
        Level set function.
    diam2 : float
        Square of the mesh diameter.
    """

    self.dt = const(domain, 0.0)
    dx = Measure("dx", domain=domain)

    self.phi_ini = Function(space)
    self.phi_prev = Function(space)

    self.uh = Function(space)
    self.w = Function(space)

    # Sign function
    sign_phi_ini = self.phi_ini / sqrt(self.phi_ini**2 + diam2)
    # Hamiltonian
    H = lambda p: sign_phi_ini * sqrt(dot(p, p))
    # Gradient of H
    GradH = lambda p: sign_phi_ini * p / sqrt(dot(p, p))

    u = TrialFunction(space)
    v = TestFunction(space)

    # Petrov-Galerkin test function
    tau = 2.0 * sqrt(1.0 / self.dt**2 + sign_phi_ini**2 / diam2)
    new_v = v + dot(grad(v), GradH(grad(phi))) / tau

    # Adams-Bashforth weak formulation
    a = u * new_v * dx
    L = phi * new_v * dx
    L += self.dt * sign_phi_ini * new_v * dx
    L += (
        (self.dt / 2.0) * (H(grad(self.phi_prev)) - 3.0 * H(grad(phi))) * new_v * dx
    )
    # Add diffusion
    L += (
        (self.dt / 2.0) * diam2 * dot(grad(self.phi_prev - 3.0 * phi), grad(v)) * dx
    )

    self.problem = create_solver(form(a), form(L), [], self.uh)

    # Explicit Euler weak formulation
    a0 = u * new_v * dx
    L0 = self.dt * sign_phi_ini * new_v * dx
    L0 += (self.phi_ini - self.dt * H(grad(self.phi_ini))) * new_v * dx
    # Add diffusion
    L0 -= self.dt * diam2 * dot(grad(self.phi_ini), grad(v)) * dx

    self.problem0 = create_solver(form(a0), form(L0), [], self.uh)

run(phi, steps, tend)

Run the iterative method over a fixed number of time steps. The solver for the first iteration is called only once, and then the main solver is called at each time step.

Parameters:

Name	Type	Description	Default
phi	Function	Level set function.	required
steps	int	Number of time steps.	required
tend	float	Final time.	required

Source code in code/formopt.py

def run(self, phi: Function, steps: int, tend: float) -> None:
    """
    Run the iterative method over a fixed number of time steps.
    The solver for the first iteration is called only once,
    and then the main solver is called at each time step.

    Parameters
    ----------
    phi : Function
        Level set function.
    steps : int
        Number of time steps.
    tend : float
        Final time.
    """

    self.dt.value = tend / steps
    self.phi_ini.interpolate(phi)

    self.problem0.solve()
    self.phi_prev.x.array[:] = self.uh.x.array[:]

    for _ in range(steps):
        self.w.x.array[:] = phi.x.array[:]
        self.problem.solve()
        phi.x.array[:] = self.uh.x.array[:]
        self.phi_prev.x.array[:] = self.w.x.array[:]