Although there are similar questions, I am also struggling with the implementation of the following term in my own code by Finite Element Method, namely,  $\nabla \phi \cdot \nabla \phi$. $\phi$ is a state variable such as scalar potential. I think the descritization of this term is much more difficult than that of advection term in N-S equation.  If possible, please tell me how to descritize and implement the term by FEM.The equation to solve is as follows,  $\frac{\partial \phi}{\partial t} = (\nabla\phi)^2 + \nabla^2 \phi.$
Implementation of nonlinear term in FEM
pde;finite element;nonlinear equations
You should rather think of what you will need in the following step, which is probably the numerical time integration of the semi-discrete equations.If you are going to use a (semi) explicit time stepping scheme, all you need is a function that for a given $\phi_0$ assembles the vector $\langle (\nabla \phi_0)^2, v \rangle$, where $v$ are your test functions. In FEM packages like FEniCS this is straight forward.If you going for an implicit scheme, you will probably apply a Newton or Picard iteration for the given iteration point $\phi_0$, you will need the forms or matrices associated with $\langle \nabla \phi_0 \nabla u, v \rangle$ where $u$ are the trial functions. This is a linear form that is easy to assemble with standard FEM methods.For atmost generality you can tensorize your problem, express $(\nabla \phi)^2$ via $D_\nabla (\phi \otimes \phi)$ and assemble the associated linear form (linear in the Kronecker product). For an example considering the convection term $(u \cdot \nabla) u$ in the Navier-Stokes equations have a look at this python function. 