In one-dimension - on a segment - a cubic interpolating function is defined by the values of the function and its derivative at both ends. This is almost the same in 2D inside a triangle, except that the values of the function and its 2 derivatives have to be defined at each triangle node.
The CubicTriInterpolator takes the value of the function at each node - provided by the user - and internally computes the value of the derivatives, resulting in a smooth interpolation. (As a special feature, the user can also impose the value of the derivatives at each node, but this is not supposed to be the common usage.)
This note is a bit technical and details how the cubic interpolation is computed.
The interpolation is based on a Clough-Tocher subdivision scheme of the triangulation mesh (to make it clearer, each triangle of the grid will be divided in 3 child-triangles, and on each child triangle the interpolated function is a cubic polynomial of the 2 coordinates). This technique originates from FEM (Finite Element Method) analysis; the element used is a reduced Hsieh-Clough-Tocher (HCT) element. Its shape functions are described in . The assembled function is guaranteed to be C1-smooth, i.e. it is continuous and its first derivatives are also continuous (this is easy to show inside the triangles but is also true when crossing the edges).
In the default case (kind ='min_E'), the interpolant minimizes a curvature energy on the functional space generated by the HCT element shape functions - with imposed values but arbitrary derivatives at each node. The minimized functional is the integral of the so-called total curvature (implementation based on an algorithm from - PCG sparse solver):
$$E(z) = \frac{1}{2} \int_{\Omega} \left( \left( \frac{\partial^2{z}}{\partial{x}^2} \right)^2 + \left( \frac{\partial^2{z}}{\partial{y}^2} \right)^2 + 2\left( \frac{\partial^2{z}}{\partial{y}\partial{x}} \right)^2 \right) dx\,dy$$
If the case kind ='geom' is chosen by the user, a simple geometric approximation is used (weighted average of the triangle normal vectors), which could improve speed on very large grids.
The triangulation to interpolate over.
Array of values, defined at grid points, to interpolate between.
Choice of the smoothing algorithm, in order to compute the interpolant derivatives (defaults to 'min_E'):
if 'min_E': (default) The derivatives at each node is computed to minimize a bending energy.
if 'geom': The derivatives at each node is computed as a weighted average of relevant triangle normals. To be used for speed optimization (large grids).
if 'user': The user provides the argument dz, no computation is hence needed.
If not specified, the Triangulation's default TriFinder will be used by calling .Triangulation.get_trifinder
.
Used only if kind ='user'. In this case dz must be provided as (dzdx, dzdy) where dzdx, dzdy are arrays of the same shape as z and are the interpolant first derivatives at the triangulation points.
Cubic interpolator on a triangular grid.
The following pages refer to to this document either explicitly or contain code examples using this.
matplotlib.tri.trirefine.UniformTriRefiner.refine_field
Hover to see nodes names; edges to Self not shown, Caped at 50 nodes.
Using a canvas is more power efficient and can get hundred of nodes ; but does not allow hyperlinks; , arrows or text (beyond on hover)
SVG is more flexible but power hungry; and does not scale well to 50 + nodes.
All aboves nodes referred to, (or are referred from) current nodes; Edges from Self to other have been omitted (or all nodes would be connected to the central node "self" which is not useful). Nodes are colored by the library they belong to, and scaled with the number of references pointing them