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.
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.