[feature] Support for tensor fields.

[CD3]

Hello, and thank you for the project.

In physics, we often work with fields, i.e. a tensor that is defined at all points in space. In this case, the components of the tensor are a function of the coordinates. I think it would be useful to have such a concept in the library, because you could then define a metric tensor to compute scalar products in non-euclidean space.

I have a small class that implements this on top of Fastor's Tensor class

template <typename Sig, size_t... Rest>
class TensorField
{
 private:
  template <class S>
  struct get_value_type {
  };
  template <class R, class... Args>
  struct get_value_type<R(Args...)> {
    using type = R;
  };

 public:
  using value_type = typename get_value_type<Sig>::type;
  using function_type = std::function<Sig>;
  using element_type = std::variant<std::monostate, value_type, function_type>;

  TensorField() = default;

  template <typename... Index>
  element_type& get(Index... idx)
  {
    return _storage(idx...);
  }

  template <typename... Args>
  Fastor::Tensor<value_type, Rest...> operator()(Args... args) const
  {
    Fastor::Tensor<value_type, Rest...> l_T;
    for (int i = 0; i < Fastor::pack_prod<Rest...>::value; ++i) {
      l_T.data()[i] = std::visit(
          [&args...](auto&& arg) -> value_type {
            using TT = std::decay_t<decltype(arg)>;
            if constexpr (std::is_same_v<TT, std::monostate>) {
              return 0;
            } else if constexpr (std::is_same_v<TT, value_type>) {
              return arg;
            } else if constexpr (std::is_same_v<TT, function_type>) {
              return arg(std::forward<Args>(args)...);
            }
          },
          this->_storage.data()[i]);
    }

    return l_T;
  }

 private:
  Fastor::Tensor<element_type, Rest...> _storage;
};

Here is an example of using the class to integrate the path length along a curve in cylindrical coordinates.

  // compute the length of a path along constant r,z in cylindrical coordinates.
  // i.e., the circumference of a circle
// define the metric tensor for cylindrical coordiantes
// 1   0   0
// 0  r^2 0
// 0   0   1
  TensorField<double(double, double, double), 3, 3> cylindrical_metric;
  cylindrical_metric.get(0, 0) = 1.;
  cylindrical_metric.get(1, 1) = [](double r, double t, double z) {
    return r * r;
  };
  cylindrical_metric.get(2, 2) = 1.;

// define a parametric path in cylindrical coordinates
// r = r(t) = 3
// theta = theta(t) = 2\pi t
// z = z(t) = 0
// t = [0,1]
  TensorField<double(double), 3> path;
  path.get(0) = 3.;
  path.get(1) = [](double t) { return 2 * M_PI * t; };

  size_t N = 100;
  double dt = 1. / (N - 1);
  double sum = 0;
  enum { I, J };
  for (size_t i = 0; i < N; ++i) {
    auto r = path(dt * i);
    auto dr = path(dt * (i + 1)) - r;
    auto g = cylindrical_metric(r(0), r(1), r(2));
    sum += sqrt(
        Fastor::einsum<Fastor::Index<I>, Fastor::Index<I, J>, Fastor::Index<J>>(
            dr, g, dr)
            .toscalar());
  }
  std::cout << "Result: " << sum << std::endl;
  std::cout << " Exact: " << 2 * M_PI * 3 << std::endl;
  std::cout << "% Diff: " << 100 * (sum - 2 * M_PI * 3) / (2 * M_PI * 3) << "%"
            << std::endl;

which prints

 Result: 19.04
  Exact: 18.8496
 % Diff: 1.0101%

Is this something you would be willing to add to the library? If so, I would be happy to prepare a pull request.

[romeric]

This has been requested before. If you have time, please feel free to prepare a PR. I would be happy to merge it under experimental directory