Optimize product of symmetric tensors

[bangerth]

Internally, SymmetricTensor<2,dim> stores its elements as a Tensor<1,n_independent_components>, and SymmetricTensor<4,dim> stores its elements as a Tensor<2,n_independent_components>, where n_independent_components=dim*(dim-1)/2. In other words, we use something like Voigt notation.

But operator* does not make use of this, but laboriously applies index translations:

template <int rank_, int dim, typename Number>
template <typename OtherNumber>
DEAL_II_HOST DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
  typename internal::SymmetricTensorAccessors::
    double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
    SymmetricTensor<rank_, dim, Number>::operator*(
      const SymmetricTensor<2, dim, OtherNumber> &s) const
{
  // need to have two different function calls
  // because a scalar and rank-2 tensor are not
  // the same data type (see internal function
  // above)
  return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
                                                                        s.data);
}



template <int rank_, int dim, typename Number>
template <typename OtherNumber>
DEAL_II_HOST DEAL_II_CONSTEXPR inline
  typename internal::SymmetricTensorAccessors::
    double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
    SymmetricTensor<rank_, dim, Number>::operator*(
      const SymmetricTensor<4, dim, OtherNumber> &s) const
{
  typename internal::SymmetricTensorAccessors::
    double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type tmp;
  tmp.data =
    internal::perform_double_contraction<dim, Number, OtherNumber>(data,
                                                                   s.data);
  return tmp;
}

We should really make this more efficient -- all it takes is a product of two rank-2 and rank-1 tensors :-)

[bangerth]

On second look, the perform_double_contraction() functions are actually about as good as it gets. What I stated above isn't true: It isn't just a matrix-vector or vector-vector product because one has to take into account that some terms of the sum appear twice, but are of course only implemented once.