File Mesh2DGrid.hpp
↰ Parent directory (src/interpolation)
Machinery to do interpolation on a semi-tensor mesh.
Definition (src/interpolation/Mesh2DGrid.hpp)
Detailed Description
Mark Krumholz
2024-06-19
Copyright (c) 2026 Mark Krumholz. All rights reserved.
The classes in this module and a few others solve a generic mathematical problem. We have a function whose value is specified at a set of data points that lie on a semi-tensor grid of the form
(x_{ij}, y_j), i = 0 … N-1, j = 0 … M-1
The points are non-decreasing, so x_{i,j} <= x_{i+1,j} and y_j <= y_{j+1}, but need not be strictly increasing. There also need not be any particular ordering of x_{i,j} and x_{i,j+1}, which implies that the outer edge of the mesh need not be convex, and that edges can be degenerate, i.e., (x_{ij}, y_j) = (x_{i+1,j}, y_j) = (x_{i,j+1}, y_{j+1}), so the ij - i+1,j and ij - i,j+1 edges are identical. Nor is the degeneracy guaranteed to be the same in all rows, i.e., we may have x_{i,j} == x_{i+1,j} but x_{i,j+1} != x_{i+1,j+1}.
Given a set of function values f_{ij} defined at the grid points, we wish to define an interpolating function f(x,y) that will give us the value of f at any point in the interior of the semi-regular grid. We wish to extract this in two forms. First, we want to be able to specify a value (x,y) that lies in the interior of the grid and get back a value of f. Second, we want to be able to specify a value x and get back a function f_x(y) that interpolates in the y direction at fixed x, or specify a value of y and get back a function f_y(x) that interpolates as fixed y.
In the context of stellar evolution, the grid consists of a set of evolutionary tracks. Each track is defined by an initial mass y_j = log m_j, and by a series of evolutionary points, each characterized by a time x_ij = log t_ij. The times t_ij are not the same for different initial masses, so in general x_{ij} != x_{i+1,j}. For each point in the evolutionary tracks, we have one or more stellar properties f_{ij}, and we want to be able to interpolate these properties, either by specifying a mass and an age directly, or by specifying an age and getting back an isochrone a function that returns stellar properties as a function of mass at fixed age.
To solve this problem, we note that the points (x_ij, y_j) define a series of quadrilaterals in the (x, y) plane. These quadrilaterals have parallel edges in the x direction, but non-parallel edges in the y direction. Graphically, the structure is something like this:
(x_{i,j+1}, y_{j+1}) (x_{i+1}, y_{i+1,j+1}) o——————-o / / / o—————————o (x_{ij}, y_j) (x_{i+1,j}, y_j})
Note that “degenerate” quadrilaterals are allowed, whereby either the top or bottom edge has length zero, and thus the shape is actually a triangle.
For convenience we will refer to the edges at constant y as “ribs”, and the edges that are not at constant y as “spines”. Let index (i,j) refer to the spine that lies above point (x_{ij}, y_j). The slope of this spine is
m_{ij} = (y_{j+1} - y_j) / (x_{i,j+1} - x_{ij}).
We can also define a distance coordinate along the spines. There are two possible ways to parmaeterize the coordinates in this direction. One is just by the y coordinate, but another choice is to define the distances to the grid points along the cell edges the define the spines recursively by
s_{i0} = 0 s_{i,j+1} = s_{i,j} + sqrt[ (x_{i,j+1} - x_{ij})^2 + (y_{j+1} - y_j)^2 ].
This can be generalized naturally to give a coordinate s for an arbitrary point (x, y) that lies on one of the spines:
s(x,y) = s_{i,jmax} + sqrt[ (x - x_{i,jmax})^2 + (y - y_{jmax})^2 ]
where jmax = max(j) s.t. y > y_j and the condition that (x, y) be on the spine is equivalent to requiring that
y - y_{jmax} = m_{i,jmax} (x - x_{i,jmax})
for some i.
Next, we define spline fits along both the horizontal and vertical edges. Specifically, consider the set of points along a rib
p_j = (x_{0j}, f_{0j}), (x_{1,j}, f_{1,j}), … (x_{N-1,j}, f_{N-1,j})
These points p_j can be used to define a spline function P_j(x) that accepts as an argument any x in [ x_{0j}, x_{N-1,j} ], and returns a function value at x. The order and form of the spline is arbitrary, but we will use a true spline, i.e., one that is constrained to go through each data point, rather than a basis spline.
We can similarly define the set of points in along the ribs. If we use the variable s defined above as our coordinate along the spines, this is:
q_i = (s_{i0}, f_{i0}), (s_{i1}, f_{i1}), … (s_{i,M-1}, f_{i,M-1})
and their corresponding spine spline functions Q_i(s(y)), where
s(y) = s(x(y), y) x(y) = x_{i,jmax} + (y - y_{jmax}) / m_{i,jmax}
Alternatively, if we use the y position as our coordinate across the spines the results are the same, except that we simply set s = y everywhere.
With these definitions in place, we can now define the procedures for interpolating to a single point or alone a line of fixed x. First consider generating an interpolating function along a line of fixed x. We find all locations where a line of fixed x intersects the mesh; these locations can be alone either spines or ribs. Let y_k be the set of intersection points, ordered so that y_{k+1} > y_k. For each intersection point, we define a corresponding interpolated function value f_k by evaluating either P(x) (if the point corresponds to intersecting a horizontal cell edge) or Q(s(y)) (if it corresponds to intersecting a vertical cell edge). Given the set (y_k, f_k), we can define a new interpolating function F(y) that allows us to evaluate the function at any y.
Again, depending on our choice of coordinate along the spines, we can also just set s = y in the above. The reason for using s = y as an option is that this option, together with an overall linear interpolation scheme, makes it possible to enforce strict monotonicity, in the sense that if some quantity is monotonic in x along all tracks, and we wish all interpolants to be similarly monotonic, then we can guarantee that this will be so by setting s = y and using linear interpolation. Any other choice does not strictly guarantee monotonic behavior. However, the price of this choice is an overall reduction in the level of accuracy in the interpolation.
The case of interpolating to a single point is analogous, except that we need not find all the intersection points y_k. We just find as many as needed on each side to have enough for the type of interpolating function we wish to use.
Finally, the case of creating a function that interpolates in x at fixed y is the same as well, except that the intersection points are all vertical.
The code to do all of this is broken up into two classes:
Mesh2DGrid this handles all the geometry of defining a grid, finding intersections of lines through it, etc. It does not know anything about interpolation, and just handles geometry. Mesh2dInterpolator this class wraps Mesh2DGrid, and handles all the interpolation machinery.
Includes
../utils/ThreadVec.hppcassertcmathcstddefcstdintcstdliblimitsmdspanspantupleutilityvector