the Lagrangian density is $L=T-U$, the difference between the kinetic and potential energy density, as we're used to in all of mechanics.
The kinetic energy density is $T = \rho v(x,y,z)^2/2$ where one has to calculate the density $\rho$ properly. The potential energy is more general and complicated,
$$U = \frac{1}{2} C_{ijkl}u_{ij}u_{kl}$$
where the tensor $C$ with four indices is called the elastic. The tensors $u$ are obtained from the displacement vectors as
$$ u_{ij} = \frac{1}{2} (\partial_i u_j + \partial_j u_i) $$
Note that the density $\rho$ in the kinetic energy also depends on the displacement vectors, if you want to re-express it via the mass density at rest (without displacement).
For isotropic materials,
$$ C_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk} ) $$
where $\lambda$ and $\mu$ are known as Lamé constants. For more general (but linear) crystals, however, $C$ is a general tensor symmetric under the $ij$ exchange, $kl$ exchange, and the exchange of $ij$ and $kl$ as pairs.
Best Answer
To construct a general spin 2 Lagrangian, you might do the following.
To begin with, we notice that Wigners classification tells you that spin-2 particles must be described by a rank-2 tensor field, i.e. $h_{\mu\nu}$.
Unitarity demands that we have a theory with at most two derivatives, and Lorentz invariance tells us that a theory containing spin two fields must have either zero or two derivatives.
Starting at second order, we would write $$ \mathcal{L} = a_1h^{\mu\nu}\partial_\mu\partial_\nu h + a_2 h^{\mu\nu}\partial_\lambda\partial_\nu h^\lambda_\nu + a_3 h\partial_\nu\partial_\mu h^{\mu\nu} + a_4 h^{\mu\nu}\partial^2 h_{\mu\nu} + a_5 h\partial^2 h + a_6h_{\mu\nu}h^{\mu\nu} + a_7 h^2 $$ We also note that since this Lagrangian is part of an action, there are many more possible terms we could write down, however the ones above form a minimal basis, since all others can be related via integration by parts.
If we now further demanding that our theory is gauge invariant, we would yield the Fierz-Pauli Lagrangian, fixing the coefficients to be $$ \left\{a_1 = 0, a_2 = x,a_3 = -x, a_4 = -\frac{x}{2}, a_5 = \frac{x}{2}, a_6 = a_7 = 0\right\} $$ Where $x$ is a free parameter. However, if gauge invariance isn't a requirement then you can be content with the Lagrangian density above. Depending on what is required of your theory, you can fix some coefficients by demanding Newtonian gravity like limits, or make specific demands of the spin-1 and spin-0 modes of $h_{\mu\nu}$.
Additionally, you can also play the same game at third order if you want to consider self-interacting spin-2 fields, and construct a third order Lagrangian in exactly the same was as before, by considering \begin{align*} \mathcal{L} &= a_1 h^{\alpha \beta } \partial_{\alpha }h_{\delta \lambda } \partial_{\beta }h^{\delta \lambda } + \frac12 a_2 \left(h^{\alpha \beta } \partial_{\alpha }h_{\delta \lambda } \partial^{\lambda }h^{\delta }{}_{\beta } + h^{\alpha \beta } \partial_{\beta }h_{\delta \lambda } \partial^{\lambda }h^{\delta }{}_{\alpha }\right)+a_3 h^{\alpha \beta } \partial_{\lambda }h_{\alpha \delta } \partial^{\lambda }h^{\delta }{}_{\beta }\\ &~~+ a_4 h^{\alpha \beta } \partial_{\lambda }h_{\alpha \delta } \partial^{\delta }h^{\lambda }{}_{\beta } + \frac{1}{2} a_5 \left(h^{\alpha \beta } \partial_{\delta }h^{\delta \lambda } \partial_{\alpha }h_{\beta \lambda }+h^{\alpha \beta } \partial_{\delta }h^{\delta \lambda }\partial_{\beta }h_{\alpha \lambda }\right)+a_6 h^{\alpha \beta } \partial_{\delta }h^{\delta \lambda } \partial_{\lambda }h_{\alpha \beta }\\&~~+\frac{1}{2} a_7 \left(h^{\alpha \beta } \partial^{\lambda }h \partial_{\alpha }h_{\beta \lambda }+h^{\alpha \beta } \partial^{\lambda }h\partial_{\beta }h_{\alpha \lambda }\right) + a_8 h^{\alpha \beta } \partial^{\lambda }h \partial_{\lambda }h_{\alpha \beta } \\&~~+\frac{1}{2} a_9 \left(h^{\alpha \beta } \partial_{\beta }h \partial_{\lambda }h^{\lambda }{}_{\alpha }+h^{\alpha \beta }\partial_{\alpha }h \partial_{\lambda }h^{\lambda }{}_{\beta }\right) + a_{10} h^{\alpha \beta } \partial_{\alpha }h \partial_{\beta }h+a_{11} h^{\alpha \beta } \partial_{\lambda }h^{\lambda }{}_{\alpha } \partial_{\delta }h^{\delta }{}_{\beta }\\ &~~+a_{12} \partial^{\gamma }h^{\delta \lambda } \partial_{\gamma }h_{\delta \lambda }+a_{13} \partial^{\delta }h^{\gamma \lambda } \partial_{\gamma }h_{\delta \lambda }+a_{14} \partial_{\gamma }h^{\gamma }{}_{\delta } \partial_{\lambda }h^{\delta \lambda }+a_{15} \partial_{\delta }h \partial_{\gamma }h^{\gamma \delta }+ a_{16} h\partial_{\gamma }h \partial^{\gamma }h \end{align*}
Playing the same game with Gauge invariance again fixes all your coefficients, up to some redundancies, and you recover the weak-field GR third order Lagrangian, however again you can choose how to fix your coefficients depending on the theory under consideration.