You may not be able to express the weak form in deal.II as a mathematical formula, but you come pretty close. For elasticity, the bilinear form reads
$$
a({\mathbf \varphi}_i, {\mathbf \varphi}_j) =
\left(
\lambda \nabla\cdot {\mathbf \varphi}_i, \nabla\cdot {\mathbf \varphi}_j
\right)_\Omega
+
\left(
\mu \nabla\mathbf \varphi_i, \nabla\mathbf \varphi_j
\right)_\Omega,
+
\left(
\mu \nabla\mathbf \varphi_i, \nabla\mathbf \varphi_j^T
\right)_\Omega,
$$
and the corresponding code you would write is as follows:
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<2,dim> phi_i_grad
= fe_values[displacements].gradient (i,q_point);
const double phi_i_div
= fe_values[displacements].divergence (i,q_point);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const Tensor<2,dim> phi_j_grad
= fe_values[displacements].gradient (j,q_point);
const double phi_j_div
= fe_values[displacements].divergence (j,q_point);
cell_matrix(i,j)
+= (lambda_values[q_point] *
phi_i_div * phi_j_div
+
mu_values[q_point] *
double_contract(phi_i_grad, phi_j_grad)
+
mu_values[q_point] *
double_contract(phi_i_grad, transpose(phi_j_grad))
) *
fe_values.JxW(q_point);
}
The last couple of lines assembling the cell matrix are, in effect, a 1:1 transcription of the formula above.
More details here: https://dealii.org/developer/doxygen/deal.II/group__vector__valued.html