I am attempting the following.
I have a method to produce one dimensional Legendre polynomials, I have a function from $\Phi: \mathbb{R}^3 \rightarrow \mathbb{R}$.
The idea is to construct a three dimensional basis of polynomials so that I can project $\Phi$ onto it.
Section 3.2 of this paper.
To that effect I coded an algorithm (see bottom). I have tested the functions called by this implementation, so I am highly confident they work. I am pretty certain that my current errors come from the provided code.
Specifically, I am noticing that no matter the input function, when I try pojecting my function on the basis, the result is always a sphere. My current toy input is a torus.
I am not 100% sure what I am doing wrong and was hoping someone had some insight.
pub fn project<F>(&mut self, degree: usize, mut function: F) -> Vec<S> where F: FnMut([S; 3]) -> S, { // We use the nodes for the 4 * degree scheme in legendre interpolation. // This is a heuristic described in "Hierarchichal HP-Adaptive Sigend Distance // fields." We may want to use a different quadrature scheme in the // future. let n = 4 * degree; let nodes = legendre_nodes_lut(n); let weights = legendre_weights_lut(n); // Iterator over all tuples (i,j,k) such that i + j + k <= degree. let bounded_tuple_iterator = BoundedTupleSum::<3>::new(degree); let total_basis_count = bounded_tuple_iterator.len(); // tensor product extension of the legendre basis, our nodes are now // X_I = x_i, x_j, x_k, where x_l is a one dimensional legendre root. // Our weights are w_I = w_i * w_j * w_k where w_l are the one dimensional // Gauss-Legendre quadrature weights. let mut values = vec![S::default(); n * n * n]; for multi_index in IterDimensions::new([n, n, n]) { let mut weight = S::from(1.).unwrap(); for i in multi_index { weight *= S::from(weights[i]).unwrap(); } let point = multi_index.map(|i| S::from(nodes[i]).unwrap()); let val = function(point); // Array based representation of the grid. We are just writing to the cell // correspoinding to the multi-index (i, j, k). values[multi_index[2] * n * n + multi_index[1] * n + multi_index[0]] = val * weight; } let mut coefficients = vec![S::default(); total_basis_count]; for multi_index in IterDimensions::new([n, n, n]) { let vw = values[multi_index[2] * n * n + multi_index[1] * n + multi_index[0]]; let quadrature_point = multi_index.map(|i| S::from(nodes[i]).unwrap()); let basis_eval = self.eval(quadrature_point, degree); for power_indices in bounded_tuple_iterator.clone() { let index = BoundedTupleSum::<3>::to_index(&power_indices); coefficients[index] += vw * basis_eval[index]; } } coefficients } ``` 
