Lagrange point interpolation

Point interpolation evaluates field values and field first derivative value at any point of the space (2D or 3D). The field value must be known at the nodes of a Cartesian grid. No regularity of the grid spacing is implied.

This module is intended to provide flexible interpolations routines that works in any configuration and for any order. For more optimized routines, see other modules such as node_interpolation.

To evaluate \( f \) at some point \( (x,y,z) \), the module evaluates Lagrange polynomials \( L_i \) of order \( p \) using a sub-mesh of size \( p \times p \) (in 2D) size \( p \times p \times p \) (in 3D) placed over \( (x,y,z) \). The value \( f(x,y,z) \) is computed as follows:

\[ \tilde{f}(x,y,z) = \sum_{ijk} f(x_i,y_j,z_k) L_i(x) L_j(y) L_k(z), \]

where the sum is performed over the sub-mesh. The first partial derivative along the first axis is computed as follows:

\[ \frac{\partial\tilde{f}}{\partial x}(x,y,z) = \sum_{ijk} f(x_i,y_j,z_k) \frac{\partial L_i}{\partial x}(x) L_j(y) L_k(z), \]

and likewise for the first partial derivatives along the other axes.

The placement of the sub-mesh can be centered around \( (x,y,z) \) or forced to some place. Interpolation order \( p \) cannot exceed \( \min\{nx,ny\} \) (in 2D) \( \min\{nx,ny,nz\} \) (in 3D).

Application Programming Interface

Point interpolation is divided in two parts: in the first prepare part, interpolation coefficients are computed and stored, and in the second exec part, field (derivative) values are evaluated. It is interesting to notice that the prepare part depends only on the mesh positions, on the point position, and on interpolation options and thus can be done once for several evaluations.

A. Preparation

Import the following modules is the current execution unit:

use mod_point_interpolation_lagrange
use mod_point_interpolation_prepare

Declare the following data variable to store interpolation coefficients:

type(interpolation_data) :: data

A.1. Memory allocation

To interpolate to the order p, on a mesh of size nx, ny, nz:

call interpolation_init(data, p, nx, ny, nz)

A.2. Sub-mesh evaluation

There is two ways to define the sub-mesh.

Centered — Just call

call interpolation_set_stencil_centered(data, xxi, yyi, zzi, xo, yo, zo)

where xxi, yyi, zzi are rank-1 array containing the mesh positions and xo, yo, zo are the coordinates of the interpolation point. The sub-mesh may be not quite centered when xo, yo, zo is close to the mesh boundaries, especially when p is high.

Oriented — To force the sub-mesh to start from the node of indexes i_pi, j_pi, k_pi and to go towards the direction tx, ty, tz, use:

call interpolation_set_stencil_oriented(data, i_pi, j_pi, k_pi, tx, ty, tz)

A.3. Computation of interpolation or derivation coefficients

Finish the preparation part by calling the appropriate method:

Interpolation

call interpolation_set_coefficients(data, xxi, yyi, zzi, xo, yo, zo)

Derivation along x

call interpolation_set_coefficients_dx(data, xxi, yyi, zzi, xo, yo, zo)

Derivation along y

call interpolation_set_coefficients_dy(data, xxi, yyi, zzi, xo, yo, zo)

Derivation along z

call interpolation_set_coefficients_dz(data, xxi, yyi, zzi, xo, yo, zo)

From this point, data is ready to be used in the exec part.

B. Interpolations, derivations

Import the following modules is the current execution unit:

use mod_point_interpolation_exec

To evaluate the field at xo, yo, zo and store it in uo, do

call interpolation_exec(data, ui, uo)

where ui is a rank-3 array that contains field values at the mesh nodes.

C. Finalization

To destroy the interpolation_data instance:

call interpolation_finalize(data)