How to discretize the trilinear form $m(w;,u,v) := \int_\Omega w^{2p} u \ v$ to get a matrix for FEM methods I need to solve the following elliptic non-linear problem using the finite elements method with polynomials of degree $1$ and $2$:
\begin{cases}
  -\Delta u + (u)^{2p} u = f  \qquad & \text{in} \ \ \Omega \\
   u = 0 \qquad & \text{su} \ \ \partial \Omega
\end{cases}
where $\Omega$ is a "standard domain" enclosed by a polygonal curve $\partial \Omega$.
The variational formulation of the problem is in the form
\begin{cases}
\text{Find } u \in H^1_0(\Omega) \\
\int_\Omega \nabla u \cdot \nabla v  + \int_\Omega u^{2p} u \ v = \int_\Omega f \ v \qquad & \forall v \in H^1_0(\Omega)
\end{cases}
Introducing the two (bi-tri) linear forms
$$
a(u,v) := \int_\Omega \nabla u \cdot \nabla v \\
m(w;,u,v) := \int_\Omega w^{2p} u \ v
$$
the problem can be rewritten as
\begin{cases}
\text{Find } u \in H^1_0(\Omega) \\
a(u,v)  + m(u;u,v) = \langle f, v\rangle \qquad & \forall v \in H^1_0(\Omega)
\end{cases}
Now I have to discretize the problem in the algebric form to get a linear system to solve. Since we have the $u^{2p}$ term we will need an iterative part to find the solution we need.
In conclusion we have to solve
\begin{cases}
\text{Find } u_h^n \in V_h \\
a(u_h^{n+1},v_h)  + m(u_h^n;u_h^{n+1},v_h) = \langle f, v_h\rangle \qquad & \forall v_h \in V_h(\Omega)
\end{cases}
where the spaces $V_h$ are the spaces of polynomials of degree $1$ or $2$.
In the construction of the matrix $A$ relative to the form $a(\cdot, \cdot)$ I had no problem (it also has been studied deeply in my course).
The problems arrived when I tried to build the matrix $M$ for the trilinear for $m$. What I thought was to define
$$M_{ij} := m(w;\varphi_i, \varphi_j) = \sum_{T\in \text{Mesh}} \int_T w^{2p} \varphi_i \varphi_j
$$
and the function $w$ would have been the $u_h^{old}$ I found in the last itervative cicle (starting for example with $u_h^0 = 0$). However how could I compute the integral with the discretized function $w$ of which I know only his values in the degrees of freedom? For example, if I'm using the polynomials of degree $1$, I know the solution $u_h^{old}$ just in the vertices of every triangle $T$..
Is it also true that the matrix $M_{ij}$ should be defined like that to solve the problem in the correct way?
Any observation, hint or enlightment would be very appreciated. Thank you.
 A: The short answer is that the polynomials of degree one 1 on each element of the triangular mesh are uniquely determined by their values in vertices of the triangulation. The same holds for polynomials of degree 2 and values in vertices and midpoints of the edges. So knowing the values is in principle the same as knowing the functions.
The usual basis for the space of $P(1)$ elements is given by functions $\phi_i$ such that
$$
\phi_i(x_j) = \delta_{ij}
$$
where $x_j$ are the vertices of the mesh and $\delta_{ij}$ is Kronecker delta. Your fucntion $u_h^{old}$ is then
$$
u_h^{old}(x) = \sum_j v_j \phi_j(x)
$$
where the sum goes over the vertices and $v_j$ are the values on the vertices. Similar relation holds for the higher order elements. For practical application, you typically want to use some library that does this computations for you.
I can recommned Chatper 3 of the book by Brenner, Scott:
The Mathematical Theory of Finite Element Methods for more (probably excesive) details on this.
