When people say use Legendre polynomial as basis polynomial in fem, are they using the polynomials themselves or the integral of those polynomials?
I'm asking this because I have this poisson equation: $\nabla^2u=f$, with $u(-1)=u(1)=0$. I think for the boundary condition to hold, we need the basis polynomials to vanish on the boundary. But surely that's not true for Legendre polynomials? That's why I think they are probably using the integrals, which do vanish on the boundary.
Does the above sound right? Or do people actually use Legendre polynomials? If so, what does the problem look like?
The Legendre Polynomials are indeed difficult to use, especially in Continuous Galerkin schemes, where you would not only need to satisfy the boundary conditions, but also inter-element smoothness. You could use them without a mesh and do a spectral element method, but by far the most common use of Legendre polynomials as basis functions will be with Discontinous Galerkin (DG) Schemes.
Even with DG you still need to tackle the boundary conditions. The most common way to do this is to enforce it weakly. Arnold et al (Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems 2002 SIAM J. Numer. Analysis) describe enforcing it through a penalty term for a weak form: $$ \int_\Omega\nabla u\cdot \nabla v\;d\Omega + \int_{\partial U} \mu(u - g)v\;ds = \int_\Omega fv \; d\Omega \quad\forall v\in H^1(\Omega) $$ Here $u=g$ on $\partial\Omega$ is your boundary condition (for your case $g=0$) and $\mu$ is some very large number (if you implement this try $10^8$)
An alternate weak form described right after (this time assuming $g=0$) is:
$$ \int_\Omega\nabla u\cdot \nabla v\;d\Omega - \int_{\partial U} \frac{\partial u}{\partial n}u\;ds - \int_{\partial U} \frac{\partial v}{\partial n}u\;ds + \int_{\partial U} \mu uv\;ds= \int_\Omega fv \; d\Omega \quad\forall v\in H^1(\Omega) $$
The Arnold et al. reference described earlier also has methods for calculating the numerical flux which is very important for DG methods.
The Dirichlet boundary condition can be enforced strongly but that usually involves converting to a nodal basis and doing some extra bookkeeping with your degrees of freedom. As for the reason you would use Legendre polynomials: since the polynomials are orthogonal with respect to the $L^2$ inner product, the mass matrix is diagonal which can improve calculation efficiency.