LeanderSchlegel
commented
2 years ago Thank you for your hint, that is a good point! I wrote a small testscript and if I did it correctly, found out that the current implementation indeed works for both, Vector3f and Vector3d as well as Grid1f and Grid1d. Like you already assumed, Vector3d are getting downcasted to Vector3f precision. I am not sure if floats are cast to doubles in the interpolation and we just do not see it with the test, however, I think this would be no problem. I think your idea to improve the implementation of the genericization by using the template is definitely a good idea, however, for now I would think the current implementation is sufficient, because if I see it correctly, all functions in GridTools.h only work for float precision grids (although one could set values on double precision grids using setValue).
`#Datatype Tests from crpropa import *
origin = Vector3d(0.) n = 100 spacing = 1 * pc
init scalar grids
grid_1f = Grid1f(origin, n, spacing) grid_1d = Grid1d(origin, n, spacing)
init vector grids
grid_3f = Grid3f(origin, n, spacing) grid_3d = Grid3d(origin, n, spacing)
fill values
grid_1f.setValue(0,0,0,0.12345678910111213141516171819) grid_1d.setValue(0,0,0,0.12345678910111213141516171819)
grid_3f.setValue(0,0,0,Vector3f(0.12345678910111213141516171819)) grid_3d.setValue(0,0,0,Vector3d(0.12345678910111213141516171819))
retrieve values using tricubic interpolation
grid_1f.setInterpolationType(TRICUBIC) grid_1d.setInterpolationType(TRICUBIC) grid_3f.setInterpolationType(TRICUBIC) grid_3d.setInterpolationType(TRICUBIC)
print grid_1f.getValue(0,0,0), grid_1f.interpolate(Vector3d(spacing/2.,spacing/2.,spacing/2.)) print grid_1d.getValue(0,0,0), grid_1d.interpolate(Vector3d(spacing/2.,spacing/2.,spacing/2.)) print grid_3f.getValue(0,0,0)[0], grid_3f.interpolate(Vector3d(spacing/2.,spacing/2.,spacing/2.))[0] print grid_3d.getValue(0,0,0)[0], grid_3d.interpolate(Vector3d(spacing/2.,spacing/2.,spacing/2.))[0]
note: taking the [0] component, because printing a Vector3d instead of one entry isolated, seems to round its entries.
`
As far as I can tell, there are currently (as of #1) two implementations of tricubicInterpolate, one for double and one for Vector3f. Do we want to make a generic implementation that covers Vector3d, float, and everything else people might throw at it?
It's possible that the current method for genericization might cause floats to be cast to doubles, and perhaps Vector3d's to Vector3f's, so that the implementation currently works for these two types, too. (But maybe not in ways that we want.) This might be something to check. We should also consider whether there is a better way to implement this genericization (using the template syntax more directly), which might also make it easier to reason about which implementation is selected. However, it might be best to leave this for the normal PR.
Open questions: