Confusion about bilinear form for elasticity equation in deal.ii tutorial

Question

I'm learning how to solve vector-valued problems with deal.II library. In particular, I'm looking at the following introduction from the official website https://www.dealii.org/current/doxygen/deal.II/group__vector__valued.html#VVAlternative

Here they write the bilinear form as $$a(\mathbf{u}, \mathbf{v})=\left( \lambda \nabla\cdot {\mathbf u}, \nabla\cdot {\mathbf v} \right)_\Omega + 2 \sum_{i,j} \left( \mu \frac 12[\partial_i u_j + \partial_j u_i], \frac 12[\partial_i v_j + \partial_j v_i] \right)_\Omega$$

and then they say $$=\left( \lambda \nabla\cdot {\mathbf u}, \nabla\cdot {\mathbf v} \right)_\Omega + 2 \sum_{i,j} \left( \mu \varepsilon(\mathbf u), \varepsilon(\mathbf v) \right)_\Omega$$

but I can't understand why that double sum over indices $i,j$ in the last equality! I think it should not be there, also because there's no dependence on $i,j$.

---

EDIT: So the last summand is $$\sum_{i,j} (\mu \varepsilon(\mathbf{u})_{ij}, \varepsilon(\mathbf{v})_{ij})_{\Omega}$$

I'm wondering now how this double sum is translated, precisely, in the snippet

for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
  for (unsigned int i=0; i<dofs_per_cell; ++i)
    {
      const SymmetricTensor<2,dim> phi_i_symmgrad
        = fe_values[displacements].symmetric_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 SymmetricTensor<2,dim> phi_j_symmgrad
            = fe_values[displacements].symmetric_gradient (j,q_point);
          const double phi_j_div
            = fe_values[displacements].divergence (j,q_point);
          cell_matrix(i,j)
            +=  (phi_i_div * phi_j_div *
                 lambda_values[q_point]
                 +
                 2 *
                 (phi_i_symmgrad * phi_j_symmgrad) *
                 mu_values[q_point]) *
                fe_values.JxW(q_point);
        }
    }

Here's my attempt:

Applying quadrature formula: $$\sum_q \mu(x_q) \underbrace{\sum_{i,j} \varepsilon(\mathbf{u}(x_q))_{ij} \varepsilon(\mathbf{v}(x_q))_{ij}}_{\star}$$

Now I recognise that the $\star$ term is precisely $$\varepsilon(\mathbf{u}) : \varepsilon(\mathbf{v})$$

Therefore, I suspect that phi_i_symmgrad * phi_j_symmgrad is precisely the scalar product between the two symmetrized gradients, where the * operator has been properly overloaded since both terms are tensors of rank 2.

Is that correct?

Answer 1

Its not the summation that is wrong, but the lack of indices inside it. Below this expression on their site they define: $$\epsilon(\mathbf{u})=\frac{1}{2}([\nabla\mathbf{u}]+[\nabla\mathbf{u}]^{\mathrm{T}})$$

So $\epsilon(\mathbf{u})$ is a matrix formed from the symmetrized gradient of $\mathbf{u}$. But to get the sum from the earlier expression, we want to sum over the elements of this matrix. So it should really be $\epsilon(\mathbf{u})_{ij}$.

For the translation from formula to code, I think you have the right idea. I had initially mistaken the $i$ and $j$ for loops as the summations from the formula, but these seem to rather be contributions from different degrees of freedom of the system. I had missed that phi_i_symm was explicitly defined as a SymmetricTensor and so * must be overridden to the scalar product for this type in order for the code to make sense dimensionally.