Project a polynomial representation in elements in one dimension to its two halves in that direction.
For each parent element the projection on the two respective child elements (half intervals) are computed for one dimension.
Preliminary data layout and interface planning. It might be that we should rather split the index into the direction in which we perform the operation and all the other directions normal to that. For a dense matrix this may allow the compiler to detect the matrix multiply. However, here for the triangular matrix it is not so sure, whether this would be possible.
After discussions with Stephan Walter, it looks like the separate indices would most likely be better. Maybe, using explicit shaped arrays and therby allowing more dimensions in the input, while keeping the interface to two dimensions for all cases (the normal direction and all independent degrees of freedom). For vectorization on x86 it also is necessary to have a stride-1 access only in reading and writing. The rotation of data might not be the best option because of this. Instead, it may be that we need to have different routines for each direction. Or, maybe, we need to use the elements as first index and vectorize over those.
As we need to perform this operation in all dimensions, it would be good to shift the indices around. When doing this, we can stick to the same implementation for all directions, without the need to put any logic in here to decide on the current direction. In 3D we would end up with this chain: (x,y,z) -> split_element for Z -> (z,x,y) -> split_element for Y -> (y,z,x) -> split_element for X -> (x,y,z) Thus, the logic is that we perform the split on the last dimension, and cycle the indices in the output.
We can generalize this to arbitrary dimensions. In 2D it would look like this: (x,y) -> split_element for Y -> (y,x) -> split_element for X -> (x,y) And in 1D, we just need to perform one transformation: (x) -> split_element for X -> (x)
As we allow for a changed number of polynomial degrees in the input and output, we need to take care of different lengths for each direction. Thus, we need: the dimensionality, two 1D arrays with the length of this dimensionality to provide the number of degrees of freedom for each direction (once for the input, and once for the output).
We need: nDofs in the direction where the transformation is to be done and the nDofs for all normal directions.
Number of dimensions of the polynomial data.
Number degrees of freedom for each direction in parent_Data.
The first index of parent_data needs to have a length equal to the product of all inLen components. The splitting operation will be done in the last dimension.
Number degrees of freedom for each direction in child_Data.
The first index of child_data needs to have a length equal to the product of all outLen components. The data will be cyclicly exchanged. Thus, the last dimension in parent_data corresponds to the first in one in child_data and all other components are shifted once to the right.
Polynomial representation in the parent elements.
The first index are the degrees of freedom in elements, the second index are the elements. In the first index the shape of data has to be in the form (inLen(1), inLen(2), ... , inLen(nDims)). The splitting operation is performed on the last dimension in that data.
Computed projection of the polynomial representation in the child elements.
Again, the first index refers to the degrees of freedom, while the second index are the elements. There need to be twice as many elements as in the parent_data. Left childs are stored in iChild = (iParent2 - 1), and the right childs in iParent2.
In the first index the shape of the data has to be in the form (outLen(1), outLen(2), ... , outLen(nDims)), the data is rotated in comparison to parent_data and the splitted direction has to be the first one in child_data (while it was the last in parent_data), and all other dimensions are shifted by one to the right.
Whether to ignore high modes that exceed the target maximal polynomial degree.
This can be used as a simple lowpass filter that cuts off the highest modes in the parent elements prior to mapping to child elements.