4.5 Declaring your own stencil operators

You can declare your own stencil operators using the macro BZ_DECLARE_STENCIL_OPERATOR1. For example, here is the declaration of Laplacian2D:

BZ_DECLARE_STENCIL_OPERATOR1(Laplacian2D, A)
    return -4*A(0,0) + A(-1,0) + A(1,0) + A(0,-1) + A(0,1);
BZ_END_STENCIL_OPERATOR

To declare a stencil operator on 3 operands, use the macro BZ_DECLARE_STENCIL_OPERATOR3. Here is the declaration of div:

BZ_DECLARE_STENCIL_OPERATOR3(div,vx,vy,vz)
  return central12(vx,firstDim) + central12(vy,secondDim)
    + central12(vz,thirdDim);
BZ_END_STENCIL_OPERATOR

The macros aren’t magical; they just declare an inline template function with the names and arguments you specify. For example, the declaration of div could also be written

template<class T>                              
inline typename T::T_numtype div(T& vx, T& vy, T& vz)   
{
  return central12(vx,firstDim) + central12(vy,secondDim)
                                + central12(vz,thirdDim);
}

The template parameter T is an iterator type for arrays.

You are encouraged to use the macros when possible, because it is possible the implementation could be changed in the future.

To declare a difference operator, use this syntax:

BZ_DECLARE_DIFF(central12,A) {
  return A.shift(1,dim) - A.shift(-1,dim);
}

The method shift(offset,dim) retrieves the element at offset in dimension dim.

Stencil operator declarations cannot occur inside a function. If declared inside a class, they are scoped by the class.