Why does FEPointEvaluation::EvaluatorTypeTraits use Number and not Tensor<0, dim, Number>?

[nfehn]

When using VectorTools::point_values(), it would be nice to write application code with the return type std::vector<Tensor<rank, dim, Number>>, where rank is a parameter of the application code. However, this does not seem to be possible. I think the reason behind is line 131 in https://www.dealii.org/developer/doxygen/deal.II/fe__point__evaluation_8h_source.html, which reads

using value_type = Number;

and not

using value_type = Tensor<0, dim, Number>;

Wouldn't it be more consistent to write everything with Tensor?

[drwells]

I've run into a similar problem and I had to write my own conversion stubs:

    std::vector<Tensor<1, 1>>
    convert(std::vector<double> &input)
    {
      std::vector<Tensor<1, 1>> result(input.size());
      for (std::size_t i = 0; i < input.size(); ++i)
        result[i][0] = input[i];

      return result;
    }

    template <int dim>
    std::vector<Tensor<1, dim>>
    convert(std::vector<Tensor<1, dim>> &input)
    {
      return std::move(input);
    }
  } // namespace

  // ...
  {
    using VectorType = LinearAlgebra::distributed::Vector<double>;
    auto result =
      VectorTools::point_values<n_components, dim, spacedim, VectorType>(
        *mapping,
        *dof_handler,
        vector,
        evaluation_points,
        remote_point_evaluation);

    return convert(result);
  }

which is not ideal. If we didn't have the specialization, though, it would be Tensor<1, 1, Number>, not the zero tensor (which only exists for defining things recursively). If I recall correctly, we added this overload to make things match FEValues and most other things in the library, where scalars are just Number and not 'tensors with one component'.

That being said - we should try to work out if there is a way to make this generic or at least publish a workaround in the documentation.

[peterrum]

I am not the main author of FPE but I would guess the reason for the special treatment is to be consistent with FEEvalution where we have special implementations for the 1D and the scalar case (in both cases you normally work with doubles instead of tensors). This is a bit inconsistent with the rest of the library but was done - as far as I am aware - for performance reasons.

[kronbichler]

Indeed, the original reason is that this infrastructure imitates what gets done in FEEvaluation in

dealii/include/deal.II/matrix_free/fe_evaluation.h

Line 771 in 5d4eff2

using value_type = Tensor<1, n_components_, VectorizedArrayType>;

and

dealii/include/deal.II/matrix_free/fe_evaluation.h

Line 843 in 5d4eff2

using value_type = VectorizedArrayType;

and

dealii/include/deal.II/matrix_free/fe_evaluation.h

Line 1001 in 5d4eff2

using value_type = Tensor<1, dim, VectorizedArrayType>;

This in turn imitates the types defined by FEValues for the scalar case

dealii/include/deal.II/fe/fe_values.h

Line 154 in 5d4eff2

using value_type = double;

and the vector-valued case

dealii/include/deal.II/fe/fe_values.h

Line 682 in 5d4eff2

using value_type = dealii::Tensor<1, spacedim>;

The reason behind these different choices is that it is more natural to work with scalar types (e.g. a double variable) than with tensors, especially when doing somewhat more complicated mathematical formulas. Performance is not the biggest concern, as most compilers optimize Tensor to essentially same degree as double. (FEEvaluation has some special implementations in a few places to circumvent some unnecessary tensor initializations and copies, but I believe that we made many relevant Tensor operations always_inline might obviate most of them.) Given that we keep expanding the use of these data structures, one could make a point for always just working with Tensor<1, n_components, Number>, which covers all cases with a single interface. Up to now, we've worked hard to equip classes like FEEvaluation with ways to convert between scalar and Tensor<1, 1, Number> to make sure most of the intuitive use cases are covered. As suggsted by @drwells, we might think of ways to convert between the types to ensure we are as flexible as possible.

[nfehn]

If one does not want to switch to Tensor inside the calculations for performance reasons, would it be an option to partially specialize VectorTools::point_values() for the case of Tensor<0,dim> or Tensor<1,1> and then do the conversion described by @drwells internally in VectorTools::point_values(), so that application code can be written more generically? Put differently, one would not change the low-level implementation, but instead provide/apply the utility function for conversion, so that not every user has to write the same code again.