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);
            +=  (phi_i_div * phi_j_div *
                 2 *
                 (phi_i_symmgrad * phi_j_symmgrad) *
                 mu_values[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?

  • 1
    $\begingroup$ The sum symbol is a typo on that page. I submitted a pull request -- thanks for pointing it out! $\endgroup$ – Wolfgang Bangerth Jun 14 at 16:12
  • $\begingroup$ @WolfgangBangerth Thanks for the check. Is my interpretation of the code (in the EDIT) right? $\endgroup$ – bobinthebox Jun 15 at 6:35
  • $\begingroup$ Yes, precisely. phi_i_symmgrad * phi_j_symmgrad has the necessary operator overload that does the double-contraction of symmetric tensors: dealii.org/developer/doxygen/deal.II/… $\endgroup$ – Wolfgang Bangerth Jun 15 at 18:16

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.

  • $\begingroup$ I updated my question. You're totally right about this. Do you agree with my edit? @Tyberius $\endgroup$ – bobinthebox Jun 13 at 22:01
  • $\begingroup$ I don't agree with your last edit, as phi_i_symm is not a scalar, but a rank-2 tensor, as you can see from its definition right after the outermost loop. Right? $\endgroup$ – bobinthebox Jun 14 at 8:00
  • $\begingroup$ @bobinthebox I guess I didn't look carefully enough at the code the first time. I edited my answer again, mostly just to say I agree with your assessment. I was going to check if there was additional info in your link to support this interpretation, but I'm getting 404 errors now whenever I visit the dealii site. $\endgroup$ – Tyberius Jun 14 at 14:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.