I am trying to implement a Finite Element Solution for the Stabilized Eikonal Equation :
$$ |\nabla u| = 1 + \Gamma \Delta u, \quad \text{ where } \quad u = \text{ distance function }$$
$$ \text{such that } \quad u_{\partial \Omega} = 0 \quad \text{ and }\quad \frac{\partial u}{\partial n} \bigg|_{\partial \Omega} = 0$$ and $$ \Gamma = \text{Some Parameter Depending on the mesh or problem definition}$$
With, $ v = \text{Test Function and } u = \text{Trial Function}$, the weak-form of the equation can be written as:
$\int_\Omega \left( v|\nabla u| \right)dx = \int_\Omega \left[ vdx + v(\Gamma \Delta u )dx \right] $
Integrating by parts and applying boundary conditions:
$$\int_\Omega \left( v|\nabla u| \right)dx = \int_\Omega \left[ vdx - \Gamma (\nabla u \cdot \nabla v)dx \right] $$
$$ F(u; v) = \int_\Omega \left( v|\nabla u| \right)dx - vdx + \Gamma (\nabla u \cdot \nabla v)dx $$
$F(u; v)$ is a Nonlinear function of u. And as a result, I am trying to employ the Newton's(Newton-Rhapson method) for the solving the above non-linear problem.
Considering the example as shown in the section Solving the Nonlinear Variational Problem Directly, the Gateaux derivative $DF(u^k; \delta u, v)$ of $F(u;v)$ at $ u = u\_ $ in the direction of $\delta u$ can be written as:
$$ \begin{align*} J = \lim_{\epsilon\rightarrow 0}{d\over d\epsilon} F_i(u_{-} + \epsilon\delta u; v) & = \lim_{\epsilon\rightarrow 0}{d\over d\epsilon} \left \{ \int_\Omega \left( v|\nabla (u_{-} + \epsilon \delta u)| \right)dx - vdx + \Gamma \nabla (u_{-} + \epsilon \delta u) \cdot \nabla v)dx \right \} \\ & = \lim_{\epsilon\rightarrow 0} \int_\Omega \frac{v \nabla (u_{-} + \epsilon \delta u)}{|\nabla (u_{-} + \epsilon \delta u)|} \cdot \nabla \delta u dx - 0 + \Gamma \nabla \delta u \cdot \nabla v dx \\ & = \int_\Omega \frac{v \nabla \delta u \cdot \nabla u_{-} }{|\nabla u_{-}|} + \Gamma \nabla \delta u \cdot \nabla v \end{align*} $$
Hence, forming $$ [ J ] \delta u = -F(u_{-}; v) \quad \text{where, } u_{-} = \text{Initial Solution for iteration} \\ \text{Next iteration Value is given by } = u_{-} + \omega \delta u \quad \forall \omega \in (0, 1] $$
Clearly, the Jacobian from the above format consist of Nonlinear terms. I have tried multiple ways to compute the expression and have no idea if I am right !!!! .
Can anyone please give me a hint, as, how to calculate the expression for the jacobian of $F(u; v)$ ??