libmesh_not_implemented();
  }

  TypeNTensor(const TypeVector<T> &) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N)))
  {
    libmesh_not_implemented();
  }

  TypeNTensor(const TypeTensor<T> &) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N)))
  {
    libmesh_not_implemented();
  }

  TypeNTensor(const TypeNTensor<N, T> & TN) : _coords(TN._coords) {}

  /**
   * Constructor from the flattened coefficients, ordered as \p _coords.

    libmesh_assert_equal_to(vec.size(), std::size_t(int_pow(LIBMESH_DIM, N)));
  }

  TypeNTensor & operator=(const TypeNTensor<N, T> & TN)
  {
    _coords = TN._coords;
    return *this;
  }

  operator TypeVector<T>() const

  /**
   * \returns A writable proxy for the \f$ i^{th} \f$ slice of the tensor.
   */
  TypeNTensor<N - 1, T> slice(const unsigned int i)
  {
    libmesh_assert_less(i, LIBMESH_DIM);
    TypeNTensor<N - 1, T> slice;
    unsigned int slice_size = int_pow(LIBMESH_DIM, N - 1);
    slice._coords.resize(slice_size);
    for (unsigned int j = 0; j < slice_size; j++)
      slice._coords[j] = _coords[i * slice_size + j];
    return slice;
  }

  /**

   * Add to this tensor.
   */
  template <typename T2>
  const TypeNTensor<N, T> & operator+=(const TypeNTensor<N, T2> & other)
  {
    unsigned int size = int_pow(LIBMESH_DIM, N);
    for (unsigned int i = 0; i < size; i++)
      _coords[i] += other._coords[i];
    return *this;
  }

   */
  template <typename T2>
  TypeNTensor<N, typename CompareTypes<T, T2>::supertype>
  operator-(const TypeNTensor<N, T2> & other) const
  {
    TypeNTensor<N, typename CompareTypes<T, T2>::supertype> subtract;
    unsigned int size = int_pow(LIBMESH_DIM, N);
    subtract._coords.resize(size);
    for (unsigned int i = 0; i < size; i++)
      subtract._coords[i] = _coords[i] - other._coords[i];
    return subtract;
  }

  /**

  typename std::enable_if<
    ScalarTraits<Scalar>::value,
    TypeNTensor<N,typename CompareTypes<T, Scalar>::supertype>>::type
  operator * (const Scalar factor) const
  {
    TypeNTensor<N, typename CompareTypes<T, Scalar>::supertype> multiplied;
    unsigned int size = int_pow(LIBMESH_DIM, N);
    multiplied._coords.resize(size);
    for (unsigned int i = 0; i < size; i++)
      multiplied._coords[i] = factor * _coords[i];
    return multiplied;
  }

  /**

   * \returns The Frobenius norm of the tensor, i.e. the square-root of
   * the sum of the elements squared.
   */
  auto norm() const -> decltype(std::norm(T())) { return std::sqrt(norm_sq()); }

  /**
   * \returns The Frobenius norm of the tensor squared, i.e. the sum of the

   */
  auto norm_sq() const
  {
    unsigned int size = int_pow(LIBMESH_DIM, N);
    Real sum = 0.;
    for (unsigned int i = 0; i < size; i++)
      sum += TensorTools::norm_sq(_coords[i]);
    return sum;
  }

  /**
   * Set all entries of the tensor to 0.
   */
  void zero()
  {
    unsigned int size = int_pow(LIBMESH_DIM, N);
    for (unsigned int i = 0; i < size; i++)
      _coords[i] = T(0);
  }

  /**
   * \returns \p true if two tensors are equal, \p false otherwise.