need a method to check the raylet, whether they are valid and aligned, passing through corresponding triangle, and blah blah.

Using different angular space discretization M, calculate the L^{\inf} error in terms of M, especially for point boundary source.

The Henyey Greenstein function is actually a conv-prod.

In 2d, FFT suffices to get it in O(nlog n) time, supposing there are n angles.
in 3d, it will consume, O(n^3 (log n)^2) time, supposing there are both n angles in the Euler angles.

which is superior than 2d's O(n^2) and 3d's O(n^4). But only approximately just one order better.

Now the test shows
hashset : 0.46s
set: 0.55s.

expansion of Greenstein.

needs to avoid the edges that are parallel to the ray directions.

This can be well-handled by taking a perturbation on the edge set.

I thought for a while, the storage for the directional solution is a little too much. But it is not yet known whether it is.

Using 128 directions maybe seem a lot, but it is not for many cases. Like localized source. Boundary sources need to be interpolated, thus, smooth source can be approximated better.

Another issue is the failure of interpolation. We interpolate from the physical space only on boundary, but it will cause some issue if the domain is not smooth, i.e. lack of input directions.

On a circle domain, this is not a problem. Of course, we need fine mesh on boundary.

rte solver in 2D only has 1 directional dependence, it is easy to handle the integration over velocity space.

For the meshes, unstructured mesh can be used directly from femm. Or just copy the relevant part of code and paste here.

For ray tracing, the geometry information has to be stored in prior. Storage efficiency needs to be considered.

It seems not hard when cascading ordering. But the ordering is not fast at all. It requires ordering for each direction.