Compile time array indexing using expression templates -- constexpr?

I want to generate a function like this:

double apply_stencil(const double *u, const int i, const width,
                     const int *offset, const double *weight)
{
  double t=0;
  for(int j=0; j<width; j++)
    t += u[i+offset[j]] * weight[j];
  return t;
}

But I want to make sure that the width, the offsets, and possibly even the weights are compile/time constants.

That is possible to achieve by defining a type:

template <typename S>
double apply_stencil(const double *u, const int i)
{
  double t=0;
  for(int j=0; j<S::width; j++)
    t += u[i+S::offset[j]] * S::weight[j];
  return t;
}

// then:
struct Stencil {
  const static double weight[];
  const static int offset[];
  const static unsigned int width = 3;
};

const double Stencil::weight[] = {1.0, 2.0, 1.0};
const int Stencil::offset[]    = {-1, 0, 1};

However, this is not very pretty. I want a user to be able to specify the Stencil in their application code, and then call my apply_stencil function from my header file (this is really a simplification of something much more complicated).

Ideally, I would like to have the thing specified using expression templates, like so:

const Symbol u;
const Stencil<3> s (1*u[-1] + 2*u[0] + 1*u[1]);

Which uses this infrastructure:

struct Accessor {
  int offset;
};

struct Symbol
{
  Accessor operator[](int i) const {
    return Accessor{i};
  }
};

struct WeightedTerm
{
  double weight;
  int offset;
  WeightedTerm()
    : weight(1), offset(0) {}

  WeightedTerm(double a, Accessor u)
    : weight(a), offset(u.offset) {}
};

WeightedTerm operator*(double a, Accessor u) {
  return WeightedTerm(a,u);
}

template <int n>
struct Sum
{
  WeightedTerm terms[n];

  Sum(WeightedTerm x, WeightedTerm y) {
    terms[0] = x;
    terms[1] = y;
  }

  Sum(Sum<n-1> x, WeightedTerm y) {
    for(int i = 0; i < n-1; ++i)
      terms[i] = x.terms[i];

    terms[n-1] = y;
  }
};

Sum<2> operator+(WeightedTerm x, WeightedTerm y) {
  return Sum<2>(x,y);
}

template <int n>
Sum<n+1> operator+(Sum<n> x, WeightedTerm y) {
  return Sum<n+1>(x,y);
}

template <int width>
struct Stencil
{
  double weights[width];
  int offsets[width];
  Stencil(const Sum<width> s) {

    for(int j = 0; j < width; ++j) {
      weights[j] = s.terms[j].weight;
      offsets[j] = s.terms[j].offset;
    }
  };
};

This looks nice, but now, the arrays are not necessarily compile time known. If I write it like I do above with literals in the expression, I have verified that the compiler can make the correct optimizations. But I want to make find a way to guarantee that they are always compile time constants.

I presume I could encode the offset as a template parameter in Accessor and WeightedTerm, but I can't see how I could do that and keep the desired expression syntax since operator() takes the offset as a regular argument.

So, the question is if there is a way to achieve this? I have a feeling that constexpr could be of use here which I am a bit unfamiliar with.