update to describe the current state of this thorn

git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@159 df1f8a13-aa1d-4dd4-9681-27ded5b42416

1 files changed, 35 insertions, 1626 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index a2b97d9..c41c596 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -76,7 +76,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % The author of the documentation
-\author{Jonathan Thornburg, incorporating code from Thomas Radke} 
+\author{Thomas Radke} 
 
 % The title of the document (not necessarily the name of the Thorn)
 \title{Thorn Guide for the {\bf LocalInterp} Thorn}
@@ -107,10 +107,12 @@
 
 % Add an abstract for this thorn's documentation
 \begin{abstract}
-This thorn provides processor-local interpolation of 1-D, 2-D, and 3-D
-data arrays.  It supports both the older \verb|CCTK_InterpLocal()| API,
-and the newer and more generic \verb|CCTK_InterpLocalUniform()| API.
-This thorn guide describes how to use the various interpolators.
+This thorn does processor-local interpolation of N-dimensional data
+arrays.  In general there may be many input arrays (all defined on the
+same uniform Cartesian grid) all being interpolated to the same set
+of interpolation points.  At present this thorn only supports the
+\verb|CCTK_InterpLocal()| API, but in the near future it will probably be
+enhanced to support the newer \verb|CCTK_InterpLocalUniform()| API.
 \end{abstract}
 
 % The following sections are suggestive only.
@@ -120,205 +122,18 @@ This thorn guide describes how to use the various interpolators.
 
 \section{Introduction}
 
-At present this thorn provides 2~interpolators (we hope to add other
-interpolators in the future):
-\begin{description}
-\item[Uniform Cartesian]
-	This is the local interpolator which used to live
-	in the \thorn{PUGHInterp} thorn.  It was written in C by
-	Thomas Radke in early 2001, drawing on earlier Fortran code
-	by Paul Walker.  It supports the \verb|CCTK_InterpLocal()| API.
-	It provides Lagrange polynomial interpolation for all
-	combinations of~1, 2, and 3~dimensions, and orders~1,
-	2, and~3.
-\item[Generalized Polynomial]
-	This interpolator was written in C (plus Maple to generate
-	the coefficients) by Jonathan Thornburg in winter 2001-2002.
-	It supports the \verb|CCTK_InterpLocalUniform()| API.
-	It provides two variants of Lagrange polynomial interpolation in
-	1~dimension (orders~1-6) and~2 and 3~dimensions (orders~1-4),
-	and Hermite polynomial interpolation in
-	1~dimension (orders~2-4) and~2 and 3~dimensions (orders~2 and~3).
-	This interpolator supports a number of options, described below.
-\end{description}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Basic Terminology}
-\label{localinterp/sect-basic-terminology}
-
-Within Cactus, each interpolator has a character-string name;
-this is mapped to a Cactus \defn{interpolator handle} by
-\verb|CCTK_InterpHandle()|.  For any given interpolator handle,
-there may be a separate interpolator defined for each of the
-interpolation APIs (both the processor-local ones provided by this
-thorn, and the global ones provided by driver-specific interpolator
-thorns such as \thorn{PUGHInterp}).
-
-Terminology for interpolation seems to differ a bit from one author
-to another.  Here we refer to the \defn{point-centered} interpolation
-of a set of \defn{input arrays} (defining data on a uniformly or
-nonuniformly spaced \defn{grid} of \defn{data points}) to a set of
-\defn{interpolation points} (specified by a corresponding set of
-\defn{interpolation coordinates}), with the results being stored
-in a corresponding set of \defn{output arrays}.  Alternatively,
-we may refer to \defn{cell-centered} interpolation, using a grid
-of \defn{data cells} and a set of \defn{interpolation cells}.
-(This latter terminology is common in hydrodynamics interpolators.)
-
-At present all the interpolators do polynomial interpolation, and
-allow the interpolation of multiple input arrays (to the same set of
-interpolation points) in a single interpolator call, using the basic
-algorithm:
+This thorn provides processor-local interpolation, using the interpolation
+operators
 \begin{verbatim}
-for each interpolation point
-{
-choose an interpolation molecule position
-  somewhere near the interpolation point
-        for each output array
-        {
-        compute an interpolating polynomial which approximates
-          the input data at the molecule points
-        output = polynomial(interpolation point)
-        }
-}
+  "first-order uniform cartesian"
+  "second-order uniform cartesian"
+  "third-order uniform cartesian"
 \end{verbatim}
-In the future, we may add other interpolators where the choice of
-molecule is data-dependent (and thus may vary from one input array to
-the next), and/or where the entire input grid is used in each interpolation.
-
-We define the \defn{order} of the interpolation to be the order of
-the fitted polynomial.  That is, in our terminology linear interpolation
-is order~1, quadratic is order~2, cubic is order~3, etc.  An order~$n$
-interpolator thus has $O\big((\Delta x)^{n+1})\big)$ interpolation errors
-for generic smooth input data.
-Section~\ref{localinterp/sect-multi-dim-interp} explains how the
-interpolating polynomial is defined for multi-dimensional interpolation.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{The Non-Smoothness of Interpolation Errors}
-
-Because the interpolating polynomial generally changes if
-the interpolation point moves from one grid cell to another, unless
-something special is done the interpolating function isn't smooth,
-\ie{} its 1st~derivative is generically {\em discontinuous\/},
-with $O(\Delta x^n)$ jump discontinuities each time the interpolating
-polynomial changes.  Correspondingly, the interpolation error is
-generically a ``bump function'' which is zero at each grid point
-and rises to a local maximum in each grid cell.  This is the case,
-for example, for tensor-product Lagrange polynomial interpolation.
-
-[For maximum-degree Lagrange polynomial interpolation, the
-interpolation error may not be zero at the grid points, and the
-interpolant itself may not be continuous there.  
-Section~\ref{localinterp/sect-multi-dim-interp} explains this
-in detail.]
-
-This thorn also provides Hermite polynomial interpolation, which
-guarantees a smooth ($C^1$) interpolant and (for smooth input data)
-interpolation error.  Unfortunately, this comes at the cost of a
-larger molecule size than the same-order Lagrange interplator, and
-a much larger interpolation error if the interpolation molecules are
-off-centered.  (By default, the Hermite interpolator doesn't allow
-off-centered molecules, though this can be changed via the
-\verb|boundary_off_centering_tolerance[]| and
-\verb|excision_off_centering_tolerance[]| parameter-table entries.)
-
-For most purposes Lagrange polynomial interpolation is better;
-only use Hermite polynomial interpolation if you need a smooth
-interpolant.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{More Terminology}
-
-As given in the Function Reference section of the Cactus User's Guide,
-\verb|interp_coords|, \verb|input_arrays|, and \verb|output_arrays| are
-actually all pointers to arrays of \verb|void *| pointers, since we
-support a number of different Cactus data types.  Internally, the
-interpolator casts these \verb|void *| pointers to \verb|CCTK_REAL *|
-or whatever the correct Cactus data types are.  But when describing
-how the interpolator accesses the various arrays, for simplicity we
-usually gloss over this casting, \ie{} we pretend that \verb|interp_coords|,
-\verb|input_arrays|, and \verb|output_arrays| are pointers to arrays
-of \verb|CCTK_REAL *| pointers.  (This may become clearer once you
-read the next example.)
 
