One of the major problems with unstructured meshes is the storage overhead.
In the case of the octree, this amounts to
having to store both the octree and the lookup table. In order to
reduce the storage costs associated with the octree, we compress both
the octree and the lookup table. The sorted, unique, linear
octree can be easily compressed by retaining only the offset of the first
element and the level of subsequent octants. Storing the offset for
each octant requires a storage of three integers (12 bytes) and
a byte for storing the level. Storing only the level represents
a 12x compression as opposed to storing the offset for every
octant.
In addition we need the neighborhood data to be prepared for each element/node. In case of an elemental iteration, we need to store the information pertaining to the 7 additional nodes (in 3D) defining the current hexahedral element. In the case of a node-based iteration, we will need to store the information about 26 neighboring nodes, connected to this node via adjacent elements. These numbers are 3 and 8 in 2D and are illustrated in the following figure.
Indexing required for node-based iterations (left) and elements based iterations (right). Red signifies the current node/element, and the neighboring nodes are represented by blue dots.
We now focus on compressing the element to node mapping, which requires eight integers for each element. In order to devise the best compression scheme, we first estimate the distribution of the node indices. The following lemmas help us analyze the distribution of the indices of the nodes of a given element.
Lemma : The Morton ids of the nodes of a given element are greater than or equal to the Morton id of the element.
Proof : Let the anchor of the given element be $( x, y, z )$ and let it's size be $h$. In that case the anchors of the 8 nodes of the element are given by, $( x, y, z ), ( x + h, y, z ), ( x, y + h, z ), ( x + h, y + h, z ) \cdots$ . By the definition of the Morton ids all of these except $( x, y, z )$ are greater than the Morton id of the element. The node at $( x, y, z )$ is equal to the Morton id of the element. □
It is important to note that because of the direct mapping between nodes and elements by our scheme, it is only necessary to store the indices of the 7 additional nodes. In the following sections we refer to these nodes when we refer to the 7 nodes of a given element.
Corollary 1 : Given a sorted list of Morton ids corresponding to the combined list of elements and nodes of a balanced linear octree, the indices of the 8 nodes of a given element in this list are strictly greater than the index of the element. Moreover, if the nodes are listed in the Morton order, the list of indices is monotonically increasing. If we store offsets in the sorted list, then these offsets are strictly positive.
Corollary 2 : If the indices are stored as offsets from the index of the given element, then the offsets are not sensitive to the distribution of octants whose Morton id is less than that of the given element.
Based on these observations we can estimate the expected range of offsets. Let us assume a certain balanced octree, $O$, with $n$ octants (elements and hanging-nodes) and with maximum possible depth $D_{\max}$. Consider an element in the octree, $o_i$, whose index is $i$, $0 \leqslant i < n$. The offset of the anchor of this element is either $i$ (if the anchor is not hanging) or $n_0 < i$. The indices for the remaining 7 nodes do not depend on octants with index less than $i$. In addition since the indices of the 7 nodes are monotonically increasing, we can store offsets between two consecutive nodes. That is, if the indices of the 8 nodes of an element, $o_i$, are [ ( n_0, n_1, n_2, n_3, n_4, n_5, n_6, n_7 ), ] we only need to store, [ ( n_0-i, n_1 - n_0, n_2 -n_1, n_3 - n_2, n_4 - n_3, n_5 - n_4, n_6 - n_5, n_7 - n_6 ). ]
To efficiently store these offsets, we need to estimate how large these offsets can be. We start with a regular grid, i.e., a balanced octree with all octants at $D_{\max}$. Note that any octree that can be generated at the same value of $D_{\max}$ can be obtained by applying a series of local coarsening operations to the regular grid. Since we only store the offsets it is sufficient to analyze the distribution of the offset values for one given direction, say for a neighbor along the $x$-axis. The expression for all other directions are similar.
For $D_{\max} = 0$, there is only one octant and correspondingly the offset is $1$. If we introduce a split in the root octant, $D_{\max}$ becomes 1, the offset increases by 2 for one octant. This corresponds to the boundary introduced by the previous split. On introducing further splits, the offset is going to increase for those octants that lie on the boundaries of the original splits, and the general expression for the maximum offset can be written as [ \mbox{offset} = 1 + \sum_{i = 1}^{D_{\max}} 2^{ d\cdot i - 1}, ] for a $d$-tree. In addition, a number of other large offsets are produced for intermediate split boundaries. Specifically for a regular grid at maximum depth $D_{\max}$, we shall have $2^{d\cdot( D_{\max} - x )}$ octants with an offset of $1 + \sum^x_{i = 1} 2^{d\cdot i - 1}$ . As can be clearly seen from the expression, the distribution of the offsets is geometric. With the largest number of octants having small offsets.
Coming back to the case of general balanced octrees, we observe that any of these can be obtained from a regular grid by a number of coarsening operations. The only concern is if the coarsening can increase the offset for a given octant. The coarsening does not affect octants that are greater than the current octant (in the Morton order). For those which are smaller, the effect is minimal since every coarsening operation reduces the offsets that need to be stored.
Golomb coding is a form of entropy encoding that is
optimal for geometric distributions, that is, when small values are
vastly more common than large values. Since, the distribution of the
offsets is geometric, we expect a lot of offsets with small
values and fewer occurrences of large offsets. The Golomb coding uses
a tunable parameter $M$ to divide an input value into two parts:
$q$, the result of a division by $M$, and $r$, the remainder.
In our implementation, the remainder is stored as a
byte, and the quotient as a short. On an average, we observe
one large jump in the node indices, and therefore the amortized
cost of storing the compressed lookup table, is 8 bytes for
storing the remainders, 2 bytes for the quotient, one byte
for storing a flag to determine which of the 8 nodes need to use a
quotient, and one additional byte for storing additional element
specific flags. Storing the lookup explicitly would require 8 ints, and therefore we obtain a 3x compression in storing the
lookup table.
Details on the octree-based meshing as well as its compression can be found in,
Low-constant Parallel Algorithms for Finite Element Simulations using Linear Octrees, Hari Sundar, Rahul S. Sampath, Santi S. Adavani, Christos Davatzikos, George Biros, Super Computing 2007.