FEM approach to obtain the projection of a function onto the mesh?

I need to explain the background of the problem I'm working on, and the questions themselves are a couple of paragraphs down. I am reading this book and implementing the algorithms along (so far in chapter one, assemble the mass matrix and the load vector in 1D). As I go on, I'm writing unit tests. The one for the assembly of the load vector is failing, either because of a programming bug, or because I'm failing to understand the math. Please let me know if it's the latter.

At his point the book isn't laying out any differential equation to solve. It is just defining the mass matrix and load vector, on the mesh/interval $$I$$, in the context of approximation and projection. In this context the mass matrix is defined by the book as $$M_{ij} = \int_I \phi_j \phi_i \,dx$$ and the load vector for a function $$f(x)$$ as $$b_i = \int_I f \phi_i \,dx$$ where $$\phi_i$$ are the "hat functions," basis of the subspace of piecewise linear functions on the mesh of $$I$$.

So the whole point of these definitions in this chapter of the book (again, not yet laying out any differential equation to solve) was supposedly that the equation $$\mathbf{M \xi = b}$$ holds for the coordinates $$\xi_i$$ of the $$L^2$$-projection of $$f(x)$$ onto the mesh.

QUESTION 1: is this so? Have I understood the book correctly so far?

On this assumption I based my unit tests. First, I understand that $$\sum_i \xi_i \phi_i(x) \approx f(x)$$ Moreover, I understand that if I choose a mesh $$I$$ and a function $$f$$ which happens to be already piecewise linear on $$I$$, the approximation should become identity.

Therefore I made up the following example: $$\mathbf{x} = \{ -10, -4, 8 \} \\ f(x) = \left| 2 x + 8 \right|$$ which gives me the following mass matrix and load vector: $$\mathbf{M} = \begin{bmatrix} 2 & 1 & 0 \\ 1 & 6 & 2 \\ 0 & 2 & 4 \end{bmatrix} \\ \mathbf{b} = \{ 36, 0, 144 \}$$ Since the projection should equal the original function, because the latter was already in the subspace (*): $$\mathbf{M \cdot \xi \stackrel{*}{=} M} \cdot \{ f\left(\mathbf{x}(1)\right), f\left(\mathbf{x}(2)\right), f\left(\mathbf{x}(3)\right) \}$$ However this gives $$\mathbf{M} \cdot \{ 12, 0, 24 \} = \{ 24, 60, 96 \} \neq \mathbf{b}$$ Not even approximate.

QUESTION 2: what have I done or understood wrong?

The calculation of the load vector is done numerically. The book uses the Trapezoidal rule, and proposes the following Matlab code:

function b = LoadAssembler1D(x,f)
n = length(x)-1;
b = zeros(n+1,1);
for i = 1:n
h = x(i+1) - x(i);
b(i) = b(i) + f(x(i))*h/2;
b(i+1) = b(i+1) + f(x(i+1))*h/2;
end


QUESTION 3: is this implementation correct, assuming the Trapezoidal approximate integration?

This code computes for this mesh and $$f(x)$$ the following load vector:

>> LoadAssembler1D( [ -10, -4, 8 ], inline('abs(2 * x + 8)', 'x') )
ans =
36
0
144


where I got the discrepancy above. What's wrong? Is it the implementation, or something else?

Is my mesh perhaps so coarse that the numerical integration may be to blame for this whole error, even if the code is correct based on it? I had the idea that the Trapezoidal rule would be exact for this simple case, but now I realize that my load vector is made up of integrals of piecewise quadratic, not linear, functions...

Or is there something else wrong that I've missed?

• The mesh is too coarse. I think you are just getting a huge integration error while integrating the second degree polynomials ($f \phi_i$) with a first degree quadrature. Normally, most FE solvers (especially for 2D, 3D problems) will assemble the matrices pulling back to a reference element, where Gauss quadrature with higher precision is used. – PierreCarre Feb 6 at 15:30
• For instance you are considering $\int_{-4}^8 \frac{1}{12}(2x+8)(x+4) \,dx \approx \frac{12}{2}(0+24) = 144$, which is not even close to the real value. – PierreCarre Feb 6 at 15:41
• Indeed; updated the answer... thank you! – J.P. Feb 6 at 19:38

You have understood correctly. The $$L^2$$ projection, $$f_h$$, is the solution of the problem $$\langle f_h - f, v_h \rangle =0, \forall v_h \in V_h$$. Writing $$f_h$$ in terms of the basis of $$V_h$$, $$f_h = \sum_i \xi_i \phi_i$$ and choosing $$v_h= \phi_1, \cdots \phi_n$$ you get the linear system you mentioned. So yes, the solution of $$M \xi = b$$ should give you the $$L^2$$ projection.
$$b_i = \int_{-10}^8 f(x) \phi_i(x)\,dx$$
yields $$b = (24, 60, 96)$$ (see the note below), and the solution of the system $$M \xi = b$$ is $$\xi = (12,0,24)$$, as expected.
$$\textbf{Note}$$: for instance $$b_3 = \int_{-4}^8 \frac{1}{12}(2x+8)(x+4) dx = 96.$$