template <typename T, int D> struct vect {
  array<T, D> elts;

  vect() {
    for (int d = 0; d < D; ++d)
      elts[d] = 0;
  }
  vect(const array<T, D> &arr) : elts(arr) {}
  static vect unit(int dir) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = d == dir;
    return r;
  }
  const T &operator[](int d) const { return elts[d]; }
  T &operator[](int d) { return elts[d]; }
  friend vect operator+(const vect &x) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = +x.elts[d];
    return r;
  }
  friend vect operator-(const vect &x) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = -x.elts[d];
    return r;
  }
  friend vect operator+(const vect &x, const vect &y) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = x.elts[d] + y.elts[d];
    return r;
  }
  friend vect operator-(const vect &x, const vect &y) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = x.elts[d] - y.elts[d];
    return r;
  }
  friend vect operator*(T a, const vect &x) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = a * x.elts[d];
    return r;
  }
  friend vect operator*(const vect &x, T a) {
    vect r;
    for (int d = 0; d < D; ++d)
      r.elts[d] = x.elts[d] * a;
    return r;
  }
  friend ostream &operator<<(ostream &os, vect &x) {
    os << "[";
    for (int d = 0; d < D; ++d) {
      if (d > 0)
        os << ",";
      os << x.elts[d];
    }
    os << "]";
    return os;
  }
};
////////////////////////////////////////////////////////////////////////////////

  }
  inline T &restrict operator()(const vect<int, dim> &I) const {
    return ptr[offset(I[0], I[1], I[2])];

  vect<int, dim> I;
  vect<int, dim> DI(int d) const { return vect<int, dim>::unit(d); }

          const PointDesc p{i, j, k, x, y, z, idx, dj, dk};

          vect<int, dim> I{{i, j, k}};
          const PointDesc p{i, j, k, x, y, z, idx, dj, dk, I};