This is the README file for the main directory of the Melkes package
of interpolation routines, so named after ref. [0].

The Melkes routines may be used for polynomial interpolation based on
data given on the vertices of a box in an arbitrary number of
dimensions.  The data may be function values only, or function values
and gradients, or function values, gradients and Hessians.  Routines
are supplied to evaluate only the interpolant, or the interpolant and
its gradient, or the interpolant, its gradient and its Hessian.

The interpolating polynomial for Melkes interpolation may be described
as a reduced Hermite polynomial --- reduced by comparison with the
full tensor product form.  A precise mathematical description follows;
note that this describes the extension to an arbitrary number of
dimensions of what is called "Type-II" interpolation in Melkes's
original paper [0].

Let p be a nonnegative integer, which will characterize the order of
the interpolation, and let d denote the number of dimensions.
Consider a box in R^d.  The required data for Melkes interpolation
with parameter p are function values and derivatives of total order
less than p, all given at each vertex of the box.  Let x stand for a
generic point in R^d, and let m be a d-vector of nonnegative integers.
Then the monomial x^m is included in the expansion of the
interpolating polynomial if and only if the sum of floor(m(k)/2),
where k ranges over the d dimension indices, is less than p.  Here,
floor(t0) means t0 rounded down to integer.

Thus, Melkes interpolation with parameter p includes in the expansion
all monomials of total order less than 2*p, and has 2*p order of
accuracy for smooth functions.  A q-th order derivative of the
interpolant will provide 2*p-q order of accuracy.  It may be verified
that for each case (characterized by p and d) the number of data items
(function and derivative values at vertices) is indeed equal to the
number of monomials in the expansion.

Let us work out the lowest-order cases in complete detail.

The case p=0 is trivial: there are no data, and there are no monomials
in the expansion.  The interpolating function is identically zero,
which is zeroth order accurate.  (We did not implement this case in
the package).

The case p=1 corresponds to multilinear interpolation.  It requires
function values only, and the terms in the interpolating polynomial
are characterized by having degree equal to 0 or 1 in each variable.
This interpolation is second order accurate for function values,
first order for gradients, and zeroth order accurate (that is, not
convergent) for second derivatives.

The case p=2 may be the most interesting one in practice.  As data it
requires function values and gradients at all vertices.  The monomials
in the interpolating polynomial are characterized by having degree
less than or equal to 3 in every variable and degree exceeding 1 in at
most one variable.  The interpolation will provide 4th order accuracy
for function values.

The case p=3 is the highest order implemented in the present package.
It requires function values, gradients and Hessians at all vertices.
In addition to the monomials included in the expansion for the case
p=2, now are also included monomials having degree 2 or 3 in two
variables and degree 0 or 1 in all others, and those having degree 4
or 5 in one variable and degree 0 or 1 in all others.

It is of interest to compare the case p=2 with tensor product splines.
Both methods give fourth order of accuracy.  On each elementary box,
the local polynomial for tensor product spline interpolation is
characterized by 4^d coefficients whereas the Melkes polynomial is
characterized by (d+1)*2^d coefficients.  The Melkes polynomial is
defined by strictly local data, and so is less tied to a Cartesian
product mesh than is spline interpolation.  Especially, Melkes
interpolation will be ideally suited for use with local Cartesian mesh
refinement.  However, Melkes interpolation with p=2 has with respect
to spline interpolation the disadvantage of requiring also first
derivative data at verices.  In some cases these might be obtained,
for example, by use of one-dimensional spline interpolations.

A manuscript or paper describing the present implementation of Melkes
interpolation will follow in 2002.  In the mean time, prospective
users of the package are asked to contact the author, B. J. Braams, at
braams@courant.nyu.edu, and also to make do primarily with the
internal routine documentations.

Two additional files in this directory provide some help.  The file
INSTALL gives rudimentary installation instructions.  I need further
experience with multiple platforms, and some outside advice, before
finalizing these instructions.  The file USERDOC provides basic
instructions for use of the package.

[0] Frantisek Melkes: Reduced Piecewise Bivariate Hermite
Interpolations.  Numer. Math. 19 (1972) 326--340.