Several nice replies already stated the Pros of finite element methods being flexible and powerful, here I will give another advantage of FEM, from Sobolev space and differential geometry point of view, is that the possibility of finite element space inheriting the physical continuity condition of the Sobolev spaces where the true solution lies in.
For example, Raviart-Thomas face element for plane elasticity, and mixed method for diffusion ; Nédélec edge element for computational electromagnetics.
Normally the solution of a PDE, which is a differential $k$-form lying in the "energy $L^2$-integrable" space:
$$
H\Lambda^k = \{\omega \in \Lambda^k: \omega \in L^2(\Lambda^k), \mathrm{d} \omega \in L^2(\Lambda^k)\}
$$
where $\mathrm{d}$ is the exterior derivative, and we could build the de Rham cohomology around this space, which means we could construct an exact de Rham sequence like the following in 3D space:
$$
\mathbb{R}^3 \xrightarrow{id}\, H(\mathbf{grad},\Omega) \xrightarrow{\nabla}\,
H(\mathbf{curl},\Omega) \xrightarrow{\nabla \times}\,H(\mathrm{div},\Omega)
\xrightarrow{\nabla \cdot}\, L^2(\Omega)
$$
the range of the operator is the null space of the next operator, and there are many nice properties about this, if we could build a finite element space to inherit this de Rham exact sequence, then the Galerkin method based on this finite element space will be stable and will converge to the real solution. And we could get the stability and approximation property of the interpolation operator simply by the commuting diagram from the de Rham sequence, plus we could build the a posteriori error estimation and adaptive mesh refining procedure based on this sequence.
More about this please see Douglas Arnold's article in Acta Numerica: "
Finite element exterior calculus,homological techniques, and applications "
and a slide briefly introducing the idea