I would like to solve the $p$-Laplace dirichlet problem \begin{align*} -\Delta_p u & = 1 \quad\text{in } \Omega:=[0,1]^2 \\ u & = 0 \quad\text{on }\partial\Omega \end{align*} for some $p>1$ with linear finite elements.
This can be done using a Newton method where we iterate over the solution $u^k$ and assemble the stiffness matrix $K = (k_{ij})$ where $$k_{ij} = \int |\nabla u^k|^{p-2}\nabla\varphi_j^T\nabla\varphi_i dx + (p-2)\int |\nabla u^k|^{p-4}(\nabla (u^k)^T\nabla\varphi_j)(\nabla (u^k)^T\nabla\varphi_i)dx$$ in each iteration.
Now my problem is that I do not know how to compute $\nabla u^k$. In each step I know $u^k$, i.e. a vector with the value of $u^k$ at each node, but I don't know how to compute the gradient from this.
Does anyone know any references for this or has experience in how to program this? Thanks!
If $u^k=\sum_{i=1}^NU^k_i\varphi_i$, where $N$ is the number of degrees of freedom, and $\varphi_i,\,\,\,i=1,\ldots,N$ are the basis functions, then $\nabla u^k=\sum_{i=1}^NU^k_i\nabla\varphi_i$.