"""Library-neutral spline airfoil runtime types."""
from __future__ import annotations
from copy import deepcopy
from typing import override
import numpy as np
from buffalo_core.numeric import as_float_array, as_float_scalar
from buffalo_core.typing import FloatArray, FloatInput, FloatScalar
from buffalo_wings.airfoil.internal.base import Airfoil
from buffalo_wings.airfoil.internal.bezier import BezierCurve2D
from buffalo_wings.airfoil.internal.runtime_common import ROOT_ABS_TOLERANCE
from buffalo_wings.airfoil.internal.runtime_types import (
AirfoilSurface,
SurfaceMappedValues,
)
from buffalo_wings.airfoil.internal.schema.spline import (
SplineAirfoilProvenanceSpec,
SplineAirfoilSpec,
SplineSurfaceSpec,
)
[docs]
class SplineAirfoil(Airfoil):
"""
Bezier-backed spline airfoil runtime.
Notes
-----
This name is intentionally backend-neutral.
The first implemented spline subtype stores upper and lower Bezier
curves directly while keeping the public runtime family name generic
enough for future B-spline and NURBS-backed implementations.
"""
def __init__(
self,
*,
upper: BezierCurve2D,
lower: BezierCurve2D,
provenance: SplineAirfoilProvenanceSpec | None = None,
) -> None:
"""
Initialize one spline-backed airfoil.
Parameters
----------
upper : BezierCurve2D
Upper-surface control curve ordered from leading edge to
trailing edge.
lower : BezierCurve2D
Lower-surface control curve ordered from leading edge to
trailing edge.
provenance : dict[str, object] | None, default=None
Optional provenance metadata for derived spline airfoils.
"""
super().__init__()
self._upper = upper
self._lower = lower
self._provenance = deepcopy(provenance)
@property
def upper(self) -> BezierCurve2D:
"""
Return the upper spline curve.
Returns
-------
BezierCurve2D
Upper-surface Bezier curve.
"""
return self._upper
@property
def lower(self) -> BezierCurve2D:
"""
Return the lower spline curve.
Returns
-------
BezierCurve2D
Lower-surface Bezier curve.
"""
return self._lower
@property
def provenance(self) -> SplineAirfoilProvenanceSpec | None:
"""
Return the optional spline provenance metadata.
Returns
-------
dict[str, object] | None
Deep-copied provenance metadata when available.
"""
return deepcopy(self._provenance)
@property
def spec(self) -> SplineAirfoilSpec:
"""
Return the persisted spline schema for this airfoil.
Returns
-------
SplineAirfoilSpec
Bezier-backed spline schema definition with upper and lower
control points ordered from leading edge to trailing edge.
Provenance metadata is deep-copied into the returned schema
when present.
"""
return SplineAirfoilSpec(
representation="bezier",
upper=SplineSurfaceSpec(
control_points=[
(as_float_scalar(point[0]), as_float_scalar(point[1]))
for point in self.upper.control_points
]
),
lower=SplineSurfaceSpec(
control_points=[
(as_float_scalar(point[0]), as_float_scalar(point[1]))
for point in self.lower.control_points
]
),
provenance=self.provenance,
)
[docs]
@override
def xy_from_u(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Evaluate the spline-backed airfoil coordinates.
Parameters
----------
u : FloatInput
Airfoil parameter values in ``[-1, 1]``.
Negative values evaluate the lower surface from trailing edge
toward the leading edge.
Nonnegative values evaluate the upper surface from leading edge
toward the trailing edge.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(x, y)`` of ``float64`` arrays matching the normalized
shape of ``u``.
"""
u_array = as_float_array(u)
x = as_float_array(np.empty_like(u_array))
y = as_float_array(np.empty_like(u_array))
upper_mask = u_array >= 0.0
lower_mask = ~upper_mask
if np.any(upper_mask):
x_upper, y_upper = self.upper.xy(u_array[upper_mask])
x[upper_mask] = x_upper
y[upper_mask] = y_upper
if np.any(lower_mask):
x_lower, y_lower = self.lower.xy(-u_array[lower_mask])
x[lower_mask] = x_lower
y[lower_mask] = y_lower
return x, y
[docs]
@override
def xy_u(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Evaluate first derivatives with respect to the airfoil parameter.
Parameters
----------
u : FloatInput
Airfoil parameter values in ``[-1, 1]``.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(dx/du, dy/du)`` of ``float64`` arrays matching the
normalized shape of ``u``.
Notes
-----
If ``u`` matches one of :meth:`breakpoints` exactly, this method
returns the ``minus``-side derivative limit.
"""
u_array = as_float_array(u)
x_u = as_float_array(np.empty_like(u_array))
y_u = as_float_array(np.empty_like(u_array))
analytic_mask = np.ones_like(u_array, dtype=bool)
flat_u = u_array.reshape(-1)
flat_x_u = x_u.reshape(-1)
flat_y_u = y_u.reshape(-1)
flat_mask = analytic_mask.reshape(-1)
breakpoints = self.breakpoints()
for index, value in enumerate(flat_u):
breakpoint_index = self._breakpoint_index(
as_float_scalar(value),
breakpoints=breakpoints,
)
if breakpoint_index is not None:
minus, _ = self.xy_u_breakpoint(index=breakpoint_index)
flat_x_u[index] = minus[0]
flat_y_u[index] = minus[1]
flat_mask[index] = False
if np.any(analytic_mask):
upper_mask = u_array >= 0.0
lower_mask = ~upper_mask
upper_mask &= analytic_mask
lower_mask &= analytic_mask
if np.any(upper_mask):
x_upper, y_upper = self.upper.xy_u(u_array[upper_mask])
x_u[upper_mask] = x_upper
y_u[upper_mask] = y_upper
if np.any(lower_mask):
x_lower, y_lower = self.lower.xy_u(-u_array[lower_mask])
x_u[lower_mask] = -x_lower
y_u[lower_mask] = -y_lower
return x_u, y_u
[docs]
@override
def xy_uu(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Evaluate second derivatives with respect to the airfoil parameter.
Parameters
----------
u : FloatInput
Airfoil parameter values in ``[-1, 1]``.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(d^2x/du^2, d^2y/du^2)`` of ``float64`` arrays matching the
normalized shape of ``u``.
Notes
-----
If ``u`` matches one of :meth:`breakpoints` exactly, this method
returns the ``minus``-side derivative limit.
"""
u_array = as_float_array(u)
x_uu = as_float_array(np.empty_like(u_array))
y_uu = as_float_array(np.empty_like(u_array))
analytic_mask = np.ones_like(u_array, dtype=bool)
flat_u = u_array.reshape(-1)
flat_x_uu = x_uu.reshape(-1)
flat_y_uu = y_uu.reshape(-1)
flat_mask = analytic_mask.reshape(-1)
breakpoints = self.breakpoints()
for index, value in enumerate(flat_u):
breakpoint_index = self._breakpoint_index(
as_float_scalar(value),
breakpoints=breakpoints,
)
if breakpoint_index is not None:
minus, _ = self.xy_uu_breakpoint(index=breakpoint_index)
flat_x_uu[index] = minus[0]
flat_y_uu[index] = minus[1]
flat_mask[index] = False
if np.any(analytic_mask):
upper_mask = u_array >= 0.0
lower_mask = ~upper_mask
upper_mask &= analytic_mask
lower_mask &= analytic_mask
if np.any(upper_mask):
x_upper, y_upper = self.upper.xy_uu(u_array[upper_mask])
x_uu[upper_mask] = x_upper
y_uu[upper_mask] = y_upper
if np.any(lower_mask):
x_lower, y_lower = self.lower.xy_uu(-u_array[lower_mask])
x_uu[lower_mask] = x_lower
y_uu[lower_mask] = y_lower
return x_uu, y_uu
[docs]
@override
def u_from_xi(
self,
xi: FloatInput,
*,
surface: AirfoilSurface,
) -> FloatArray:
"""
Convert surface-local ``xi`` coordinates to native parameters.
Parameters
----------
xi : FloatInput
Surface-local coordinates in ``[0, 1]`` measured from the
leading edge to the trailing edge.
surface : {"lower", "upper"}
Surface to evaluate.
Returns
-------
FloatArray
Signed native spline-airfoil parameters matching ``xi`` on the
selected surface.
Notes
-----
The spline runtime uses the linear mapping ``u = +/- xi``, with the
sign determined by ``surface``.
"""
return self._u_from_xi_signed_linear(xi, surface=surface)
[docs]
@override
def xi_from_u(self, u: FloatInput) -> SurfaceMappedValues:
"""
Convert native parameters to surface-local ``xi`` coordinates.
Parameters
----------
u : FloatInput
Signed native spline-airfoil parameters in ``[-1, 1]``.
Returns
-------
SurfaceMappedValues
Surface-local ``xi`` values and upper-surface membership flags.
Notes
-----
The spline runtime uses the linear mapping ``xi = |u|``.
"""
return self._xi_from_u_signed_linear(u)
[docs]
@override
def breakpoints(self) -> list[FloatScalar]:
"""
Return the airfoil boundary and leading-edge breakpoints.
Returns
-------
list[FloatScalar]
Ordered breakpoint list ``[-1.0, 0.0, 1.0]``.
The endpoints are boundary breakpoints and ``0.0`` is the
leading-edge join between lower and upper surfaces.
"""
return [
as_float_scalar(-1.0),
as_float_scalar(0.0),
as_float_scalar(1.0),
]
@staticmethod
def _breakpoint_index(
value: FloatScalar,
*,
breakpoints: list[FloatScalar],
) -> int | None:
"""
Return the matching spline breakpoint index within tolerance.
Parameters
----------
value : FloatScalar
Native spline airfoil parameter to classify.
breakpoints : list[FloatScalar]
Ascending native breakpoint parameters.
Returns
-------
int | None
Matching breakpoint index when ``value`` lies within the root
tolerance of one breakpoint, otherwise ``None``.
"""
matches = np.flatnonzero(
np.isclose(
as_float_array(breakpoints),
value,
atol=ROOT_ABS_TOLERANCE,
rtol=0.0,
)
)
if matches.size == 0:
return None
return int(matches[0])
[docs]
def xy_u_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided first derivatives at one breakpoint.
Parameters
----------
index : int
Breakpoint index into :meth:`breakpoints`.
Returns
-------
tuple[tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]]
``((x_u_minus, y_u_minus), (x_u_plus, y_u_plus))`` for the
selected breakpoint.
Endpoint breakpoints return the same boundary value for both
sides because only one in-domain side exists.
"""
u_breakpoint = self.breakpoints()[index]
if u_breakpoint <= -1.0:
boundary = self._curve_values(self.lower.xy_u(1.0), sign=-1.0)
return boundary, boundary
if u_breakpoint >= 1.0:
boundary = self._curve_values(self.upper.xy_u(1.0))
return boundary, boundary
minus = self._curve_values(self.lower.xy_u(0.0), sign=-1.0)
plus = self._curve_values(self.upper.xy_u(0.0))
return minus, plus
[docs]
def xy_uu_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided second derivatives at one breakpoint.
Parameters
----------
index : int
Breakpoint index into :meth:`breakpoints`.
Returns
-------
tuple[tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]]
``((x_uu_minus, y_uu_minus), (x_uu_plus, y_uu_plus))`` for the
selected breakpoint.
Endpoint breakpoints return the same boundary value for both
sides because only one in-domain side exists.
"""
u_breakpoint = self.breakpoints()[index]
if u_breakpoint <= -1.0:
boundary = self._curve_values(self.lower.xy_uu(1.0))
return boundary, boundary
if u_breakpoint >= 1.0:
boundary = self._curve_values(self.upper.xy_uu(1.0))
return boundary, boundary
minus = self._curve_values(self.lower.xy_uu(0.0))
plus = self._curve_values(self.upper.xy_uu(0.0))
return minus, plus
[docs]
@override
def xy_s_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided arc-length derivatives at one breakpoint.
Parameters
----------
index : int
Breakpoint index into :meth:`breakpoints`.
Returns
-------
tuple[tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]]
``((x_s_minus, y_s_minus), (x_s_plus, y_s_plus))`` for the
selected breakpoint.
Notes
-----
This method composes the exact arc-length tangent values from the
exact native breakpoint derivatives returned by
:meth:`xy_u_breakpoint`.
If one side has zero native speed, the generic sampled fallback is
retained for that breakpoint.
"""
minus_u, plus_u = self.xy_u_breakpoint(index=index)
if np.isclose(np.hypot(*minus_u), 0.0) or np.isclose(
np.hypot(*plus_u),
0.0,
):
return super().xy_s_breakpoint(index=index)
return (
self._xy_s_from_native_breakpoint(minus_u),
self._xy_s_from_native_breakpoint(plus_u),
)
[docs]
@override
def xy_ss_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided arc-length second derivatives at one breakpoint.
Parameters
----------
index : int
Breakpoint index into :meth:`breakpoints`.
Returns
-------
tuple[tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]]
``((x_ss_minus, y_ss_minus), (x_ss_plus, y_ss_plus))`` for the
selected breakpoint.
Notes
-----
This method composes the exact arc-length curvature-vector values from
the exact native breakpoint derivatives returned by
:meth:`xy_u_breakpoint` and :meth:`xy_uu_breakpoint`.
If one side has zero native speed, the generic sampled fallback is
retained for that breakpoint.
"""
minus_u, plus_u = self.xy_u_breakpoint(index=index)
if np.isclose(np.hypot(*minus_u), 0.0) or np.isclose(
np.hypot(*plus_u),
0.0,
):
return super().xy_ss_breakpoint(index=index)
minus_uu, plus_uu = self.xy_uu_breakpoint(index=index)
return (
self._xy_ss_from_native_breakpoint(minus_u, minus_uu),
self._xy_ss_from_native_breakpoint(plus_u, plus_uu),
)
@staticmethod
def _xy_s_from_native_breakpoint(
native_u: tuple[FloatScalar, FloatScalar],
) -> tuple[FloatScalar, FloatScalar]:
"""
Return one exact arc-length tangent from native derivatives.
Parameters
----------
native_u : tuple[FloatScalar, FloatScalar]
Exact one-sided native derivative tuple ``(x_u, y_u)``.
Returns
-------
tuple[FloatScalar, FloatScalar]
Exact arc-length tangent ``(x_s, y_s)`` on the same breakpoint
side.
"""
x_u, y_u = native_u
speed = as_float_scalar(np.hypot(x_u, y_u))
return (
as_float_scalar(x_u / speed),
as_float_scalar(y_u / speed),
)
@staticmethod
def _xy_ss_from_native_breakpoint(
native_u: tuple[FloatScalar, FloatScalar],
native_uu: tuple[FloatScalar, FloatScalar],
) -> tuple[FloatScalar, FloatScalar]:
"""
Return one exact arc-length second derivative from native data.
Parameters
----------
native_u : tuple[FloatScalar, FloatScalar]
Exact one-sided native derivative tuple ``(x_u, y_u)``.
native_uu : tuple[FloatScalar, FloatScalar]
Exact one-sided native second-derivative tuple
``(x_uu, y_uu)``.
Returns
-------
tuple[FloatScalar, FloatScalar]
Exact arc-length second derivative ``(x_ss, y_ss)`` on the same
breakpoint side.
"""
x_u, y_u = native_u
x_uu, y_uu = native_uu
speed_sq = as_float_scalar(x_u**2 + y_u**2)
speed_pow4 = as_float_scalar(speed_sq**2)
projection = as_float_scalar(x_u * x_uu + y_u * y_uu)
return (
as_float_scalar(x_uu / speed_sq - x_u * projection / speed_pow4),
as_float_scalar(y_uu / speed_sq - y_u * projection / speed_pow4),
)
@staticmethod
def _curve_values(
values: tuple[FloatArray, FloatArray],
*,
sign: FloatScalar = 1.0,
) -> tuple[FloatScalar, FloatScalar]:
"""
Convert one Bezier evaluation pair into scalar breakpoint values.
Parameters
----------
values : tuple[FloatArray, FloatArray]
Pair ``(x_values, y_values)`` returned by one Bezier-curve
evaluation call.
This helper is used only for scalar breakpoint evaluations.
sign : FloatScalar, default=1.0
Sign convention applied to both returned components.
Lower-surface first derivatives use ``-1.0`` because the public
airfoil parameter runs from trailing edge to leading edge on the
lower branch.
Returns
-------
tuple[FloatScalar, FloatScalar]
Scalar pair ``(x_component, y_component)``.
"""
return (
sign * as_float_scalar(values[0]),
sign * as_float_scalar(values[1]),
)