Interpolation in Patches
Most of this will be a pretty straightforward extension to triangular coordinates of the material in
Interpolation on Edges. First, however, comes a short description of where the
N N N N N N
/ \ / \ / \ / \ / \ / \
V---V---V---V---V---V---V
\ / \ / \ / \ / \ / \ / \
V---V---V---V---V---V---V
\ / \ / \ / \ / \ / \ /
S S S S S S
patches are in the global grid, which is
a latitude-vs-longitude map of the sphere, composed of twenty-four
triangular regions (faces) that are connected as depicted to the right.
M=6:
V
• •
• • •
• • • •
• • • • •
• • • • • •
W • • • • • E
• • • • • •
• • • • •
• • • •
• • •
• •
V
The map can be folded and glued together to make a 'globe' in the
form of a hexagonal antiprism.
The faces occur in pairs that fill a diamond by sharing a base that is a line of constant
latitude, as depicted to the left. There are twelve 'equatorial' triangles with bases that are
alternately at 30º and -30º, and their partners have vertices at the corresponding 'north' and 'south'
poles. The vertices marked 'V' at the left might be 'N', 'V' or 'S' at the right, and the ones marked 'E'
or 'W' are always 'V'. Each face has M2 triangular patches that are defined
by integer-values of triangular coordinates (E,W,V) for which E + W + V = M , W = M
at 'W' , E = M at 'E'
and V = M at 'V'. The fractional parts of E , W and V can be used to define
triangular coordinates, e, w and v, with e + w + v = 1 in the patches. All of that
is to be done by the algorithm that requires interpolated data, and all the interpolator
sees is arguments (e,w,P), where e and w are just numbers, and P is typically a pointer
to a structure (patch) that contains everything else it will need.
The north- and south-facing patches are mirror images of one another, with three interior vertices and
nine exterior vertices
whose names are assigned as in the diagrams. The names aren't really points of the compass;
the first letter of two
nw ne
\ /
\ /
wn-----n-----en
\ / \ /
\ / \ /
ww-----w-----e-----ee
/ \ / \
/ \ / \
ws s es
wn n en
\ / \ /
\ / \ /
ww-----w-----e-----ee
/ \ / \
/ \ / \
ws-----s-----es
/ \
/ \
sw se
is meant to convey a relative position, as in 'nw' for north of west and 'es' for east of south, etc.
Because of the symmetry, the same interpolator works for either case, provided the data are specified in the
appropriate order, e.g.
P → 'n'w'e'ww'ee'es'nw'ne'ws's'en'wn'
or
P → 's'w'e'ww'ee'en'sw'se'wn'n'es'ws'
Thus, the algorithms for interpolation in patches can all be discussed in the context
of a north-facing patch with triangular coordinates, e , w and n.
Linear interpolation
Since e, w and n are linear functions of longitude and latitude (or x and y on the map
or in the tangent plane),
F1(e,w,n) = eFe + wFw + nFn
= Fn + (Fe - Fn)e + (Fw - Fn)w(1)
is the linear interpolation of F(x,y) that is Fw at 'w', Fe at 'e'
and Fn at 'n'. As in one dimension, the error is systematic, and it causes an unphysical
degradation of the steepness of shock-waves and other moving fronts. Corrections of it will be made up of
cubic polynomial functions of e, w and n that are zero at the interior vertices, 'e', 'w' and 'n'.
Cubic interpolation
eew:
0 0
\ /
1 0 -1
/ \
2 0 0 -4
/ \
0 1 0
eww:
0 0
\ /
-1 0 1
/ \
-4 0 0 2
/ \
0 1 0
ew(1+w):
0 0
\ / \
-2 0 0
/ \ \
-6 0 0 0
/ \
0 2 0
ew+ne+wn:
-2 -2
-1 0 -1
-2 0 0 -2
-2 -1 -2
ewn:
0 0
\ /
-1 0 -1
/ \
0---0---0---0
/ \
0 -1 0
ewn(n-1):
0 0
\ /
0---0---0
/ \
0---0---0---0
/ \
0 2 0
e:
-1 0
/ /
-1 0 1
/ / /
-1 0 1 2
/ / /
0 1 2
w:
0 -1
\ \
1 0 -1
\ \ \
2 1 0 -1
\ \ \
2 1 0
ew:
0 0
\ /
-1 0 -1
/ \
-2 0 0 -2
/ \
0 1 0
Shown above and to the right are diagrams of the values of some of the
terms in the twelve Lagrange polynomials. Consider first the last diagram,
from which it follows that Ls = ewn(n-1)/2 and its two rotations are
three of the twelve. By inspection they are quartic polynomials, and it follows
easily from the diagrams that the sum of the three is the cubic function, -ewn.
From the diagram for ew(1+w) it follows easily that Lww = -ew(1+w+(e-n)n)/6 is the Lagrange
polynomial that is '1' at the vertex 'ww'. After rotations and reflections of Lww
there remains only Ln and two rotations of it. Again from the diagrams,
ew+ne+wn+ewn+2 is '2' at the interior nodes and '0' at the exterior nodes, and from
the first diagram it follows that
Le = (e/2)(ew+ne+wn+ewn+2). Thus the polynomial of lowest degree that
has the prescribed values of F at the twelve vertices is
FL(e,w,n) = LnFn + LsFs +
LeeFee + LwwFww
+ LwFw + LenFen +
LnwFnw + LesFes
+ LeFe + LwnFwn +
LwsFws + LneFne,(2)
and every last one of the L's is a quartic function of e, w and n.
Just as in one dimensional interpolation, the use of Lagrange polynomials
introduces nonphysical waves in semi-Lagrangian algorithms for transport equations, and
more reliable results are obtained with an approximate Hermite interpolation where six of
the values of F at external vertices are used to estimate three pairs of slopes of F at the internal
vertices. The six slopes are also aligned pairwise on edges that span four vertices, so
results from cubic interpolation on edges can be used, but there are some modifications.
On the triangular patches the
vertices appear as three rotations of six vertices where F - F1
has the values indicated below:
0
/ \
/ \
/ \
Fww+Fe-2Fw---0-------0--- Fee+Fw-2Fe
\ /
\ /
\ /
Fs+Fn-Fw-Fe
On the 'ew'-edge, where n = 0, the result (edges-eq(8)) is
F3(e,w,0) = F1(e,w,0) +
a(F1-F)wweww +
a(F1-F)eeeew ,(3)
where a is an adjustable parameter, generally taken to be in the range,
0 <- a <- 1.
The problem with eq(3) is that its contribution on the
patch is
(F3-F)ew =
a(Fe-Fee+Fw-Fww)ew(1-n)/2 +
a(3Fw-Fww+Fee-3Fe)ew(w-e)/2 ,(4)
Now the second term and its rotations provide the three antisymmetric cubic terms that appear in
Lagrange polynomials, but the symmetric quadratic terms are not isolated. If only the three
rotations of eq(3) are employed, the resulting approximate interpolation cannot possibly
get all cases where F is a quadratic function right. The quick fix,
F3 = F1 +
a(F1-F)wwew(w+n/2) +
a(F1-F)eeew(e+n/2)
+
a(F1-F)esne(e+w/2) +
a(F1-F)nwne(n+w/2)
+
a(F1-F)newn(n+e/2) +
a(F1-F)wswn(w+e/2) ,(5)
eliminates the troublesome ewn's , and if F is quadratic and a=1/2, the central difference
approximations in edges-eq(7) are exact, the cubic terms in rotations
of eq(4) are zero, and F3 in eq(5) is F(e,w,n).
ew+ne+wn slopes:
1 1
\ /
n
1/2 / \ 1/2
c c
/ \
1---w---c---e---1
/ \
1 1/2 1
ewn slopes:
0 0
\ /
n
1/4 / \ 1/4
c c
/ \
0---w---c---e---0
/ \
0 1/4 0
(ew + ne + wn)A6 + 4ewn(A3 - A6/2) ,
where A6 is the average of the six coefficients in eq(5) and
A3 = a(Fn - Fs +
Fw - Fen + Fe - Fwn)/3.(6)
At the point, e = w = n = 1/3, where ewn/(ew+ne+wn)
is largest, the totally symmetric part is
A6/3 + (4/27)(A3 - A6/2) =
(A6/3)(1 + 4(A3/A6 - 1/2)/9).
Roughly speaking, the totally symmetric quadratic term is the dominant one if
|A3/A6 - 1/2| < 1/2
or
0 < A3/A6 < 1 .
If ewn is discarded the approximation has nine degrees of freedom (values of F), and
in eqs(1 and 5) the computational effort is 8'x' 22'+' per function. The coefficients
of the nine values of F can be collected to evaluate the basis, with its computational
effort of 9'x' 8'+'. Which of the two is the more efficient may well depend upon the
computational environment.
Babylonian Approximations
What remains now is to examine a few approximations that are technically less accurate than
cubic interpolation; they will be tested to see if they are adequate corrections of linear
interpolation. The criterion that sets the order in which they are presented is the minimum
number of multiplications per function (hence, the title of this section). Multiplication or
division by a power of two is not counted at all, and multiplications by other integers are
counted as additions.
Consider first the linear interpolation (three ways)
F1n(e,w) = Fn + (Fe - Fn)e +
(Fw - Fn)w ,
F1w(n,e) = Fw + (Fn - Fw)n +
(Fe - Fw)e ,
F1e(w,n) = Fe + (Fw - Fe)w +
(Fn - Fe)n ,(7)
and the three linear extrapolations,
Fxn(e,w) = Fn + eF'nw +
wF'ne ,
Fxw(n,e) = Fw + nF'ws +
eF'ww ,
Fxe(w,n) = Fe + wF'ee +
nF'es .(8)
The names of inward slopes are consistent with those of Lagrange polynomials, and the estimates
of the slopes of F (not F-F1) are rotations and reflections of
F'ww = Fw - Fww (extrapolation)
or F'ww = (Fe - Fww)/2
(central difference). The two choices correspond with a = 1 or a = 1/2,
and neither counts as a multiplication.
For superlinear (a = 1/2) and ultralinear (a = 1)
interpolation it
doesn't matter which form of the linear interpolarion is used, but the use of the three
helps to clarify the tests
that determine whether or not the extrapolations define a convex hull that
limits deviations from F1. The tests are all reflections and rotations
of the question,
is (F'ww + Fw - Fe)
> 0 or < 0 or = 0 ?
The argument that accompanies all this is: if the differences between exterior and interior slopes
(F'xv - F'1v)
are not all positive
or all negative, F has at least one inflection point on an edge of the patch, and there is
no convex hull to limit estimates based upon the extrapolations. Whether it be argued that
F1 is not systematically wrong near inflection points or that extrapolation is
systematically
wrong there, the choice when both tests fail is to use any one of the versions of eq(7). All told,
the cost of that case is 2'x' 14'+', and that would be a bargain indeed if inflection points
were less rarely found in reasonably well-behaved functions.
In the more common cases: if all
(F'xv - F'1v) > 0 the
Babylonian approximation is the least of the three extrapolations in eq(8), and if all
(F'xv - F'1v) < 0 it
is the greatest of the three. Either way, the cost is 6'x' 18'+', and it saves two multiplies
and a few adds because there is no need to evaluate F1.
Quadratic approximations
A quadratic approximation of eq(5) (c.f. eq(4)) is
F2 = Fn + (Fe - Fn)e +
(Fw - Fn)w +
(aew/2)(Fe + Fw - Fee - Fww)
+
(ane/2)(Fn + Fe - Fnw - Fes) +
(awn/2)(Fw + Fn - Fws - Fne) .(9)
This, with a cost of 5'x' 16'+' per function, is slightly more parsimonious than the previous two
approximations. It still has nine degrees of freedom, but collecting the coefficients of the
values of F would be a bad idea.
A bit further on the way to downright stingy is the average of the corrections in eq(5),
F2A = Fn + (Fe - Fn)e +
(Fw - Fn)w
+
(a/6)(ew + ne + wn)(2(Fw + Fe + Fn) -
Fww - Fee - Fes - Fnw -
Fne - Fws) .(10)
It would be very nice indeed if this 3'x' 13'+' approximation were to help, even a little, to
compensate for numerical dissipation, but wait __ there is one more.
None of the above
Consider
F1.5 = Fn + (Fe - Fn)e +
(Fw - Fn)w
+
(a/6)(b)(2(Fw + Fe + Fn) -
Fww - Fee - Fes - Fnw -
Fne - Fws) ,(11)
where b is a representative value of ew + ne + wn . The maximum value of the symmetric
quadratic correction is 1/3 at e = w = n = 1/3, and
a <- 1, so it is reasonable to try a value of ab
that is somewhat less than 1/3. While something like ab = 1/4 might work, either of 3/16 or 3/32
would cancel the 3 in ab/6, and the result would be a 2'x' 13'+' approximation. It's hard
to believe there could be a stingier correction that had much of anything to do with convexity.