Add mesh.interpolate_line(mesh, xi, axis)

CHANGELOG.md

@@ -8,6 +8,7 @@ All notable changes to this project will be documented in this file. The format
 - Add `Field(..., dtype=None)` to cast the array with field values to a specified type.
 - Add `Field.extract(dtype=None)` to cast the extracted field gradient/interpolated values to a specified type.
 - Add `hello_world()` to print the lines of a minimal-working-example to the console.
+- Add `mesh.interpolate_line(mesh, xi, axis)` to interpolate a line mesh.
 
 ### Changed
 - Change the internal initialization of `field = Field(region, values=1, dtype=None)` values from `field.values = np.ones(shape) * values` to `field = np.full(shape, fill_value=values, dtype=dtype)`. This enforces `field = Field(region, values=1)` to return the gradient array with data-type `int` which was of type `float` before.

docs/felupe/mesh.rst

@@ -55,6 +55,7 @@ This module contains meshing-related classes and functions. Standalone mesh-tool
    mesh.dual
    mesh.stack
    mesh.read
+   mesh.interpolate_line
 
 **Detailed API Reference**
 
@@ -114,4 +115,4 @@ This module contains meshing-related classes and functions. Standalone mesh-tool
    :show-inheritance:
 
 .. automodule:: felupe.mesh
-   :members: expand, translate, rotate, revolve, sweep, mirror, concatenate, runouts, triangulate, convert, collect_edges, collect_faces, collect_volumes, add_midpoints_edges, add_midpoints_faces, add_midpoints_volumes, flip, fill_between, dual, stack, merge_duplicate_points, merge_duplicate_cells, read
+   :members: expand, translate, rotate, revolve, sweep, mirror, concatenate, runouts, triangulate, convert, collect_edges, collect_faces, collect_volumes, add_midpoints_edges, add_midpoints_faces, add_midpoints_volumes, flip, fill_between, dual, stack, merge_duplicate_points, merge_duplicate_cells, read, interpolate_line

src/felupe/mesh/__init__.py

@@ -20,6 +20,7 @@
     RectangleArbitraryOrderQuad,
     Triangle,
 )
+from ._interpolate import interpolate_line
 from ._line_rectangle_cube import cube_hexa as _cube_hexa
 from ._line_rectangle_cube import line_line as _line_line
 from ._line_rectangle_cube import rectangle_quad as _rectangle_quad
@@ -58,6 +59,7 @@
     "expand",
     "flip",
     "fill_between",
+    "interpolate_line",
     "mirror",
     "merge_duplicate_points",
     "merge_duplicate_cells",

src/felupe/mesh/_interpolate.py

@@ -0,0 +1,98 @@
+# -*- coding: utf-8 -*-
+"""
+This file is part of FElupe.
+
+FElupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+FElupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+
+
+def interpolate_line(mesh, xi, axis, Interpolator=None, **kwargs):
+    r"""Return an interpolated line mesh from an existing line mesh with provided
+    interpolation points on a given axis.
+
+    Parameters
+    ----------
+    mesh : Mesh
+        A :class:`~felupe.Mesh` with cell-type ``line``.
+    xi : ndarray
+        Points at which to interpolate data.
+    axis : int
+        The axis of the points at which the data will be interpolated.
+    Interpolator : callable or None, optional
+        An interpolator class (default is :class:`scipy.interpolate.PchipInterpolator`).
+    **kwargs : dict, optional
+        Optional keyword arguments are passed to the interpolator.
+
+    Returns
+    -------
+    Mesh
+        A new line mesh with interpolated points.
+
+    Examples
+    --------
+    ..  pyvista-plot::
+        :context:
+        :force_static:
+
+        >>> import felupe as fem
+        >>> import numpy as np
+        >>>
+        >>> mesh = fem.mesh.Line(n=5).expand(n=1)
+        >>> t = mesh.x.copy()
+        >>> mesh.points[:, 0] = np.sin(np.pi / 2 * t)
+        >>> mesh.points[:, 1] = np.cos(np.pi / 2 * t)
+        >>>
+        >>> mesh.plot(style="points", color="black").show()
+
+    ..  pyvista-plot::
+        :context:
+        :force_static:
+
+        >>> mesh_new = interpolate_line(mesh, xi=np.linspace(0, 1, 101), axis=1)
+        >>>
+        >>> mesh_new.plot(style="points", color="black").show()
+    """
+
+    # line connectivity
+    # concatenate the first point of each cell and the last point of the last cell
+    line = np.concatenate([mesh.cells[:, :-1].ravel(), mesh.cells[-1, -1:]])
+
+    # independent spline variable
+    points = mesh.points[line, axis]
+    ascending = np.argsort(points)
+
+    # dependent spline variable(s)
+    mask = np.ones(mesh.dim, dtype=bool)
+    mask[axis] = False
+
+    axes = np.arange(mesh.dim)
+    values = mesh.points[line.reshape(-1, 1), axes[mask]]
+
+    # spline interpolator
+    if Interpolator is None:
+        from scipy.interpolate import PchipInterpolator as Interpolator
+
+    # create a spline
+    spline = Interpolator(points[ascending], values[ascending], **kwargs)
+
+    # evaluation points for the independent spline variable
+    points_new = np.zeros((len(xi), mesh.dim))
+    points_new[:, axis] = xi
+    points_new[:, axes[mask]] = spline(xi)
+
+    cells_new = np.repeat(np.arange(len(xi)), 2)[1:-1].reshape(-1, 2)
+
+    return type(mesh)(points_new, cells_new, cell_type="line")

tests/test_mesh.py

@@ -553,6 +553,20 @@ def test_expand():
     assert np.allclose(cube.points.max(axis=0), [1, 2, 3])
 
 
+def test_interpolate_line():
+    import felupe as fem
+
+    mesh = fem.mesh.Line(n=5).expand(n=1).expand(n=1)
+    t = mesh.x.copy()
+    mesh.points[:, 0] = np.sin(np.pi / 2 * t)
+    mesh.points[:, 1] = np.cos(np.pi / 2 * t)
+
+    xi = np.linspace(0, 1, 101)
+    mesh_new = fem.mesh.interpolate_line(mesh, xi=xi, axis=1)
+
+    assert mesh_new.npoints == len(xi)
+
+
 if __name__ == "__main__":
     test_meshes()
     test_mirror()
@@ -572,3 +586,4 @@ def test_expand():
     test_mesh_update()
     test_modify_corners()
     test_expand()
+    test_interpolate_line()