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);
          }