-We use \verb|pt|, \verb|in|, and \verb|out| as generic 0-origin integer
-subscripts into the arrays of interpolation points, input arrays, and
-output arrays respectively.  We use \verb|(i,j,k)| as a generic
-\verb|N_dims|-vector of integer subscripts into the input array
-\verb|input_arrays[in]|.  (Thus \defn{{\tt (i,j,k)} space} refers to
-the grid of data points.)  We usually only write array subscripting
-expressions for the 3-D case; the restrictions/generalizations to
-other dimensions should be obvious.
-
-For example, for 3-D interpolation, the (x,y,z) coordinates of a typical
-interpolation point are given by
-\begin{verbatim}
-x = interp_coords[0][pt]
-y = interp_coords[1][pt]
-z = interp_coords[2][pt]
-\end{verbatim}
-(Notice here that as described above, we've implicitly taken
-\verb|interp_coords| to have the C~type
-\verb|const CCTK_REAL* interp_coords[]|, glossing over the casting
-from its actual C~type of \verb|const void* interp_coords[]|.)
-
-We take \verb|axis| to be a 0-origin integer specifying a coordinate
-axis (dimension), \ie{} 0~for~$x$, 1~for~$y$, 2~for~$z$, \dots.
-We take \verb|ibndry| to be a 0-origin integer specifying the
-combination of a coordinate axis (dimension) and a minimum/maximum
-``end'' of the grid; these are enumerated in the order
-$x_{\min}$, $x_{\max}$, $y_{\min}$, $y_{\max}$, $z_{\min}$,
-$z_{\max}$, \dots.
-
-When describing Jacobians and domains of dependence, it's useful to
-introduce the notion of \defn{molecule position}, a nominal reference
-point for the interpolation molecule in \verb|(i,j,k)| coordinate
-space.  (For example, the molecule position might just be the \verb|(i,j,k)|
-coordinates of the molecule's center.)  We also introduce
-\defn{molecule coordinates} \verb|(mi,mj,jk)|, which are just
-\verb|(i,j,k)| coordinates relative to the molecule position.
-We use \verb|m| or as a generic molecule coordinate.  Thus
-(in notation which should be obvious) a generic molecule operation
-can be written
-\begin{equation}
-\verb|output| = \sum_{\tt m} \verb|coeff[posn+m]| \times \verb|input[posn+m]|
-\end{equation}
-Note that there is no requirement that the output be semantically
-located at the position \verb|posn|!  (This may become clearer once
-you look at the example in
-section~\ref{localinterp/sect-Jacobian/fixed-sized-hypercube}.)
-However, by convention we do assume that $\verb|m|=0$ is always a valid
-\verb|m| coordinate; this simplifies pointer arithmetic.
-
-When describing various entries in the parameter table in
-section~\ref{localinterp/sect-generic-options}, we use \verb|const| qualifiers
-on table entries to indicate that the interpolator treats them as
-\verb|const| variables/arrays, \ie{} it promises not to change them.
-In contrast, table entries which are not shown as \verb|const| are
-considered read-write by the interpolator; it typically uses them
-to return outputs back to the caller.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\section{The Uniform Cartesian Interpolator}
-\label{localinterp/sect-uniform-Cartesian-interp}
-
-The uniform Cartesian interpolator is the one which used to live
-in the \thorn{PUGHInterp} thorn.  It was written in C by Thomas Radke
-in early 2001, drawing on earlier Fortran code by Paul Walker.  It
-supports the \verb|CCTK_InterpLocal()| API, and thus doesn't use
-any parameter table at all.  It provides Lagrange polynomial
-interpolation of orders 1, 2, and 3, registered under the operator
-names \verb|"first-order uniform cartesian"|,
-\verb|"second-order uniform cartesian"|, and\\
-\verb|"third-order uniform cartesian"| respectively.  Each of these
-supports 1, 2, and 3-D arrays.  The code is hard-wired for
-hypercube-shaped interpolation molecules, but it would be easy
-to add additional orders and/or dimensions if desired.
+It supports 1, 2, and 3-dimensional interpolation.
+At present this thorn only supports the \verb|CCTK_InterpLocal()| API,
+but in the near future it will probably be enhanced to support the
+newer \verb|CCTK_InterpLocalUniform()| API.
 
 Although the \verb|CCTK_InterpLocal()| API supports both uniform
 and nonuniform grids for the input data, the present implementation
@@ -331,7 +146,25 @@ API, and some examples of how to use it.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\subsection{Implementation}
+\subsection{History}
+
+This interpolator was written by Thomas Radke in early 2001 (drawing
+on older code by Paul Walker).  It originally lived in the PUGHInterp
+thorn, but it turned to have very little to do with PUGH, so was moved
+here in winter 2001-2002.
+
+From winter 2001-2002 to July 2003 this thorn also contained another
+interpolator written by Jonathan Thornburg, but in July~2003 that
+interpolator was moved to {\bf AEIThorns/AEILocalInterp/} because it
+was (is) GPL and Cactus policies forbids GPL code in this arrangement.
+
+Because CVS can't delete directories, this thorn still contains a lot
+of empty-except-for-CVS-directories directory trees left over from
+Jonathan Thornburg's interpolator.  You can/should ignore these.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Implementation}
 
 Internally, this interpolator always does 3-D interpolation, inserting
 zero coordinates as appropriate for lower dimensionalities.  The
@@ -350,1430 +183,6 @@ axis.%%%
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\section{Generalized Polynomial Interpolation}
-\label{localinterp/sect-generic-options}
-
-The generalized polynomial interpolator was written in C
-(plus Maple to generate the coefficients) by Jonathan Thornburg in
-winter 2001-2002.  It provides:
-\begin{itemize}
-\item	Lagrange polynomial interpolation of orders 1--6 for 1-D arrays,
-	and 1--4 for 2- and 3-D arrays.  This interpolator actually provides
-	two slightly different flavors of Lagrange interpolation,
-	\defn{tensor product} and \defn{maximum degree}, registered
-	under the operator names\\
-	\verb|"Lagrange polynomial interpolation (tensor product)"| and\\
-	\verb|"Lagrange polynomial interpolation (maximum degree)"|
-	respectively.  Section~\ref{localinterp/sect-multi-dim-interp}
-	of this thorn guide explains the two variants in detail, but
-	for most purposes the \underline{tensor-product} variant is
-	preferable.  For convenience this is also registered under
-	the operator name\\
-	\verb|"Lagrange polynomial interpolation"|.%%%
-\footnote{%%%
-	 For backwards compatability it's {\em also\/}
-	 registered under the operator name
-	 \hbox{\tt "generalized polynomial interpolation"}.
-	 }%%%
-\item	Hermite polynomial interpolation of orders 2-4 for 1-D arrays,
-	and 2--3 for 2- and 3-D arrays, registered under the operator
-	name \verb|"Hermite polynomial interpolation"|.%%%
-\end{itemize}
-The code allows arbitrarily-shaped interpolation molecules,
-but at the moment only hypercube-shaped molecules are implemented.
-It would be easy to add additional orders and/or dimensions if desired.
-
-This interpolator supports a number of options specified by a
-\defn{parameter table} (a Cactus key-value table, manipulated by the
-\verb|Util_Table*()| APIs).  Note that all the interpolator options
-described here apply only to the current interpolator call: there is
-no visible state kept inside the interpolator from one call to another.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Interpolation Order}
-
-The only mandatory parameter for this interpolator is the interpolation
-order:
-\begin{verbatim}
-const CCTK_INT order;
-\end{verbatim}
-As noted in section~\ref{localinterp/sect-basic-terminology},
-in our terminology linear interpolation is order~1, quadratic is
-order~2, cubic is order~3, etc.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Molecule Size and Centering}
-\label{localinterp/sect-molecule-size+centering}
-
-If no grid boundaries or excised points are nearby, the interpolator
-centers the molecules around the interpolation point as much as possible.
-Table~\ref{localinterp/tab-molecule-size+centering} gives the molecule
-size and details of the centering policy for each interpolation operator
-and order.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\begin{table}
-
-\def\Linterval{\lower0.60ex\hbox{\hskip-0.20em$\big[$}}
-\def\Rinterval{\lower0.60ex\hbox{\hskip-0.25em$\big)$}}
-
-\def\TwoPointMolecule{%%%
-\raise0.5ex\hbox{%%%
-\begin{picture}(0,0)
-\put(-4,0){\circle*{1.5}}
-\put(-4,0){\Linterval}
-\put(4,0){\Rinterval}
-\put(4,0){\circle*{1.5}}
-\end{picture}%%%
-}}%%%
-
-\def\ThreePointMolecule{%%%
-\raise0.5ex\hbox{%%%
-\begin{picture}(0,0)
-\put(-8,0){\circle*{1.5}}
-\put(-4,0){\Linterval}
-\put(0,0){\circle*{1.5}}
-\put(4,0){\Rinterval}
-\put(8,0){\circle*{1.5}}
-\end{picture}%%%
-}}%%%
-
-\def\FourPointMolecule{%%%
-\raise0.5ex\hbox{%%%
-\begin{picture}(0,0)
-\put(-12,0){\circle*{1.5}}
-\put(-4,0){\circle*{1.5}}
-\put(-4,0){\Linterval}
-\put(4,0){\Rinterval}
-\put(4,0){\circle*{1.5}}
-\put(12,0){\circle*{1.5}}
-\end{picture}%%%
-}}%%%
-
-\def\FivePointMolecule{%%%
-\raise0.5ex\hbox{%%%
-\begin{picture}(0,0)
-\put(-16,0){\circle*{1.5}}
-\put(-8,0){\circle*{1.5}}
-\put(-4,0){\Linterval}
-\put(0,0){\circle*{1.5}}
-\put(4,0){\Rinterval}
-\put(8,0){\circle*{1.5}}
-\put(16,0){\circle*{1.5}}
-\end{picture}%%%
-}}%%%
-
-\def\SixPointMolecule{%%%
-\raise0.5ex\hbox{%%%
-\begin{picture}(0,0)
-\put(-20,0){\circle*{1.5}}
-\put(-12,0){\circle*{1.5}}
-\put(-4,0){\circle*{1.5}}
-\put(-4,0){\Linterval}
-\put(4,0){\Rinterval}
-\put(4,0){\circle*{1.5}}
-\put(12,0){\circle*{1.5}}
-\put(20,0){\circle*{1.5}}
-\end{picture}%%%
-}}%%%
-
-\def\SevenPointMolecule{%%%
-\raise0.5ex\hbox{%%%
-\begin{picture}(0,0)
-\put(-24,0){\circle*{1.5}}
-\put(-16,0){\circle*{1.5}}
-\put(-8,0){\circle*{1.5}}
-\put(-4,0){\Linterval}
-\put(0,0){\circle*{1.5}}
-\put(4,0){\Rinterval}
-\put(8,0){\circle*{1.5}}
-\put(16,0){\circle*{1.5}}
-\put(24,0){\circle*{1.5}}
-\end{picture}%%%
-}}%%%
-
-\renewcommand{\arraystretch}{1.5}
-\begin{flushleft}
-\begin{tabular}{@{}lcc@{\quad\hspace{25mm}}c@{\hspace{25mm}}}
-						&	& Molecule
-	&								\\[-1ex]
-Interpolation Operator				& Order	& Size
-	& \csmash{Molecule Centering Policy}				\\
-\hline %-------------------------------------------------------
-\verb|"Lagrange polynomial interpolation"|	& 1	& 2-point
-	& \TwoPointMolecule						\\
-\verb|"Lagrange polynomial interpolation"|	& 2	& 3-point
-	& \ThreePointMolecule						\\
-\verb|"Lagrange polynomial interpolation"|	& 3	& 4-point
-	& \FourPointMolecule						\\
-\verb|"Lagrange polynomial interpolation"|	& 4	& 5-point
-	& \FivePointMolecule						\\
-\verb|"Lagrange polynomial interpolation"|	& 5	& 6-point
-	& \SixPointMolecule						\\
-\verb|"Lagrange polynomial interpolation"|	& 6	& 7-point
-	& \SevenPointMolecule						\\
-\hline %-------------------------------------------------------
-\verb|"Hermite polynomial interpolation"|	& 2	& 4-point
-	& \FourPointMolecule						\\
-\verb|"Hermite polynomial interpolation"|	& 3	& 6-point
-	& \SixPointMolecule						\\
-\verb|"Hermite polynomial interpolation"|	& 4	& 6-point
-	& \SixPointMolecule						\\
-\hline %-------------------------------------------------------
-\end{tabular}
-\end{flushleft}
-\caption[Molecule Sizes]
-	{
-	This table gives the molecule size and centering for each
-	interpolation operator and order.  $\bullet$~shows the input
-	grid points, and $[ \quad )$~shows the interval of interpolation
-	points (relative to the grid) for which the interpolator
-	will use the molecule shown.
-	For example, with grid points at integer coordinates, if
-	the the interpolator is using 5-point molecules, it will
-	use a molecule containing the grid points $\{ -2, -1, 0, 1, 2 \}$
-	for all interpolation points in the range $[-0.5, 0.5)$.
-	}
-\label{localinterp/tab-molecule-size+centering}
-\end{table}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Handling of Grid Boundaries}
-
-Near grid boundaries and/or excised points the interpolator can either
-off-center the interpolation molecules, or refuse to interpolate
-(returning a \verb|CCTK_INTERP_ERROR_POINT_OUTSIDE|\,%%%
-\footnote{%%%
-	 For backwards compatability the error code
-	 {\tt CCTK\_INTERP\_ERROR\_POINT\_X\_RANGE}
-	 is also defined; it's a synonym for
-	 {\tt CCTK\_INTERP\_ERROR\_POINT\_OUTSIDE}.
-	 }%%%
-{} or \verb|CCTK_INTERP_ERROR_POINT_EXCISED| error code).  This behavior
-is controlled by the (optional) parameter-table entries
-\begin{verbatim}
-const CCTK_INT N_boundary_points_to_omit[2*N_dims]
-const CCTK_REAL boundary_off_centering_tolerance[2*N_dims]
-const CCTK_REAL boundary_extrapolation_tolerance[2*N_dims]
-const CCTK_REAL excision_off_centering_tolerance[2*N_dims] # not implemented yet
-const CCTK_REAL excision_extrapolation_tolerance[2*N_dims] # not implemented yet
-\end{verbatim}
-The elements of these arrays are matched up with the grid boundaries
-in the order\\
-$[x_{\min}, x_{\max}, y_{\min}, y_{\max}, z_{\min}, z_{\max}, \dots]$.
-(For \verb|excision_*_tolerance[]| the minimum/maximum are interpreted
-as the corresponding coordinate decreasing/increasing when moving
-out of the active region of the grid into the excised region.)
-We use \verb|ibndry| as a generic 0-origin index into these arrays.
-
-%%%%%%%%%%%%%%%%%%%%
-
-\subsubsection{Omitting Boundary Points}
-
-For any given grid boundary \verb|ibndry|, the interpolator doesn't use
-input data from the\\
-\verb|N_boundary_points_to_omit[ibndry]| grid points closest to the
-boundary.  In other words, these points are effectively omitted from
-the input grid.  For example:
-\begin{itemize}
-\item	Setting this parameter to an array of all 0s means the
-	interpolator may use data from all grid points.
-	This is the default if this parameter isn't specified.
-\item	Setting this parameter to an array of all 1s means the
-	interpolator may not use input data from the outermost
-	row/plane of grid points on any face.
-\end{itemize}
-
-More generally, this option may be useful for controlling the extent
-to which the interpolator uses input data from ghost and/or symmetry
-zones.
-
-%%%%%%%%%%%%%%%%%%%%
-
-\subsubsection{Off-Centering Molecules near Grid Boundaries
-	       and/or Excision Regions}
-
-The (optional) parameter-table entries
-\begin{verbatim}
-const CCTK_REAL boundary_off_centering_tolerance[2*N_dims]
-const CCTK_REAL boundary_extrapolation_tolerance[2*N_dims]
-const CCTK_REAL excision_off_centering_tolerance[2*N_dims] # not implemented yet
-const CCTK_REAL excision_extrapolation_tolerance[2*N_dims] # not implemented yet
-\end{verbatim}
-control the local interpolator's willingness to off-center
-(change away from the default centering) its interpolation molecules
-near boundaries and excised points.
-
-To describe how these parameters work, it's useful to define two
-regions in the data coordinate space:
-\begin{description}
-\item[The \defn{valid-data region}]
-	is the entire grid minus any ommited-on-boundary
-	or excised points; to make this a region we take the
-	Lego-block bounding box of the non-excised grid points.
-\item[The \defn{default-centering region}]
-	is that region in interpolation-point space where the
-	interpolator can use the default molecule centering
-	(described in table~\ref{localinterp/tab-molecule-size+centering}),
-	\ie{} where the default-centering molecules require data
-	only from the data-valid region.
-\end{description}
-Figure~\ref{localinterp/fig-valid-data+default-centering} shows an
-example of these regions.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\begin{figure}
-
-% "X" to mark excised point
-\def\ExcisedPoint{%%%
-\begin{picture}(0,0)
-\put(-2,-2){\line(1,1){4}}
-\put(2,-2){\line(-1,1){4}}
-\end{picture}%%%
-}%%%
-
-\def\Linterval{\lower0.60ex\hbox{\hskip-0.20em$\big[$}}
-
-\def\xmin{\begin{picture}(0,0)\put(4,0){\vector(-1,0){4}}\end{picture}}
-\def\xmax{\begin{picture}(0,0)\put(-4,0){\vector(1,0){4}}\end{picture}}
-\def\ymin{\begin{picture}(0,0)\put(0,4){\vector(0,-1){4}}\end{picture}}
-\def\ymax{\begin{picture}(0,0)\put(0,-4){\vector(0,1){4}}\end{picture}}
-
-\begin{center}
-\begin{picture}(150,150)
-%%% grid
-\multiput(0,0)(10,0){16}{\multiput(0,0)(0,10){16}{\circle*{1.5}}}
-%%% excised points
-\multiput(70,60)(10,0){2}{\ExcisedPoint}
-\multiput(60,70)(10,0){4}{\ExcisedPoint}
-\multiput(60,80)(10,0){4}{\ExcisedPoint}
-\multiput(70,90)(10,0){2}{\ExcisedPoint}
-%%% outer-boundary valid-data region
-\multiput(0,0)(2,0){75}{\line(1,0){1}}
-\multiput(150,0)(0,2){75}{\line(0,1){1}}
-\multiput(150,150)(-2,0){75}{\line(-1,0){1}}
-\multiput(0,150)(0,-2){75}{\line(0,-1){1}}
-% ... and movement directions
-\put(75,10){\ymin}
-\put(140,75){\xmax}
-\put(75,140){\ymax}
-\put(10,75){\xmin}
-%%% outer-boundary default-centering regions...
-% ... for 3-point molecules
-\put(5,5){\line(1,0){140}}
-\put(145,5){\line(0,1){140}}
-\put(145,145){\line(-1,0){140}}
-\put(5,145){\line(0,-1){140}}
-% ... for 4-point molecules
-\put(10,10){\line(1,0){130}}
-\put(140,10){\line(0,1){130}}
-\put(140,140){\line(-1,0){130}}
-\put(10,140){\line(0,-1){130}}
-%%% excised-point valid-data region
-\multiput(60,60)(0,-2){5}{\line(0,-1){1}}
-\multiput(60,50)(2,0){15}{\line(1,0){1}}
-\multiput(90,50)(0,2){5}{\line(0,1){1}}
-\multiput(90,60)(2,0){5}{\line(1,0){1}}
-\multiput(100,60)(0,2){15}{\line(0,1){1}}
-\multiput(100,90)(-2,0){5}{\line(-1,0){1}}
-\multiput(90,90)(0,2){5}{\line(0,1){1}}
-\multiput(90,100)(-2,0){15}{\line(-1,0){1}}
-\multiput(60,100)(0,-2){5}{\line(0,-1){1}}
-\multiput(60,90)(-2,0){5}{\line(-1,0){1}}
-\multiput(50,90)(0,-2){15}{\line(0,-1){1}}
-\multiput(50,60)(2,0){5}{\line(1,0){1}}
-%%% excised-point default-centering regions...
-% ... for 3-point molecules
-\put(55,55){\line(0,-1){10}}
-\put(55,45){\line(1,0){40}}
-\put(95,45){\line(0,1){10}}
-\put(95,55){\line(1,0){10}}
-\put(105,55){\line(0,1){40}}
-\put(105,95){\line(-1,0){10}}
-\put(95,95){\line(0,1){10}}
-\put(95,105){\line(-1,0){40}}
-\put(55,105){\line(0,-1){10}}
-\put(55,95){\line(-1,0){10}}
-\put(45,95){\line(0,-1){40}}
-\put(45,55){\line(1,0){10}}
-% ... for 4-point molecules
-\put(50,50){\line(0,-1){10}}
-\put(50,40){\line(1,0){50}}
-\put(100,40){\line(0,1){10}}
-\put(100,50){\line(1,0){10}}
-\put(110,50){\line(0,1){50}}
-\put(110,100){\line(-1,0){10}}
-\put(100,100){\line(0,1){10}}
-\put(100,110){\line(-1,0){50}}
-\put(50,110){\line(0,-1){10}}
-\put(50,100){\line(-1,0){10}}
-\put(40,100){\line(0,-1){50}}
-\put(40,50){\line(1,0){10}}
-% ... and movement directions for 4-point molecules
-\put(50,45){\xmax}
-\put(75,40){\ymax}
-\put(100,45){\xmin}
-\put(105,50){\ymax}
-\put(110,75){\xmin}
-\put(105,100){\ymin}
-\put(100,105){\xmin}
-\put(75,110){\ymin}
-\put(50,105){\xmax}
-\put(45,100){\ymin}
-\put(40,75){\xmax}
-\put(45,50){\ymax}
-\end{picture}
-\end{center}
-
-\caption[Example of valid-data and default-centering regions]
-	{
-	This figure shows an example of the valid-data and
-	default-centering regions for a 2-D grid.
-	$\times$~marks the excised points, and the dashed lines
-	show the boundaries of the valid-data region.
-	The solid lines show the boundaries of the default-centering
-	regions for 3-point and 4-point molecules.  The arrows on the
-	default-centering region for 4-point molecules show which
-	elements of the \hbox{\tt boundary\_*\_tolerance[ibndry]}
-	and \hbox{\tt excision\_*\_tolerance[ibndry]} arrays apply
-	to each line segment of the region boundary:\\
-\centerline{%%%
-\renewcommand{\arraystretch}{1.25}
-\begin{tabular}{c@{\qquad}r}
-$\raise0.50ex\hbox{\xmin \hskip4mm}$	& $x_{\min}$ ({\tt ibndry = 0})	\\
-$\raise0.50ex\hbox{\hskip4mm \xmax}$	& $x_{\max}$ ({\tt ibndry = 1})	\\
-$\lower1mm\hbox{\,\ymin\,}$		& $y_{\min}$ ({\tt ibndry = 2})	\\
-$\raise3mm\hbox{\,\ymax\,}$		& $y_{\max}$ ({\tt ibndry = 3})	%%%\\
-\end{tabular}%%%
-}
-	}
-\label{localinterp/fig-valid-data+default-centering}
-\end{figure}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-Near a boundary or excised regions, the interpolator will off-center
-its molecules (if and) only if the interpolation point is both
-\begin{itemize}
-\item	within ($\le$) \verb|*_off_centering_tolerance[ibndry]|
-	grid spacings of the default-centering region, {\bf and}
-\item	within ($\le$) \verb|*_extrapolation_tolerance[ibndry]|
-	grid spacings of the valid-data region.
-\end{itemize}
-where we use ``\verb|*|'' to denote \verb|boundary| or \verb|excision|
-as appropriate.
-
-There are four cases for these parameters:
-\begin{description}
-\item[\mathversion{bold}
-      $\text{\tt *\!\_off\_centering\_tolerance[ibdnry]} = 0.0$ and
-      $\text{\tt *\!\_extrapolation\_tolerance[ibndry]} = 0.0$:]\mbox{}\\
-	No off-centering is allowed: the interpolator reports an
-	error (\verb|CCTK_ERROR_INTERP_POINT_OUTSIDE| or
-	\verb|CCTK_ERROR_INTERP_POINT_EXCISED| return code
-	as appropriate) if any interpolation point is outside
-	the default-centering region.
-\item[\mathversion{bold}
-      $\text{\tt *\!\_off\_centering\_tolerance[ibdnry]} > 0.0$ and
-      $\text{\tt *\!\_extrapolation\_tolerance[ibndry]} = 0.0$:]\mbox{}\\
-	The interpolator allows interpolation points to be up to ($\le$)
-	\verb|*_off_centering_tolerance[ibndry]| grid spacings outside
-	the default-centering region in this direction.
-	If an interpolation point is beyond this limit in any
-	direction, then the interpolator reports an error.
-\item[\mathversion{bold}
-      $\text{\tt *\!\_off\_centering\_tolerance[ibdnry]} = \infty$%%%
-\footnotemark{}%%%
-{}    and
-      $\text{\tt *\!\_extrapolation\_tolerance[ibndry]} > 0.0$:]\mbox{}\\
-\footnotetext{%%%
-	     In practice $999.0$ is a conventional value here.
-	     The actual value doesn't matter so long as it's
-	     larger than any possible molecule radius.
-	     }%%%
-	The interpolator allows interpolation points to be up to ($\le$)
-	\verb|*_extrapolation_tolerance[ibndry]| grid spacings outside
-	the valid-data region in this direction.
-	If an interpolation point is beyond this limit in any
-	direction, then the interpolator reports an error.
-\item[\mathversion{bold}
-      $\text{\tt *\!\_off\_centering\_tolerance[ibdnry]} = 0.0$ and
-      $\text{\tt *\!\_extrapolation\_tolerance[ibndry]} > 0.0$:]\mbox{}\\
-	In practice the default-centering region is always a
-	(normally proper) subset of the valid-data region, so this
-	case is nonsensical: the positive value for
-	\verb|*_extrapolation_tolerance[ibndry]| has no effect
-	because of the $\verb|*_off_centering_tolerance[ibndry]| = 0.0$
-	setting.  The interpolator gives a warning for this case.
-\end{description}
-
-If any of these table entries aren't specified, the defaults are
-\begin{verbatim}
-boundary_off_centering_tolerance[ibndry] = 999.0   # effectively infinity
-boundary_extrapolation_tolerance[ibndry] = 1.0e-10
-\end{verbatim}
-In other words, the interpolation points may be anywhere within the
-valid-data region or up to $10^{-10}$ grid spacing outside it.  (The
-$10^{-10}$ ``fudge factor'' helps to avoid spurious errors in case
-floating-point roundoff moves an interpolation point which was supposed
-to be just on the boundary of the valid-data region, slightly outside it.)
-
-[In the near future these defaults may be changed:  For Lagrange
-polynomial interpolation the new defaults would be
-\begin{verbatim}
-boundary_off_centering_tolerance[ibndry] = 999.0   # effectively infinity
-boundary_extrapolation_tolerance[ibndry] = 1.0e-10
-\end{verbatim}
-while for Hermite polynomial interpolation they would be
-\begin{verbatim}
-boundary_off_centering_tolerance[ibndry] = 1.0e-10
-boundary_extrapolation_tolerance[ibndry] = 0.0
-\end{verbatim}
-This would leave Lagrange interpolation unchanged, while for Hermite
-interpolation the defaults would forbid any significant off-centering
-of the interpolation molecules.]
-
-%%%%%%%%%%%%%%%%%%%%
-
-\subsubsection{Out-of-Range Interpolation Point Error Handling}
-
-If the interpolator finds an interpolation point which is ``out of range''
-as described above, it sets the following parameter-table entries to
-give more information about the out-of-range point:
-\begin{verbatim}
-/* which interpolation point is it? */
-/* (value gives 0-origin  pt  for the offending point) */
-CCTK_INT error_pt;
-
-/* in which coordinate axis and direction is the point out of range? */
-/* (value gives 0-origin  ibndry  for the offending point) */
-CCTK_INT error_ibndry;
-
-/* in which coordinate axis is the point outside the grid or excised? */
-/* (value gives 0-origin  axis  for the offending point) */
-CCTK_INT error_axis;
-
-/* in which direction is the point out of range? */
-/* (value is -1 for min, +1 for max) */
-CCTK_INT error_direction;
-\end{verbatim}
-
-The interpolator then returns a \verb|CCTK_ERROR_INTERP_POINT_OUTSIDE|
-error code.
-
-Note that if the point is out of range in multiple axes, it's undefined
-which of them will be reported.  Also, if there are multiple out-of-range
-points, user code shouldn't make any assumptions about which of them
-will be reported -- the only safe assumption is that {\em some\/}
-out-of-range point will be reported.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Multi-Dimensional Interpolation}
-\label{localinterp/sect-multi-dim-interp}
-
-In multiple dimensions, there are two plausible definitions for the
-generic interpolating polynomial of a given order (degree) $n$.  For
-convenience I'll describe these for the 2-D case, but the generalization
-to any number of dimensions should be obvious:
-
-\begin{itemize}
-\item	The generic ``\defn{tensor product}'' polynomial of degree $n$ is
-	\begin{equation}
-	f(x,y) = \sum_{{0 \le i \le n} \atop {0 \le j \le n}} a_{ij} x^i y^j
-					\label{localinterp/eqn-polynomial-2d-TP}
-	\end{equation}
-\item	The generic ``\defn{maximum degree}'' polynomial of degree $n$ is
-	\begin{equation}
-	g(x,y) = \sum_{0 \le i+j \le n} a_{ij} x^i y^j
-					\label{localinterp/eqn-polynomial-2d-MD}
-	\end{equation}
-\end{itemize}
-
-Because it has $(n+1)^2$ coefficients, the tensor-product polynomial $f$
-can (and does) pass through all the $(n+1)^2$ input data points in a
-square molecule of size $n+1$.  This implies that the interpolation
-error vanishes at the input grid points, and that the overall
-interpolating function is continuous (up to floating-point roundoff
-errors).  However, $f$ does have the slightly peculiar property of
-having terms up to $x^n y^n$ despite being of ``degree'' $n$.
-(For example, the ``linear'' interpolant for $n=1$ would have $xy$ terms,
-even though those are formally quadratic in the independent variables
-$x$ and $y$.)
-
-In contrast, the maximum-degree polynomial $g$ has only
-$\frac{1}{2} (n+1)(n+2)$ coefficients, so for generic input data
-(\ie{} input data which isn't actually sampled from a polynomial of
-the form~$(\ref{localinterp/eqn-polynomial-2d-MD})$) $g$ can't pass
-through all the $(n+1)^2$ input data points in a square molecule of
-size $n+1$.  The interpolator actually does a least-squares fit of
-the polymomial $g$ to the input data in the molecule, so in general
-the molecule won't pass through {\em any\/} of the data points!
-Moreover, each time the interpolation point crosses a grid square
-(for odd~$n$) or the center lines of a grid square (for even~$n$),
-the set of points used in the molecule changes 
-(\cf{}~section~\ref{localinterp/sect-molecule-size+centering}),
-so the interpolant generally has a jump discontinuity!
-For these reasons, the tensor-product
-choice~$(\ref{localinterp/eqn-polynomial-2d-TP})$ is generally
-preferable.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Molecule Family}
-\label{localinterp/sect-molecule-family}
-
-An interpolator may support different families/shapes of interpolation
-molecules.  Hypercube-shaped molecules are the simplest and most common
-case, but one could also imagine (say) octagon-shaped molecules in 2-D,
-or some generalization of this in higher numbers of dimensions.
-
-The (optional) parameter
-\begin{verbatim}
-/* set with Util_TableSetString() */
-const char molecule_family[];
-\end{verbatim}
-may be used to set or query the molecule family.
-
-If this key is present in the parameter table, the interpolator sets
-the molecule family/shape based on the value specified.
-If this key {\em isn't\/} present in the parameter table, then the
-interpolator sets it to the molecule family being used.
-
-At present only hypercubed-shaped molecules are implemented; these
-are specified by \verb|molecule_family| set to \verb|"cube"|.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Non-Contiguous Input Arrays}
-\label{localinterp/sect-non-contiguous-inputs}
-
-Sometimes the input ``arrays'' used by the interpolator might not
-be contiguous in memory.  For example, we might want to do 2-D interpolation
-within a plane of a 3-D grid array, but the plane might or might not
-be contiguous in memory.  (Another example would be that the input
-arrays might be members of a compact group.)
-
-The following (optional) parameter-table entries specify non-contiguous
-input arrays:
-\begin{verbatim}
-const CCTK_INT input_array_offsets[N_input_arrays];
-
-/* the next 3 table entries are shared by all input arrays */
-const CCTK_INT input_array_strides       [N_dims];
-const CCTK_INT input_array_min_subscripts[N_dims];
-const CCTK_INT input_array_max_subscripts[N_dims];
-\end{verbatim}
-
-In general, the interpolator accesses the input using the generic
-subscripting expression
-\begin{verbatim}
-input_array[in][offset + i*stride_i + j*stride_j + k*stride_k]
-\end{verbatim}
-where
-\begin{verbatim}
-offset = input_array_offsets[in]
-(stride_i,stride_j,stride_k) = input_array_strides[]
-\end{verbatim}
-and where \verb|(i,j,k)| run from \verb|input_array_min_subscripts[]|
-to \verb|input_array_max_subscripts[]| inclusive
-(n.b.~this is an {\em inclusive\/} range, \ie{} 
-$\verb|min| \le \verb|(i,j,k)| \le \verb|max|$).
-
-The defaults are that each input array is contiguous in memory,
-\ie{} \verb|input_array_offsets[]| = 0,
-\verb|stride| determined from \verb|input_array_dims[]|
-              in the usual Fortran manner,
-\verb|input_array_min_subscripts[]| = all 0s, and
-\verb|input_array_max_subscripts[]| = \verb|input_array_dims[]|-1.
-If the stride and max subscript are both specified explicitly, then the
-explicit \verb|input_array_dims[]| argument to
-\verb|CCTK_InterpLocalUniform()| is ignored.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Derivatives}
-\label{localinterp/sect-derivatives}
-
-If we view the input data as being samples of a smooth function,
-then instead of estimating values of that function at the interpolation
-points, the interpolator can instead or additionally estimate values
-of various derivatives of that function at the interpolation points.
-(We don't currently implement it, but one could also imagine
-interpolating more general molecule operations such as Laplacians.)
-
-The following (optional) parameter-table entries are used to specify
-this:
-\begin{verbatim}
-const CCTK_INT operand_indices[N_output_arrays];
-const CCTK_INT operation_codes[N_output_arrays];
-\end{verbatim}
-The semantics here are that
-\begin{verbatim}
-output_array[out] = op(input_array[operand_indices[out]])
-\end{verbatim}
-where \verb|op| is an operator specified by the \verb|operation_codes[out]|
-value as described below.
-
-Note that \verb|operand_indices[]| doesn't directly name the inputs,
-but rather gives their (0-origin) subscripts in the list of inputs.
-This allows for a more efficient implementation in the (common) case
-where some of the input arrays have many different operations applied
-to them.  (It's most efficient to group all operations on a given
-input array together in the \verb|operand_indices| and
-\verb|operation_codes| arrays, as in the example in
-section~\ref{localinterp/sect-example-derivatives}.)
-
-Negative \verb|operation_codes[out]| values  are reserved for future
-use.  An \verb|operation_codes[out]| value which is $\ge 0$ is taken
-as a decimal integer encoding a coordinate partial derivative: each
-decimal digit means to take the coordinate partial derivative along
-that (1-origin) axis; the order of the digits in a number is ignored.
-Table~\ref{localinterp/tab-derivative-codes} summarizes these resulting
-derivative codes.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\begin{table}
-\begin{center}
-\begin{tabular}{cl}
-\verb|operation_codes[out]|
-	& What it means							\\
-\hline %---------------------------------------------------------------
-0	& interpolate the input array itself (no derivative)		\\[1ex]
-1	& interpolate $\partial \big/ \partial x^1$ of the input array	\\
-2	& interpolate $\partial \big/ \partial x^2$ of the input array	\\
-3	& interpolate $\partial \big/ \partial x^3$ of the input array	\\[1ex]
-12 or 21& interpolate $\partial^2 \big/ \partial x^1 \partial x^2$
-	  of the input array						\\
-13 or 31& interpolate $\partial^2 \big/ \partial x^1 \partial x^3$
-	  of the input array						\\
-23 or 32& interpolate $\partial^2 \big/ \partial x^2 \partial x^3$
-	  of the input array						\\
-11	& interpolate $\partial^2 \big/ \partial (x^1)^2$ of the input array \\
-22	& interpolate $\partial^2 \big/ \partial (x^2)^2$ of the input array \\
-33	& interpolate $\partial^2 \big/ \partial (x^3)^2$ of the input array \\
-\hline %---------------------------------------------------------------
-\end{tabular}
-\end{center}
-\caption[Derivative Codes]
-	{
-	This table gives the codes in {\tt operation\_codes[out]}
-	for each possible 1st or 2nd~derivative in 3-D; for 1-D or
-	2-D codes referring to nonexistent coordinates are invalid.
-	}
-\label{localinterp/tab-derivative-codes}
-\end{table}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-At present we do {\em not\/} have any \verb|#define|s for the
-operation codes in the Cactus header files.  However, you can avoid
-most of the software-engineering problems of having ``magic numbers''
-for the operation codes, by using the macro
-\begin{verbatim}
-#define DERIV(op)   op
-\end{verbatim}
-to mark all such \verb|operation_codes[]| values in your code.
-There's an example of this in
-section~\ref{localinterp/sect-example-derivatives}.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Jacobian and Domain of Dependence}
-\label{localinterp/sect-Jacobian}
-
-The Jacobian of the interpolator is defined as
-\begin{equation}
-\frac{\partial \hbox{\tt output\_array[out][pt]}}
-     {\partial \hbox{\tt input\_array[in][(i,j,k)]}}
-						\label{localinterp/eqn-Jacobian}
-\end{equation}
-We may want to know the Jacobian itself, and/or ``just'' the set of
-\verb|(i,j,k)| for which this is nonzero (\ie{} the Jacobian's sparsity
-structure, or equivalently the domain of dependence of the output result,
-or equivalently the interpolation molecule size and shape) for a given
-\verb|out|, \verb|in|, and \verb|pt|.%%%
-\footnote{%%%
-	 For something like a spline interpolator the Jacobian
-	 is generally nonzero for all {\tt (i,j,k)}, but
-	 exponentially small for most of them; in this case
-	 the Jacobian-query API would probably specify a cutoff
-	 to return an approximate Jacobian with reasonable sparsity.
-	 }%%%
-
-The complexity of doing this depends (strongly!) on the structure
-of the Jacobian, and in particular on the answers to the following
-questions:
-\begin{itemize}
-\item	What molecule family is being used?
-\item	Does the interpolation molecule size and/or shape depend
-	on where the interpolation points are in the grid?%%%
-\footnote{%%%
-	 We can always make the interpolation molecules be
-	 the ``same size and shape'' by padding them with
-	 zero entries to bring them all up to the size of
-	 the largest molecule used anywhere in the grid,
-	 but this would be very inefficient if many molecules
-	 were padded in this way.
-	 }%%%
-\item	If this interpolator supports computing derivatives
-	as described in section~\ref{localinterp/sect-derivatives},
-	does the interpolation molecule size and/or shape depend
-	on \verb|operation_codes[]|?
-\item	Does the interpolation molecule size and/or shape depend
-	on the actual floating-point values being interpolated?
-	(Examples of this might include ENO (essentially nonoscillatory)
-	and/or TVD (total-variation diminishing) interpolators for
-	hydrodynamics calculations.)
-\item	Do the actual floating-point values of the Jacobian depend
-	on the actual floating-point values being interpolated?
-	Equivalently, is the interpolation nonlinear?
-\end{itemize}
-
-Because the different cases differ so much in the amount of information
-required to describe the Jacobian, it's hard to define a single API
-which covers all cases without burdening the simpler cases with
-excessive complexity.  Instead, the interpolator defines a
-Jacobian-structure--query API to determine which case applies,
-together with (conceptually) several different APIs for the
-different cases.  (At the moment only a single case is implemented.)
-
-%%%%%%%%%%%%%%%%%%%%
-
-\subsubsection{Determining the Jacobian's structure}
-\label{localinterp/sect-Jacobian/structure}
-
-The following parameter-table entries may be used to query which
-of the different Jacobian-structure cases applies:
-
-The parameter \verb|molecule_family| may may be used to query what
-molecule family is being used.  This is described in detail in
-section~\ref{localinterp/sect-molecule-family}.
-
-If the interpolation molecule size and/or shape vary with the
-interpolation coordinates, the interpolator sets the parameter
-\begin{verbatim}
-CCTK_INT MSS_is_fn_of_interp_coords;
-\end{verbatim}
-to 1.  Otherwise (\ie{} if the interpolation molecule size and shape
-are independent of the interpolation coordinates) it should set this
-parameter to 0.
-
-If the interpolator supports computing derivatives as described in
-section~\ref{localinterp/sect-derivatives}, {\em and\/}
-if the interpolation molecule's size and/or shape varies with
-\verb|operation_codes[]|, the interpolator sets the parameter
-\begin{verbatim}
-CCTK_INT MSS_is_fn_of_which_operation;
-\end{verbatim}
-to 1.  Otherwise (\ie{} if the interpolator doesn't support computing
-derivatives, or if the interpolator does support computing derivatives
-but the interpolation molecule size and shape are independent of the
-\verb|operation_code[]| values), the interpolator sets this parameter
-to 0.  Note that this query tests whether the molecule size and/or shape
-depend on \verb|operation_codes[]| in general, independent of whether
-there are in fact any distinct values (or even any values at all) passed
-in \verb|operation_codes[]| in this particular interpolator call.  In
-other words, this is a query about the basic design of the interpolator,
-not about this particular call.
-
-If the interpolation molecule's size and/or shape varies with the
-actual floating-point values of the input arrays, the interpolator
-sets the parameter
-\begin{verbatim}
-CCTK_INT MSS_is_fn_of_input_array_values;
-\end{verbatim}
-to 1.  Otherwise (\ie{} if the interpolation molecule size and shape
-are independent of the input array values; this is a necessary, but
-not sufficient, condition for the interpolation to be linear), the
-interpolator sets this parameter to 0.
-
-If the actual floating-point values of the
-Jacobian~$(\ref{localinterp/eqn-Jacobian})$
-(for a given \verb|out|, \verb|in|, and \verb|pt|) depend on the actual
-floating-point values of the input arrays (\ie{} if the interpolation
-is nonlinear), the interpolator sets the parameter
-\begin{verbatim}
-CCTK_INT Jacobian_is_fn_of_input_array_values;
-\end{verbatim}
-to 1.  Otherwise (\ie{} if the interpolation is linear) the interpolator
-sets this parameter to 0.
-
-The current implementation always sets
-\begin{verbatim}
-MSS_is_fn_of_interp_coords = 0
-MSS_is_fn_of_which_operation = 0
-MSS_is_fn_of_input_array_values = 0
-Jacobian_is_fn_of_input_array_values = 0
-\end{verbatim}
-
-%%%%%%%%%%%%%%%%%%%%
-
-\subsubsection{Fixed-Size Hypercube-Shaped Molecules}
-\label{localinterp/sect-Jacobian/fixed-sized-hypercube}
-
-The simplest case (and the only one for which we have defined an
-API at present) is when the molecules are hypercube-shaped and of
-(typically small) fixed size, independent of the interpolation
-coordinates and the actual floating-point values in the input arrays
-(though presumably depending on the interpolation order and on
-\verb|operation_code|).  In other words, this case applies if
-(and only if) the Jacobian structure information described in
-section~\ref{localinterp/sect-Jacobian/structure}
-returns
-\begin{verbatim}
-MSS_is_fn_of_interp_coords = 0
-MSS_is_fn_of_which_operation = 0
-MSS_is_fn_of_input_array_values = 0
-Jacobian_is_fn_of_input_array_values = 0
-\end{verbatim}
-(These are precisely the values set by the current implementation.)
-In the rest of this section we describe the query API for this case.
-
-The following parameters may be used to query the molecule size:
-\begin{verbatim}
-CCTK_INT const molecule_min_m[N_dims];
-CCTK_INT const molecule_max_m[N_dims];
-\end{verbatim}
-The semantics of these are that if both of these keys are present
-(the values and datatypes don't matter), then the interpolator will
-(re)set the values to give the (inclusive) minimum and maximum
-\verb|m|~molecule coordinates.  (Note that either both of these
-keys should be present, or neither should be present.  This
-simplifies the logic in the interpolator slightly.)
-
-The following parameter may be used to query the molecule positions:
-\begin{verbatim}
-CCTK_INT* const molecule_positions[N_dims];
-\end{verbatim}
-The semantics of this is that the caller should set
-\verb|molecule_positions[]| to an array of \verb|N_dims| pointers to
-(caller-supplied) arrays of \verb|N_interp_points| \verb|CCTK_INT|s
-each.  If this key exists, then the interpolator will store the
-molecule positions in the pointed-to arrays.
-
-The following parameters may be used to query the
-Jacobian~$(\ref{localinterp/eqn-Jacobian})$ itself:
-\begin{verbatim}
-CCTK_REAL* const Jacobian_pointer[N_output_arrays];
-const CCTK_INT   Jacobian_offset [N_output_arrays];   # optional, default=all 0
-
-/* the next 3 table entries are shared by all Jacobians */
-const CCTK_INT Jacobian_interp_point_stride;
-const CCTK_INT Jacobian_m_strides[N_dims];
-const CCTK_INT Jacobian_part_stride;                  # optional, default=1
-\end{verbatim}
-If \verb|Jacobian_pointer| is present in the table, then
-\verb|Jacobian_interp_point_stride| and \verb|Jacobian_m_strides|
-must also be present.  For each \verb|out| where
-\verb|Jacobian_pointer[out] != NULL|,%%%
-\footnote{%%%
-	 That is, setting {\tt Jacobian\_pointer[out] = NULL}
-	 supresses the Jacobian query.  If (as is often the
-	 case) the Jacobian is independent of {\tt out}, then
-	 it's a good idea to avoid unnecessary computations by
-	 supressing the queries in this manner for any remaining
-	 output arrays.
-	 }%%%
-{} the interpolator would then store the
-Jacobian~$(\ref{localinterp/eqn-Jacobian})$ in
-\begin{verbatim}
-Jacobian_pointer[out][offset
-                      + pt*Jacobian_interp_point_stride
-                      + mi*stride_i + mj*stride_j + mk*stride_k
-                      + part*Jacobian_part_stride]
-\end{verbatim}
-where
-\begin{verbatim}
-offset = Jacobian_offset[out]
-(stride_i,stride_j,stride_k) = Jacobian_m_strides[]
-\end{verbatim}
-and where \verb|part| is 0 for real values and the real parts of complex
-values, and 1 for the imaginary parts of complex values.
-
-By appropriately setting the various stride parameters, this allows
-a fairly wide variety of possible storage layouts for the Jacobian.
-
-An example may help to clarify this:  Suppose we have a 1-D grid
-with 11~grid points, with integer subscripts~0 through~10 inclusive,
-and interpolation coordinates given by \verb|coord_origin = 0.0|
-and \verb|coord_delta = 0.1|.
-Suppose further that we're doing Lagrange polynomial interpolation,
-with \verb|order = 2| and hence
-(by table~\ref{localinterp/tab-molecule-size+centering})
-using 3-point molecules.
-Finally, suppose that we query the Jacobian molecule positions for
-the \verb|N_interp_points=14| interpolation points 0.0, 0.04, 0.06,
-0.10, 0.14, 0.16, 0.20, 0.80, 0.84, 0.86, 0.90, 0.94, 0.96, 1.00.
-Then the queries might return
-\begin{verbatim}
-interp_molecule_min_m = -1
-interp_molecule_max_m = +1
-                                         /* interp_x      molecule     */
-interp_molecule_positions[0][ 0] =  1    /*   0.00     [0.0, 0.1, 0.2] */
-interp_molecule_positions[0][ 1] =  1    /*   0.04     [0.0, 0.1, 0.2] */
-interp_molecule_positions[0][ 2] =  1    /*   0.06     [0.0, 0.1, 0.2] */
-interp_molecule_positions[0][ 3] =  1    /*   0.10     [0.0, 0.1, 0.2] */
-interp_molecule_positions[0][ 4] =  1    /*   0.14     [0.0, 0.1, 0.2] */
-interp_molecule_positions[0][ 5] =  2    /*   0.16     [0.1, 0.2, 0.3] */
-interp_molecule_positions[0][ 6] =  2    /*   0.20     [0.1, 0.2, 0.3] */
-interp_molecule_positions[0][ 7] =  8    /*   0.80     [0.7, 0.8, 0.9] */
-interp_molecule_positions[0][ 8] =  8    /*   0.84     [0.7, 0.8, 0.9] */
-interp_molecule_positions[0][ 9] =  9    /*   0.86     [0.8, 0.9, 1.0] */
-interp_molecule_positions[0][10] =  9    /*   0.90     [0.8, 0.9, 1.0] */
-interp_molecule_positions[0][11] =  9    /*   0.94     [0.8, 0.9, 1.0] */
-interp_molecule_positions[0][12] =  9    /*   0.96     [0.8, 0.9, 1.0] */
-interp_molecule_positions[0][13] =  9    /*   1.00     [0.8, 0.9, 1.0] */
-\end{verbatim}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Smoothing}
-
-The way the generalized polynomial interpolator is implemented it's
-easy to also do Savitzky-Golay smoothing.%%%
-\footnote{%%%
-	 See section~14.8 of Numerical Recipes
-	 (2nd edition) for a general discussion of
-	 Savitzky-Golay smoothing.
-	 }%%%
-{}  This is best described by way of an example:  Suppose we're doing
-1-D cubic interpolation, using
-(by table~\ref{localinterp/tab-molecule-size+centering})
-4-point molecules.  In other words, at each interpolation point we
-use a cubic interpolation polynomial fitted to 4~surrounding data points.
-For Savitzky-Golay smoothing, we would instead {\em least-squares fit\/}
-a cubic polynomial to some {\em larger\/} number of surrounding data
-points.  This combines interpolation with smoothing, so there's less
-amplification of noise in the input data in the interpolation outputs.
-
-The optional input
-\begin{verbatim}
-const CCTK_INT smoothing;
-\end{verbatim}
-specifies how much (how many points) to enlarge the interpolation
-molecule for this.  The default is 0 (no smoothing).  1 would mean to
-enlarge the molecule by 1 point (\eg{} to use a 5-point molecule instead
-of the usual 4-point one for cubic interpolation).  2 would mean to
-enlarge by 2 points, (\eg{} to use a 6-point molecule for cubic
-interpolation).  Etc etc.
-
-Note that in $>1$~dimension, the maximum-degree Lagrange interpolation
-already uses more data points than the number of free parameters in
-the interpolating polynomials, \ie{} it already does some Savitzky-Golay
-smoothing.  For example, in 2-D a generic cubic polynomial
-$f(x,y) = \sum_{i+j \le 3} c_{ij} x^i y^j$ has 10 free parameters
-$c_{ij}$, which we least-squares fit to the 16 data points in the
-$4 \times 4$ molecule.
-
-Savitzky-Golay smoothing is basically free apart from the increase in
-the molecule size, e.g. a \verb|smoothing|=2 cubic interpolation has
-exactly the same cost as any other 6-point--molecule interpolation.
-
-The current implementation has all the framework for smoothing, but
-only the \verb|smoothing=0| case is implemented at the moment.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Supressing Interpolation}
-
-Sometimes you may want to call the interpolator, but not actually
-do any interpolation for some or all of the arrays:
-\begin{itemize}
-\item	you may be doing a query
-\item	you may be just checking values in the parameter table
-\item	you may want to interpolate varying subsets of some set
-	of arrays
-\item	you may have no interpolation points, yet still want to
-	do the interpolator call for other reasons
-\end{itemize}
-
-In such cases there are several possible ways you can supress some
-or all of the interpolations:
-\begin{description}
-\item[\mathversion{bold}
-      Set $\text{\tt N\_interp\_points} = 0$:]\mbox{}\\
-	This supresses all interpolation.  In this case
-	\verb|interp_coords| may also be \verb|NULL| if desired.
-\item[\mathversion{bold}
-      Set $\text{\tt N\_input\_arrays} = 0$ and
-      $\text{\tt N\_output\_arrays} = 0$:]\mbox{}\\
-	This suppresses all interpolation.  However, note that
-	some parameter-table entries like \verb|operand_indices|
-	and \verb|operation_codes| are specified as being arrays
-	whose size depends on \verb|N_output_arrays|, and it's an
-	error for these table entries to be present with the wrong
-	size (in this case, with any nonzero size).
-\item[\mathversion{bold}
-      Set $\text{\tt output\_arrays[out]} = \text{\tt NULL}$:]\mbox{}\\
-	This supresses the interpolation for this particular
-	output array.  This is probably the cleanest way of
-	selectively turning individual arrays on and off.
-	If $\verb|output_arrays[out]| = \verb|NULL|$ and
-	the input arrays aren't needed for any queries you're doing,
-	then it's useful to also set the corresponding
-	$\verb|input_arrays[in]| = \verb|NULL|$, to supress the
-	interpolator fetching data (which in this case will never be used)
-	from the input arrays.
-\end{description}
-
-Note that the combination
-$\verb|input_arrays[in]| = \verb|NULL|$ and
-$\verb|output_arrays[out]| \ne \verb|NULL|$ is an error: conceptually
-you're asking for an interpolation (the output array pointer is
-non-\verb|NULL|) but not supplying any input data (the input array
-pointer is \verb|NULL|).
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{Implementation}
-
-This interpolator's basic design is to use separate specialized code
-for each combination of
-\begin{verbatim}
-   (N_dims, molecule_family, order, smoothing)
-\end{verbatim}
-\ie{} in practice for each distinct choice of interpolation molecule.
-Maple is used to generate all the interpolation coefficients.
-The C preprocessor is then used to generate all the specialized code
-from a single master ``template'' file.  The template code uses
-\verb|#ifdef|s to handle lower dimensions with no extra overhead,
-\eg{} 1-D/2-D interpolation is basically just as efficient as in
-a purpose-built 1-D/2-D interpolator, with no extra overhead imposed
-by the interpolator also supporting higher-dimensional interpolation.
-
-The Maple code which generates the interpolation coefficients is quite
-general-purpose, and can handle an arbitrary dimensionality and molecule
-size/shape.  Generating new coefficients can be rather time-consuming,
-though, \eg{} the current coefficients for 3-D for orders~1-4
-take about 8~cpu minutes to generate using Maple~7 on a 1.7~GHz~P4.
-
-Note that when compiling the code in the directory
-\verb|src/GeneralizedPolynomial-Uniform| of this thorn, you may get
-compiler warnings about casts discarding \verb|const| qualifiers from
-pointers in (the \verb|#include|-ed file) \verb|template.c|.  Don't
-worry -- the code is actually ok.%%%
-\footnote{%%%
-	 The warnings are really a design bug in C's
-	 const-qualifier rules when pointers point to
-	 other pointers.  See question~11.10 in the
-	 (excellent!) {\tt comp.lang.c} online FAQ
-	 ({\tt http://www.eskimo.com/~scs/C-faq/top.html})
-	 for further discussion.  gcc~2.95.3 gives the
-	 spurious warnings, but gcc~3.2.1 doesn't, nor
-	 does the Intel C compiler in any version I've
-	 tried.  Also, \Cplusplus{} has fixed the problems
-	 in C's const-qualifier rules, so the same code
-	 compiled as \Cplusplus{} shouldn't give any
-	 warnings (and doesn't on the compilers I've tried).
-	 }%%%
-{}  You may also get compiler warnings about unused variables in
-this same file; again don't worry, the code is ok.
-
-See the \verb|README| file in the source code directory
-\verb|LocalInterp/src/UniformCartesian/| for further details on the
-implementation.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{A Simple Example of {\tt CCTK\_InterpLocalUniform} Usage}
-
-Here's a simple example in C, of interpolating a \verb|CCTK_REAL| and a
-\verb|CCTK_COMPLEX| $10 \times 20$ 2-D array, at 5 interpolation points,
-using cubic interpolation.
-
-Note that since C allows arrays to be initialized only if the
-initializer values are compile-time constants, we have to declare the
-\verb|interp_coords[]|, \verb|input_arrays[]|, and \verb|output_arrays[]|
-arrays as non-\verb|const|, and set their values with ordinary (run-time)
-assignment statements.  In \Cplusplus, there's no restriction on initializer
-values, so we could declare the arrays \verb|const| and initialize them
-as part of their declarations.
-
-\begin{verbatim}
-#define N_DIMS   2
-#define N_INTERP_POINTS   5
-#define N_INPUT_ARRAYS    2
-#define N_OUTPUT_ARRAYS   2
-
-/* (x,y) coordinates of data grid points */
-#define X_ORIGIN   ...
-#define X_DELTA    ...
-#define Y_ORIGIN   ...
-#define Y_DELTA    ...
-const CCTK_REAL origin[N_DIMS] = { X_ORIGIN, Y_ORIGIN };
-const CCTK_REAL delta [N_DIMS] = { X_DELTA,  Y_DELTA  };
-
-/* (x,y) coordinates of interpolation points */
-const CCTK_REAL interp_x[N_INTERP_POINTS];
-const CCTK_REAL interp_y[N_INTERP_POINTS];
-const void* interp_coords[N_DIMS];              /* see note above */
-
-/* input arrays */
-/* ... note Cactus uses Fortran storage ordering, i.e. X is contiguous */
-#define NX   10
-#define NY   20
-const CCTK_REAL    input_real   [NY][NX];
-const CCTK_COMPLEX input_complex[NY][NX];
-const CCTK_INT input_array_dims[N_DIMS] = { NX, NY };
-const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
-        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_COMPLEX };
-const void* input_arrays[N_INPUT_ARRAYS];       /* see note above */
-
-/* output arrays */
-CCTK_REAL    output_real   [N_INTERP_POINTS];
-CCTK_COMPLEX output_complex[N_INTERP_POINTS];
-const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
-        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_COMPLEX };
-void* output_arrays[N_OUTPUT_ARRAYS];           /* see note above */
-
-int operator_handle, param_table_handle;
-operator_handle = CCTK_InterpHandle("my interpolation operator");
-if (operator_handle < 0)
-        CCTK_WARN(-1, "can't get interpolation handle!");
-param_table_handle = Util_TableCreateFromString("order=3");
-if (param_table_handle < 0)
-        CCTK_WARN(-1, "can't create parameter table!");
-
-/* initialize the rest of the parameter arrays */
-interp_coords[0] = (const void *) interp_x;
-interp_coords[1] = (const void *) interp_y;
-input_arrays [0] = (const void *) input_real;
-input_arrays [1] = (const void *) input_complex;
-output_arrays[0] = (      void *) output_real;
-output_arrays[1] = (      void *) output_complex;
-
-/* do the actual interpolation, and check for error returns */
-if (CCTK_InterpLocalUniform(N_DIMS,
-                            operator_handle, param_table_handle,
-                            origin, delta,
-                            N_INTERP_POINTS,
-                               CCTK_VARIABLE_REAL,
-                               interp_coords,
-                            N_INPUT_ARRAYS,
-                               input_array_dims,
-                               input_array_type_codes,
-                               input_arrays,
-                            N_OUTPUT_ARRAYS,
-                               output_array_type_codes,
-                               output_arrays) < 0)
-        CCTK_WARN(-1, "error return from interpolator!");
-\end{verbatim}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\subsection{An Example of Interpolating Derivatives}
-\label{localinterp/sect-example-derivatives}
-
-Consider the problem described earlier: the computation of all the
-2nd~derivatives of the 3-metric at a set of interpolation points on
-a 2-sphere.  This example shows how we could do this in C.
-
-[Since we're not using \Cplusplus, we again declare the pointer arrays
-non-\verb|const|, and ``initialize'' them with ordinary (run-time)
-assignment statements.  However, in this example, for greater clarity
-we place these assignment statements right after the declarations.
-Since C allows declarations only at the start of a \verb|{ }| block,
-not in the middle of a block, we nest the rest of the program in extra
-blocks (with the \verb|{ }| indented 2 spaces to distinguish them from
-``normal'' \verb|{ }| pairs) to allow for further declarations.  The
-reader will have to decide for herself whether this style is more or
-less ugly than separating the declarations and initializations of the
-pointer arrays.]
-\begin{verbatim}
-#define N_DIMS        3
-
-/* interpolation points */
-#define N_INTERP_POINTS 1000
-const CCTK_REAL interp_x[N_INTERP_POINTS],
-                interp_y[N_INTERP_POINTS],
-                interp_z[N_INTERP_POINTS];
-const void* interp_coords[N_DIMS];              /* see note above */
-interp_coords[0] = (const void *) interp_x;
-interp_coords[1] = (const void *) interp_y;
-interp_coords[2] = (const void *) interp_z;
-
-  {
-/* dimensions of the data grid */
-#define NX   30
-#define NY   40
-#define NZ   50
-
-/* input arrays */
-/* n.b. we use Fortran storage order: X is contiguous, Z least so */
-#define N_INPUT_ARRAYS  6
-const CCTK_REAL gxx[NZ][NY][NX], gxy[NZ][NY][NX], gxz[NZ][NY][NX],
-                                 gyy[NZ][NY][NX], gyz[NZ][NY][NX],
-                                                  gzz[NZ][NY][NX];
-
-const CCTK_INT input_array_dims[N_DIMS] = {NX, NY, NZ};
-const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
-        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                                CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                                                    CCTK_VARIABLE_REAL };
-const void* input_arrays[N_INPUT_ARRAYS];       /* see note above */
-input_arrays[0] = (const void *) gxx;
-input_arrays[1] = (const void *) gxy;
-input_arrays[2] = (const void *) gxz;
-input_arrays[3] = (const void *) gyy;
-input_arrays[4] = (const void *) gyz;
-input_arrays[5] = (const void *) gzz;
-
-  {
-/* output arrays */
-#define N_OUTPUT_ARRAYS 36
-CCTK_REAL
- dxx_gxx[N_INTERP_POINTS], dxy_gxx[N_INTERP_POINTS], dxz_gxx[N_INTERP_POINTS],
-                           dyy_gxx[N_INTERP_POINTS], dyz_gxx[N_INTERP_POINTS],
-                                                     dzz_gxx[N_INTERP_POINTS],
- dxx_gxy[N_INTERP_POINTS], dxy_gxy[N_INTERP_POINTS], dxz_gxy[N_INTERP_POINTS],
-                           dyy_gxy[N_INTERP_POINTS], dyz_gxy[N_INTERP_POINTS],
-                                                     dzz_gxy[N_INTERP_POINTS],
- dxx_gxz[N_INTERP_POINTS], dxy_gxz[N_INTERP_POINTS], dxz_gxz[N_INTERP_POINTS],
-                           dyy_gxz[N_INTERP_POINTS], dyz_gxz[N_INTERP_POINTS],
-                                                     dzz_gxz[N_INTERP_POINTS],
- dxx_gyy[N_INTERP_POINTS], dxy_gyy[N_INTERP_POINTS], dxz_gyy[N_INTERP_POINTS],
-                           dyy_gyy[N_INTERP_POINTS], dyz_gyy[N_INTERP_POINTS],
-                                                     dzz_gyy[N_INTERP_POINTS],
- dxx_gyz[N_INTERP_POINTS], dxy_gyz[N_INTERP_POINTS], dxz_gyz[N_INTERP_POINTS],
-                           dyy_gyz[N_INTERP_POINTS], dyz_gyz[N_INTERP_POINTS],
-                                                     dzz_gyz[N_INTERP_POINTS],
- dxx_gzz[N_INTERP_POINTS], dxy_gzz[N_INTERP_POINTS], dxz_gzz[N_INTERP_POINTS],
-                           dyy_gzz[N_INTERP_POINTS], dyz_gzz[N_INTERP_POINTS],
-                                                     dzz_gzz[N_INTERP_POINTS];
-const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
-   = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                                               CCTK_VARIABLE_REAL,
-       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                                               CCTK_VARIABLE_REAL,
-       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                                               CCTK_VARIABLE_REAL,
-       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                                               CCTK_VARIABLE_REAL,
-       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                                               CCTK_VARIABLE_REAL,
-       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-                                               CCTK_VARIABLE_REAL };
-void* output_arrays[N_OUTPUT_ARRAYS];           /* see note above */
-output_arrays[ 0] =      (void *) dxx_gxx;
-output_arrays[ 1]  =     (void *) dxy_gxx;
-output_arrays[ 2]   =    (void *) dxz_gxx;
-output_arrays[ 3]    =   (void *) dyy_gxx;
-output_arrays[ 4]     =  (void *) dyz_gxx;
-output_arrays[ 5]      = (void *) dzz_gxx;
-output_arrays[ 6] =      (void *) dxx_gxy;
-output_arrays[ 7]  =     (void *) dxy_gxy;
-output_arrays[ 8]   =    (void *) dxz_gxy;
-output_arrays[ 9]    =   (void *) dyy_gxy;
-output_arrays[10]     =  (void *) dyz_gxy;
-output_arrays[11]      = (void *) dzz_gxy;
-output_arrays[12] =      (void *) dxx_gxz;
-output_arrays[13]  =     (void *) dxy_gxz;
-output_arrays[14]   =    (void *) dxz_gxz;
-output_arrays[15]    =   (void *) dyy_gxz;
-output_arrays[16]     =  (void *) dyz_gxz;
-output_arrays[17]      = (void *) dzz_gxz;
-output_arrays[18] =      (void *) dxx_gyy;
-output_arrays[19]  =     (void *) dxy_gyy;
-output_arrays[20]   =    (void *) dxz_gyy;
-output_arrays[21]    =   (void *) dyy_gyy;
-output_arrays[22]     =  (void *) dyz_gyy;
-output_arrays[23]      = (void *) dzz_gyy;
-output_arrays[24] =      (void *) dxx_gyz;
-output_arrays[25]  =     (void *) dxy_gyz;
-output_arrays[26]   =    (void *) dxz_gyz;
-output_arrays[27]    =   (void *) dyy_gyz;
-output_arrays[28]     =  (void *) dyz_gyz;
-output_arrays[29]      = (void *) dzz_gyz;
-output_arrays[30] =      (void *) dxx_gzz;
-output_arrays[31]  =     (void *) dxy_gzz;
-output_arrays[32]   =    (void *) dxz_gzz;
-output_arrays[33]    =   (void *) dyy_gzz;
-output_arrays[34]     =  (void *) dyz_gzz;
-output_arrays[35]      = (void *) dzz_gzz;
-
-  {
-/* integer codes to specify the derivatives */
-/* (for best efficiency we group all operations on a given input together) */
-const CCTK_INT operand_indices[N_OUTPUT_ARRAYS]
-   = { 0, 0, 0, 0, 0, 0,
-       1, 1, 1, 1, 1, 1,
-       2, 2, 2, 2, 2, 2,
-       3, 3, 3, 3, 3, 3,
-       4, 4, 4, 4, 4, 4,
-       5, 5, 5, 5, 5, 5 };
-#define DERIV(x)   x
-const CCTK_INT operation_codes[N_OUTPUT_ARRAYS]
-   = { DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
-       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
-       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
-       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
-       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
-       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33) };
-
-int operator_handle, param_table_handle;
-operator_handle = CCTK_InterpHandle("my interpolation operator");
-if (operator_handle < 0)
-        CCTK_WARN(-1, "can't get interpolation handle!");
-
-param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
-if (param_table_handle < 0)
-        CCTK_WARN(-1, "can't create parameter table!");
-if (Util_TableSetInt(param_table_handle, 3, "order") < 0)
-        CCTK_WARN(-1, "can't set order in parameter table!");
-if (Util_TableSetIntArray(param_table_handle,
-                          N_OUTPUT_ARRAYS, operand_indices,
-                          "operand_indices") < 0)
-        CCTK_WARN(-1, "can't set operand_indices array in parameter table!");
-if (Util_TableSetIntArray(param_table_handle,
-                          N_OUTPUT_ARRAYS, operation_codes,
-                          "operation_codes") < 0)
-        CCTK_WARN(-1, "can't set operation_codes array in parameter table!");
-
-if (CCTK_InterpLocalUniform(N_DIMS,
-                            operator_handle, param_table_handle,
-                            origin, delta,
-                            N_INTERP_POINTS,
-                               CCTK_VARIABLE_REAL,
-                               interp_coords,
-                            N_INPUT_ARRAYS,
-                               input_array_dims,
-                               input_array_type_codes,
-                               input_arrays,
-                            N_OUTPUT_ARRAYS,
-                               output_array_type_codes,
-                               output_arrays) < 0)
-        CCTK_WARN(-1, "error return from interpolator!");
-  }
-  }
-  }
-\end{verbatim}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\section{Acknowledgments}
-
-Thanks to Thomas Radke for contributing the former \thorn{PUGHInterp}
-\verb|CCTK_InterpLocal()| interpolator,
-to Ian Hawke and Erik Schnetter for helpful comments on the documentation,
-to Tom Goodale and Thomas Radke for many useful design discussions,
-to Erik Schnetter for bug reports,
-and to all the Cactus crew for a great infrastructure!
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 % Do not delete next line
 % END CACTUS THORNGUIDE