"""Base airfoil class and core runtime state."""
from __future__ import annotations
from abc import abstractmethod
from typing import Literal, 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 scipy.optimize import root_scalar
from .curve import Curve
from .parameter_validation import (
surface_is_upper,
validate_curve_u,
validate_xi,
)
from .runtime_common import (
DOMAIN_ABS_TOLERANCE,
LOWER_SURFACE_BRACKET,
ROOT_ABS_TOLERANCE,
ROOT_MAX_ITERATION,
U_LEADING_EDGE,
U_MAX,
U_MIN,
UPPER_SURFACE_BRACKET,
X_MIN,
)
from .runtime_types import (
AirfoilCamberResult,
AirfoilSurface,
SurfaceMappedValues,
)
from .schema import AirfoilDefinitionSpec
CamberCurveSpacing = Literal["uniform", "cosine"]
_MIN_CAMBER_SAMPLES = 2
_DEFAULT_CAMBER_SAMPLES = 81
[docs]
class Airfoil(Curve):
"""
Base class for airfoil specific geometries.
Airfoils are represented in a normalized local section frame with the
leading edge at ``(0, 0)`` and the nominal trailing-edge midpoint at
``(1, 0)``.
Concrete airfoil families are expected to follow this unit-chord
convention.
Note that some airfoils might have points that are
before the leading edge point because of their definition.
The airfoil coordinates and their derivatives can be queried using a
parameterization, ``u``, that uses ``[-1, 0]`` for the lower surface,
running from trailing edge to leading edge, and ``[0, 1]`` for the upper
surface, running from leading edge to trailing edge.
This is called the curve parameter, and it provides a smooth
parameterization across the full airfoil.
Coordinates can also be queried using the :py:class:`Curve` arc-length
parameterization interface, ``s``.
The final way that coordinates can be queried is through a chord-like
parameterization for the upper and lower surfaces, ``xi``.
For this both surfaces leading to trailing edges are mapped to ``[0, 1]``.
The arc-length parameterization uses surface distance measured from the
lower trailing edge to the upper trailing edge.
Arc-length queries are more expensive because the mapping from surface
distance to the curve parameter is not available in closed form
for general airfoil shapes.
Airfoils inherit the :class:`Curve` breakpoint contract.
The ordinary derivative evaluators return the ``minus``-side value when a
query lands exactly on a reported breakpoint, while the paired
``*_breakpoint`` methods expose both one-sided values explicitly.
"""
def __init__(self) -> None:
"""
Initialize the cached airfoil arc-length state.
Notes
-----
Subclasses should call this initializer so cached geometric values
remain synchronized with the current airfoil shape.
"""
super().__init__()
self._lower_surface_length_cache: FloatScalar | None = None
self._upper_surface_length_cache: FloatScalar | None = None
#
# Spec handling
#
@property
def spec(self) -> AirfoilDefinitionSpec:
"""
Return the schema definition used to create this airfoil.
Returns
-------
AirfoilDefinitionSpec
Serialized airfoil definition that can recreate this runtime
object.
Raises
------
NotImplementedError
If the concrete airfoil type does not preserve its source spec.
Notes
-----
For schema-backed runtime families that participate in the current
round-trip contract, this property preserves the original supported
schema form exactly rather than normalizing it to a merely equivalent
definition.
Placeholder or not-yet-constructable families may still raise
``NotImplementedError`` until their schema contract is defined.
"""
raise NotImplementedError(
f"{type(self).__name__} does not provide a source spec."
)
[docs]
def to_spec(self) -> AirfoilDefinitionSpec:
"""
Return the schema definition needed to recreate this airfoil.
Returns
-------
AirfoilDefinitionSpec
Serialized airfoil definition that can recreate this runtime
object.
Notes
-----
For runtime families covered by the current schema round-trip
contract, this returns the same schema content as :attr:`spec`.
"""
return self.spec
#
# Cached geometric properties
#
@property
@override
def length(self) -> FloatScalar:
"""
Return the full airfoil surface length.
Returns
-------
FloatScalar
Total airfoil surface length measured from the lower trailing
edge to the upper trailing edge.
"""
return as_float_scalar(
self._lower_surface_length() + self._upper_surface_length()
)
def _lower_surface_length(self) -> FloatScalar:
"""
Return the cached lower-surface arc length.
Returns
-------
FloatScalar
Arc length measured from the lower trailing-edge breakpoint at
``u=-1`` to the leading edge at ``u=0``.
"""
if self._lower_surface_length_cache is None:
self._lower_surface_length_cache = as_float_scalar(
self.arc_length(-1.0, 0.0)
)
return as_float_scalar(self._lower_surface_length_cache)
def _upper_surface_length(self) -> FloatScalar:
"""
Return the cached upper-surface arc length.
Returns
-------
FloatScalar
Arc length measured from the leading edge at ``u=0`` to the
upper trailing-edge breakpoint at ``u=1``.
"""
if self._upper_surface_length_cache is None:
self._upper_surface_length_cache = as_float_scalar(
self.arc_length(0.0, 1.0)
)
return as_float_scalar(self._upper_surface_length_cache)
#
# Airfoil reference geometry
#
[docs]
def chord(self) -> FloatScalar:
"""
Return the airfoil chord length.
Returns
-------
FloatScalar
Distance between the leading-edge reference and trailing-edge
midpoint reference.
"""
xle, yle = self.leading_edge()
xte, yte = self.trailing_edge()
return as_float_scalar(np.hypot(xte - xle, yte - yle))
[docs]
def leading_edge(self) -> tuple[FloatScalar, FloatScalar]:
"""
Return the leading-edge location.
Returns
-------
tuple[FloatScalar, FloatScalar]
``(x, y)`` location of the leading-edge reference point.
"""
x_le, y_le = self.xy_from_u(0.0)
return as_float_scalar(x_le), as_float_scalar(y_le)
[docs]
def trailing_edge(self) -> tuple[FloatScalar, FloatScalar]:
"""
Return the midpoint of the trailing-edge points.
Returns
-------
tuple[FloatScalar, FloatScalar]
``(x, y)`` location of the trailing-edge midpoint reference.
"""
xl, yl = self.xy_from_u(-1.0)
xu, yu = self.xy_from_u(1.0)
return (
as_float_scalar(0.5 * (xl + xu)),
as_float_scalar(0.5 * (yl + yu)),
)
[docs]
def camber_curve(
self,
*,
num_points: int = _DEFAULT_CAMBER_SAMPLES,
spacing: CamberCurveSpacing = "cosine",
) -> AirfoilCamberResult:
"""
Return a camber-curve representation for this airfoil.
Parameters
----------
num_points : int, default=81
Number of shared surface samples to use when an approximate
camber line must be derived from the airfoil geometry.
spacing : {"uniform", "cosine"}, default="cosine"
Spacing rule used for the shared surface-local sample locations
in the approximate extraction path.
Returns
-------
AirfoilCamberResult
Exact or approximate camber-curve result for this airfoil.
Raises
------
ValueError
If ``num_points`` or ``spacing`` is invalid for the approximate
extraction path.
"""
if num_points < _MIN_CAMBER_SAMPLES:
msg = "num_points must be at least 2."
raise ValueError(msg)
xi = self._camber_curve_xi_samples(
num_points=num_points,
spacing=spacing,
)
upper_x, upper_y = self.xy_from_xi(xi, surface="upper")
lower_x, lower_y = self.xy_from_xi(xi, surface="lower")
x_mid = as_float_array(0.5 * (upper_x + lower_x))
y_mid = as_float_array(0.5 * (upper_y + lower_y))
x_mid[0] = 0.0
x_mid[-1] = 1.0
from .approximate_camber import ApproximateCamberCurve
return AirfoilCamberResult(
curve=ApproximateCamberCurve(
x_samples=x_mid,
y_samples=y_mid,
),
mode="approximate",
source_family=type(self).__name__,
)
#
# Curve airfoil differential geometry
#
[docs]
def dydx(self, u: FloatInput) -> FloatArray:
"""
Return the surface slope at curve parameter locations.
Parameters
----------
u : FloatInput
Airfoil parameters.
Returns
-------
FloatArray
Surface slope values ``dy/dx`` evaluated at ``u``.
"""
x_u, y_u = self.xy_u(u)
with np.errstate(divide="ignore", invalid="ignore"):
return y_u / x_u
[docs]
def d2ydx2(self, u: FloatInput) -> FloatArray:
"""
Return the second surface derivative at curve parameter locations.
Parameters
----------
u : FloatInput
Airfoil parameters.
Returns
-------
FloatArray
Second derivative values ``d^2y/dx^2`` evaluated at ``u``.
"""
x_u, y_u = self.xy_u(u)
x_uu, y_uu = self.xy_uu(u)
with np.errstate(divide="ignore", invalid="ignore"):
return (x_u * y_uu - y_u * x_uu) / x_u**3
#
# Surface parameterization API
#
[docs]
@abstractmethod
def u_from_xi(
self, xi: FloatInput, *, surface: AirfoilSurface
) -> FloatArray:
"""
Convert one-surface ``xi`` coordinates to curve 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
Curve parameters matching ``xi`` on the selected surface.
Notes
-----
Concrete airfoil families define this mapping because the surface
coordinate ``xi`` is airfoil-specific and need not be a sign-only
transformation of the curve parameter ``u``.
"""
[docs]
@abstractmethod
def xi_from_u(self, u: FloatInput) -> SurfaceMappedValues:
"""
Convert curve airfoil parameters to surface-local ``xi`` values.
Parameters
----------
u : FloatInput
Curve airfoil parameters in ``[-1, 1]``.
Returns
-------
SurfaceMappedValues
Surface-local ``xi`` values and upper-surface membership flags.
Notes
-----
Concrete airfoil families define this mapping because ``xi`` need not
equal ``|u|`` for every airfoil parameterization.
"""
[docs]
def xy_from_xi(
self, xi: FloatInput, *, surface: AirfoilSurface
) -> tuple[FloatArray, FloatArray]:
"""
Return one-surface coordinates at surface-local ``xi`` locations.
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
-------
tuple[FloatArray, FloatArray]
Tuple ``(x, y)`` of ``float64`` arrays matching the normalized
shape of ``xi``.
"""
return self.xy_from_u(self.u_from_xi(xi, surface=surface))
[docs]
def slope_from_xi(
self, xi: FloatInput, *, surface: AirfoilSurface
) -> FloatArray:
"""
Return one-surface slope values at surface-local ``xi`` locations.
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
Surface slope values ``dy/dx`` on the selected surface.
"""
return self.dydx(self.u_from_xi(xi, surface=surface))
[docs]
def curvature_from_xi(
self, xi: FloatInput, *, surface: AirfoilSurface
) -> FloatArray:
"""
Return one-surface curvature values at surface-local ``xi`` locations.
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
Surface-oriented curvature values on the selected surface.
"""
curvature = self.k(self.u_from_xi(xi, surface=surface))
if not surface_is_upper(surface):
curvature = -curvature
return curvature
#
# Surface parameterization shared helpers
#
@staticmethod
def _u_from_xi_signed_linear(
xi: FloatInput,
*,
surface: AirfoilSurface,
) -> FloatArray:
"""
Convert ``xi`` to curve parameters using ``u = +/- xi``.
Notes
-----
This helper exists for airfoil families whose curve parameter uses
the same leading-edge to trailing-edge coordinate magnitude on both
surfaces and differs only by sign between lower and upper branches.
"""
xi_array = validate_xi(xi)
return xi_array if surface_is_upper(surface) else -xi_array
@staticmethod
def _xi_from_u_signed_linear(u: FloatInput) -> SurfaceMappedValues:
"""
Convert curve parameters to ``xi`` using ``xi = |u|``.
Notes
-----
This helper exists for airfoil families whose ``xi`` coordinate is
the absolute value of the signed curve parameter.
"""
u_array = validate_curve_u(u)
return SurfaceMappedValues(
value=np.abs(u_array),
is_upper=u_array >= U_LEADING_EDGE,
)
@staticmethod
def _camber_curve_xi_samples(
*,
num_points: int,
spacing: CamberCurveSpacing,
) -> FloatArray:
"""
Build shared surface-local samples for approximate camber extraction.
Parameters
----------
num_points : int
Number of sample points to generate.
spacing : {"uniform", "cosine"}
Spacing rule applied over ``[0, 1]``.
Returns
-------
FloatArray
Monotone surface-local sample coordinates.
Raises
------
ValueError
If ``spacing`` is unsupported.
"""
if spacing == "uniform":
return as_float_array(
np.linspace(
0.0,
1.0,
num_points,
dtype=np.float64,
)
)
if spacing == "cosine":
theta = np.linspace(
0.0,
np.pi,
num_points,
dtype=np.float64,
)
return as_float_array(0.5 * (1.0 - np.cos(theta)))
msg = "spacing must be either 'uniform' or 'cosine'."
raise ValueError(msg)
#
# Inverse parameter queries
#
[docs]
def u_from_x(self, x: FloatInput, *, surface: AirfoilSurface) -> FloatArray:
"""
Return curve parameters that correspond to ``x``.
Parameters
----------
x : FloatInput
Chordwise coordinates in the normalized airfoil frame.
surface : {"lower", "upper"}
Surface to solve on.
Returns
-------
FloatArray
Curve parameters on the requested surface.
Raises
------
ValueError
If any requested chordwise coordinate lies outside the reachable
x-range of the selected surface.
"""
x_array = as_float_array(x)
xmin, xmax = self._surface_x_bounds(surface=surface)
if (x_array < xmin - DOMAIN_ABS_TOLERANCE).any() or (
x_array > xmax + DOMAIN_ABS_TOLERANCE
).any():
msg = (
"Invalid x-coordinate provided. "
f"The {surface} surface supports "
f"{xmin:.6g} <= x <= {xmax:.6g}."
)
raise ValueError(msg)
upper = surface_is_upper(surface)
def fun(u_i: FloatScalar, x_target: FloatScalar) -> FloatScalar:
x_value, _ = self.xy_from_u(u_i)
return as_float_scalar(x_value) - x_target
bracket = UPPER_SURFACE_BRACKET if upper else LOWER_SURFACE_BRACKET
u_array = np.empty_like(x_array)
flat_x = x_array.ravel()
flat_u = u_array.ravel()
for index, x_value in enumerate(flat_x):
x_target = as_float_scalar(x_value)
if x_target < X_MIN:
root = root_scalar(
fun,
args=(x_target,),
x0=U_LEADING_EDGE,
x1=abs(x_target),
xtol=ROOT_ABS_TOLERANCE,
rtol=ROOT_ABS_TOLERANCE,
maxiter=ROOT_MAX_ITERATION,
)
flat_u[index] = as_float_scalar(root.root)
elif np.abs(fun(bracket[0], x_target)) < ROOT_ABS_TOLERANCE:
flat_u[index] = bracket[0]
elif np.abs(fun(bracket[1], x_target)) < ROOT_ABS_TOLERANCE:
flat_u[index] = bracket[1]
else:
root = root_scalar(
fun,
args=(x_target,),
bracket=bracket,
xtol=ROOT_ABS_TOLERANCE,
rtol=ROOT_ABS_TOLERANCE,
maxiter=ROOT_MAX_ITERATION,
)
flat_u[index] = as_float_scalar(root.root)
return u_array
def _surface_x_bounds(
self, *, surface: AirfoilSurface
) -> tuple[FloatScalar, FloatScalar]:
"""
Return the reachable x-range for one airfoil surface.
Parameters
----------
surface : {"lower", "upper"}
Surface interval to query.
Returns
-------
FloatScalar
Minimum x-coordinate reachable on the selected surface.
FloatScalar
Maximum x-coordinate reachable on the selected surface.
"""
upper = surface_is_upper(surface)
bracket = UPPER_SURFACE_BRACKET if upper else LOWER_SURFACE_BRACKET
x0 = as_float_scalar(self.xy_from_u(bracket[0])[0])
x1 = as_float_scalar(self.xy_from_u(bracket[1])[0])
return min(x0, x1), max(x0, x1)
[docs]
@override
def u_from_s(self, s: FloatInput) -> FloatArray:
"""
Return curve parameters that correspond to arc length.
Parameters
----------
s : FloatInput
Arc lengths measured from the lower trailing edge.
Returns
-------
FloatArray
Curve parameters corresponding to ``s``.
Raises
------
ValueError
When arc-length provided is larger than airfoil surface length.
"""
s_array = as_float_array(s)
total_length = self.length
if (s_array > total_length).any() or (s_array < X_MIN).any():
msg = (
"Invalid arc length provided. "
f"Valid range is 0 <= s <= {total_length:.6g}."
)
raise ValueError(msg)
lower_surface_length = self._lower_surface_length()
u_array = np.empty_like(s_array)
lower_mask = s_array <= lower_surface_length
def solve_lower(s_target: FloatScalar) -> FloatScalar:
if np.abs(s_target) < ROOT_ABS_TOLERANCE:
return U_MIN
if np.abs(s_target - lower_surface_length) < ROOT_ABS_TOLERANCE:
return U_LEADING_EDGE
def lower_residual(u: FloatScalar) -> FloatScalar:
return as_float_scalar(self.arc_length(U_MIN, u)) - s_target
root = root_scalar(
lower_residual,
bracket=LOWER_SURFACE_BRACKET,
xtol=ROOT_ABS_TOLERANCE,
rtol=ROOT_ABS_TOLERANCE,
maxiter=ROOT_MAX_ITERATION,
)
return as_float_scalar(root.root)
def solve_upper(s_target: FloatScalar) -> FloatScalar:
upper_arc_length = s_target - lower_surface_length
if np.abs(upper_arc_length) < ROOT_ABS_TOLERANCE:
return U_LEADING_EDGE
if np.abs(s_target - total_length) < ROOT_ABS_TOLERANCE:
return U_MAX
def upper_residual(u: FloatScalar) -> FloatScalar:
return as_float_scalar(
as_float_scalar(self.arc_length(U_LEADING_EDGE, u))
- upper_arc_length
)
root = root_scalar(
upper_residual,
bracket=UPPER_SURFACE_BRACKET,
xtol=ROOT_ABS_TOLERANCE,
rtol=ROOT_ABS_TOLERANCE,
maxiter=ROOT_MAX_ITERATION,
)
return as_float_scalar(root.root)
flat_s = s_array.ravel()
flat_u = u_array.ravel()
flat_lower_mask = lower_mask.ravel()
for index, s_value in enumerate(flat_s):
s_target = as_float_scalar(s_value)
if flat_lower_mask[index]:
flat_u[index] = solve_lower(s_target)
else:
flat_u[index] = solve_upper(s_target)
return u_array
#
# Caching
#
def _airfoil_changed(self) -> None:
"""
Notify airfoil that shape has changed.
This needs to be called by child classes when the airfoil geometry
has changed so any cached values can be invalidated.
Returns
-------
None
This method updates internal caches in place.
"""
self._lower_surface_length_cache = None
self._upper_surface_length_cache = None