Refinement implementation for simplices
[Refinement]

This is mainly based on Jürgen Beys <http://www.igpm.rwth-aachen.de/bey/> dissertation. The relevant part is available from <http://www.igpm.rwth-aachen.de/Download/reports/bey/simplex.ps.gz>.

Terminology

Kuhn simplex: To triangulate hypercubes we use the Kuhn triangulation. The members of this triangulation we call Kuhn simplices. The Kuhn simplices are indexed by their corresponding permutation.
Kuhn0 simplex: The Kuhn simplex corresponding to the permutation number 0.
size of a Kuhn simplex: The size of a kuhn simplex is equal to the size of the hypercube that it triangulates.
width of a Kuhn simplex: See size of a Kuhn simplex.

Describing Kuhn simplices by their permutation

A Kuhn simplex of dimension n can be described by its size s and a permutation of the vector $\vec{p}=(0,\ldots,n-1)$ . To get the coordinates of the corners $\vec{x}_0,\ldots,\vec{x}_n$ of the simplex, you do as follows:

Start at the origin $\vec{x}_0$ .
For each dimension d from 0 to n-1:
- $\vec{x}_{d+1}:=\vec{x}_d+s\cdot\vec{e}_{p_d}$ ( $\vec{e}_i$ is the unit vector in direction i.)

Number of vertices in a Kuhn0 simplex

Let N(n, x) be the number of gridpoints within an n-dimensional Kuhn0 simplex of x gridunits width.

The number of points in a 0-dimensional simplex is 1, independent of its width.

N(0, x) = 1

The recursion formula is

$N(n+1,x)=\sum^x_{i=0}N(n,i)$

We slice the n+1 dimensional simplex ortogonal to one of the dimensions and sum the number of points in the n dimensional subsimplices.

This formula is satisfied by the binomial koefficient

$N(n,x)=\left({n+x}\atop x\right)$

See Bronstein, Semendjajew, Musiol, Mühlig "Taschenbuch der Mathematik" (1999), Formula (1.35c)

Observations:

N(n, 0) = 1
N(n, x) = N(x, n)

Index of a vertex within a Kuhn0 simplex

The image shows the Kuhn0 tetrahedron of width 2 (wireframe). It is partitioned into a tetrahedron (green), a triangle (red), a line (blue), and a vertex (black), each of width 1 and each a Kuhn0 simplex.

Let us calculate the index of vertex 9, which has the coordinates $(x_0,x_1,x_2)=(2,2,2)$

First we count the number of vertices in the green tetrahedron of width (4). Then we take the green tetrahedron away. Whats left is a triangle, which extends into the (1,2)-plane.
Now we count the number of vertices in the red triangle of width (3). Again we take the counted points away and are left with a line which extends into direction 2.
We take the blue line of width , count the number of vertices (2), and throw away the counted stuff. The only thing remaining is a point, so we're done.
We add the counted stuff together and get indeed 9.

On to a more complicated example: vertex 6 with coordinates $(x_0,x_1,x_2)=(2,1,1)$ .

First count the vertices in the green tetrahedron again (width ). The result is 4.
Count the vertices in the triangle of width (vertex 4), which is just a point. The result is 1.
Count the vertices in the line of width (vertex 5), which is also just a point. The result is 1.
Add everything together and get 6.

The general algorithm for n dimensions and a vertex with coordinates $\vec{x}=(x_0,\ldots,x_{n-1})$ is as follows:

For each dimension d from 0 to n-1
- Count the vertices in the n-d dimensional simplex of width ,
Add all counts together to get the index of the vertex

In formulas it looks like this:

Let $I(n,\vec{x})$ be the index of point $\vec{x}$ in the n-dimensional Kuhn0 simplex. The coordinates measure the position of the point in gridunits and thus are integer.

$I(n,\vec{x})=\sum_{i=0}^{n-1}N(n-i,x_i-1)=\sum_{i=0}^{n-1}\left({n-i+x_i-1}\atop{n-i}\right)$

Since the coordinates of a vertex within the Kuhn0 obey the relation $x_0\geq x_1\geq\ldots\geq x_{n-1}$ , they cannot simply be swapped so the sum is somewhat ugly.

Index of a subelement within a Kuhn0 simplex

We don't know of a way to simply map a subelement of a Kuhn0 simplex to an index number. Luckily, the iterator interface only requires that we be able to calculate the next subelement.

Each subelement is a Kuhn (sub)simplex which triangulates a hypercube. We need to remember the vertex which is the origin of that hypercube and the index of the permutation of that identifies the Kuhn subsimplex. Now to get to the next subelement, we simply need to increment the permutation index, and if the overflows we reset it and increment the origin to the next vertex (we already know how to do that).

Some subelements generated this way are outside the refined Kuhn0 simplex, so we have to check for that, and skip them.

Index of a permutation

[NOTE: There may be some interesting stuff in http://en.wikipedia.org/wiki/Factoradic . I was not aware of it while writing this code, however.]

We need to index the n! permutations of the integers from 0 to n-1 and a way to calculate the permutation if given the index.

I'll discuss the permutation P, which operates on a vector $\vec{x}=(x_0,\ldots,x_{n-1})$ .

P can be made up of n transpositions, $P=T_0\cdots T_{n-1}$ . Each transposition $T_i$ exchanges some arbitrary element $x_{t_i}$ with the element $x_i$ , where $t_i\leq i$ . So we can totally describe $T_i$ by the integer $t_i$ . Thus we can describe P by the integer vector $\vec{t}=(t_0,\cdots,t_{n-1})$ , where $0\leq t_i\leq i$ .

Now we need to encode the vector $\vec{t}$ into a single number. To do that we view $t_i$ as digit i of a number p written in a base faculty notation:

$p=\sum_{i=1}^{n-1}i!t_i$

This number p is unique for each possible permutation P so we could use this as index. There is a problem though: we would like the identity permutation $\vec{x}=P\vec{x}$ to have index 0. So we define the index I of the permutation slightly differently:

$I=\sum_{i=1}^{n-1}i!(i-t_i)$

I can easily calculate the $t_i$ from I:

$i-t_i=(I/i!)\%(i+1)$

('/' is integer division and '' calculates the remainder).

Mapping between some Kuhn and the reference simplex

Step 1 moves each point by its x2 value into x1 direction. Step 2 moves each point by its new x1 value into x0 direction.

The algorithm to transform a point $\vec{x}=(x_0,\ldots,x_{n-1})$ from the reference simplex of dimension n to the Kuhn0 simplex (as illustrated in the image) is as follows:

For each dimension d from n-2 down to 0:
- $x_d:=x_d+x_{d+1}$ .

The reverse (from Kuhn0 to reference) is simple as well:

For each dimension d from 0 up to n-2:
- $x_d:=x_d-x_{d+1}$ .

Arbitrary Kuhn simplices

For arbitrary Kuhn simplices we have to take the permutation of that simplex into account. So to map from the reference simplex of n dimensions to the Kuhn simplex with the permutation P (which is described by the vector $\vec{p}=P(0,\ldots,n-1)$ ) we do:

For each dimension d from n-2 down to 0:
- $x_{p_d}:=x_{p_d}+x_{p_{d+1}}$ .

And or the reverse:

For each dimension d from 0 up to n-2:
- $x_{p_d}:=x_{p_d}-x_{p_{d+1}}$ .

Refinement implementation for simplices [Refinement]