Consider the BVP: find $u = u(x)$, for $x \in (0,1)$ that satisfies \begin{align} u'' + u u' = f, \\ u'(0) = g_n, u(1) = g_d. \end{align} To derive the weak form for this BVP, we multiply the first equation by a suitably smooth test function $\Phi = \Phi(x)$ and integrate both sides. This leads to \begin{align} - \int_\Omega \Phi ' u' + \int_\Omega \Phi u u' = \int_\Omega \Phi f + \Phi(0) g_n - \Phi(1) u'(1). \end{align} From this equation, we choose our test function space $\mathscr{T}$ to be a subspace of the Sobolev space $H^1(\Omega)$ that eliminate that $\Phi(1) u'(1)$ term, i.e., \begin{align} \mathscr{T} = \{w: w \in H^1(\Omega), w(1) = 0 \}. \end{align} Furthermore, we choose our set of trial functions $\mathscr{S}$ to be a subset of $H^1(\Omega)$ that satisfy the essential boundary condition, i.e., \begin{align} \mathscr{S} = \{v:v \in H^1(\Omega), v(1) = g_d \}. \end{align}
We then consider the WP associated with the original BVP: find $u \in \mathscr{S}$ such that \begin{align} -\int_\Omega \Phi'u'+\int_\Omega \Phi u u' = \int_\Omega \Phi f + \Phi(0)g_n, \qquad \forall \Phi \in \mathscr{T}. \end{align}
For Galerkin methods, we fix a positive integer $n$ and determine an n-dimensional subspace of $\mathscr{T}$ by specifying a basis, i.e., \begin{align} \mathscr{T}_n = \text{span}\{\phi_1, ..., \phi_n\} \subset\mathscr{T}. \end{align}
For Bubnov-Galerkin, we determine an n-dimensional subset of $\mathscr{S}$ by taking $\mathscr{S}_n = \mathscr{T}_n \cup \{g(x)\}$, where $g(x) = g_d$ handles the essential boundary. Our ``finite element solution' can then be written $u_h(x) = \sum_{j = 1}^n u_j \phi_j(x) + g(x)$.
Substituting $u_h$ into the weak equation should lead to a nonlinear system of equations ($\forall \Phi \in \mathscr{T}$ becomes $\forall \phi_i \in \mathscr{T}_n$): \begin{align} -\int \phi_i \left( \sum_{j = 1}^n u_j \phi_j + g\right)' + \int \phi_i \left( \sum_{j = 1}^n u_j \phi_j + g \right) \left( \sum_{j = 1}^n u_j \phi_j + g \right)' = \int \phi_i f + \phi_i (0). \end{align} The first term on the left-hand-side can be split up into \begin{align} - \sum_{j = 1}^n \int u_j \phi_i \phi_j' - \int \phi_i g', \end{align} which is just a bilinear form (matrix) and linear form (vector). The entire right hand side consists of just linear forms (vectors).
What do I do with the term in the middle?
If I distribute or ``foil'' the middle term, I get \begin{align} \int \phi_i \left(\sum_{j = 1}^n u_j \phi_j \right)\left(\sum_{j = 1}^n u_j \phi_j' \right) + \int \phi_i \left(\sum_{j = 1}^n u_j \phi_j g' \right) + \int \phi_i \left(g \sum_{j = 1}^n u_j \phi_j' \right) + \int \phi_i g g'. \end{align} Here the second and third terms look like bilinear forms (matrices), the fourth term looks like a linear form (vector), but what about the first term, the "trilinear form"? It seems to me that this first term would end up like $u^T A u$ since the unknown coefficients $u_j$ appear twice, but I am not too sure.
Also, my motivation for this problem is to understand similar looking nonlinearities are handled practically, because they seem to appear a lot (e.g., in Navier-Stokes).