"""General CST airfoil runtime type."""
from __future__ import annotations
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 (
BezierDemotionContinuity,
)
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.analytic import (
CstAirfoilSpec,
CstSurfaceSpec,
)
from .cst_geometry_side import CstGeometrySide
[docs]
class CstGeometry(Airfoil):
"""
General CST airfoil backed by upper and lower CST geometry sides.
Notes
-----
The runtime stores the exact schema content used to construct it so
supported CST geometry airfoils participate in the schema round-trip
contract.
This general runtime keeps the standard full-airfoil curve
parameterization ``x = |u|``.
When ``n1 < 1``, curve-parameter derivatives with respect to ``u`` can
remain singular at the leading edge because they are inherited directly
from the underlying ``dy/dx`` CST relation.
Canonical airfoil CST definitions with ``n1 = 0.5`` and ``n2 = 1.0`` use
:class:`CstAirfoil` instead, which applies a different curve
parameterization while preserving the same geometric shape.
"""
def __init__(
self,
*,
upper: CstGeometrySide,
lower: CstGeometrySide,
trailing_edge_thickness: FloatScalar = 0.0,
) -> None:
"""
Initialize a general CST airfoil from upper and lower sides.
Parameters
----------
upper : CstGeometrySide
Upper-side CST definition.
lower : CstGeometrySide
Lower-side CST definition.
trailing_edge_thickness : FloatScalar, default=0.0
Explicit trailing-edge gap as a fraction of chord.
"""
super().__init__()
half_thickness = 0.5 * trailing_edge_thickness
self._upper = CstGeometrySide(
shape=upper.shape,
n1=upper.n1,
n2=upper.n2,
delta_te=upper.delta_te + half_thickness,
)
self._lower = CstGeometrySide(
shape=lower.shape,
n1=lower.n1,
n2=lower.n2,
delta_te=lower.delta_te - half_thickness,
)
@property
def upper(self) -> CstGeometrySide:
"""
Return the upper-side CST geometry.
This property exposes the upper-side CST geometry definition.
"""
return self._upper
@property
def lower(self) -> CstGeometrySide:
"""
Return the lower-side CST geometry.
This property exposes the lower-side CST geometry definition.
"""
return self._lower
@property
def trailing_edge_thickness(self) -> FloatScalar:
"""
Return the explicit trailing-edge thickness.
This property reports the explicit trailing-edge gap as a fraction
of chord.
"""
return self.upper.delta_te - self.lower.delta_te
def _rebuild_with_side_shapes(
self,
*,
upper_shape: CstGeometrySide,
lower_shape: CstGeometrySide,
) -> CstGeometry:
"""Rebuild this CST airfoil from replacement side shapes."""
half_thickness = as_float_scalar(0.5 * self.trailing_edge_thickness)
upper = self.upper.rebuild_with_shape(
upper_shape.shape,
delta_te=as_float_scalar(self.upper.delta_te - half_thickness),
)
lower = self.lower.rebuild_with_shape(
lower_shape.shape,
delta_te=as_float_scalar(self.lower.delta_te + half_thickness),
)
return self.__class__(
upper=upper,
lower=lower,
trailing_edge_thickness=self.trailing_edge_thickness,
)
[docs]
def demote_degree(
self,
*,
count: int = 1,
continuity: BezierDemotionContinuity = "NOT_CONNECTED",
) -> CstGeometry:
"""
Lower the Bezier shape degree on both CST sides.
Parameters
----------
count : int, default=1
Number of Bezier degree-reduction steps applied to each side
shape curve.
continuity : {"NOT_CONNECTED", "C0", "C1", "C2"},
default="NOT_CONNECTED"
Symmetric endpoint continuity preserved during each side
demotion step.
Returns
-------
CstGeometry
Rebuilt CST airfoil with reduced-degree side shape curves.
Notes
-----
This operation is intentionally approximate unless the side
shape curves are exactly reducible to the requested lower
degree.
"""
return self._rebuild_with_side_shapes(
upper_shape=self.upper.demote_degree(
count=count,
continuity=continuity,
),
lower_shape=self.lower.demote_degree(
count=count,
continuity=continuity,
),
)
@property
def spec(self) -> CstAirfoilSpec:
"""
Return the schema definition used to create this airfoil.
The returned schema contains the current upper and lower side
coefficients, exponents, and trailing-edge thickness.
"""
return CstAirfoilSpec(
trailing_edge_thickness=as_float_scalar(
self.trailing_edge_thickness
),
upper=CstSurfaceSpec(
n1=as_float_scalar(self.upper.n1),
n2=as_float_scalar(self.upper.n2),
a=[as_float_scalar(value) for value in self.upper.coefficients],
),
lower=CstSurfaceSpec(
n1=as_float_scalar(self.lower.n1),
n2=as_float_scalar(self.lower.n2),
a=[as_float_scalar(value) for value in self.lower.coefficients],
),
)
[docs]
@override
def xy_from_u(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate the CST airfoil coordinates.
Parameters
----------
u : FloatInput
Signed airfoil parameter values in ``[-1, 1]``.
Negative values evaluate the lower surface and non-negative
values evaluate the upper surface.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(x, y)`` of ``float64`` arrays matching the normalized
shape of ``u``.
Notes
-----
This uses ``x = |u|`` on both surface branches.
"""
u_array = as_float_array(u)
x = as_float_array(np.abs(u_array))
y = np.empty_like(x)
upper_mask = u_array >= 0.0
lower_mask = ~upper_mask
if np.any(upper_mask):
upper_x = x[upper_mask]
y[upper_mask] = self.upper.y(upper_x)
if np.any(lower_mask):
lower_x = x[lower_mask]
y[lower_mask] = self.lower.y(lower_x)
return x, as_float_array(y)
[docs]
@override
def xy_u(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate first derivatives with respect to the airfoil parameter.
Parameters
----------
u : FloatInput
Signed airfoil parameter values in ``[-1, 1]``.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(dx/du, dy/du)`` of ``float64`` arrays.
Notes
-----
At listed breakpoints, this method returns the minus-side
derivative so array-valued evaluations remain single-valued.
For CST class exponents with singular ``dy/dx`` behavior at the
leading edge, this native derivative can remain singular.
"""
u_array = as_float_array(u)
x_u = np.empty_like(u_array)
np.sign(u_array, out=x_u)
y_u = np.empty_like(u_array)
analytic_mask = np.ones_like(u_array, dtype=bool)
breakpoints = self.breakpoints()
flat_u = u_array.ravel()
flat_x_u = x_u.ravel()
flat_y_u = y_u.ravel()
flat_mask = analytic_mask.ravel()
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):
analytic_u = u_array[analytic_mask]
analytic_x = np.abs(analytic_u)
upper_mask = analytic_u >= 0.0
lower_mask = ~upper_mask
analytic_y_u = np.empty_like(analytic_u)
if np.any(upper_mask):
upper_x = analytic_x[upper_mask]
upper_y_x = self.upper.y_x(upper_x)
analytic_y_u[upper_mask] = upper_y_x
if np.any(lower_mask):
lower_x = analytic_x[lower_mask]
lower_y_x = self.lower.y_x(lower_x)
analytic_y_u[lower_mask] = -lower_y_x
y_u[analytic_mask] = analytic_y_u
return as_float_array(x_u), as_float_array(y_u)
[docs]
@override
def xy_uu(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate second derivatives with respect to the airfoil parameter.
Parameters
----------
u : FloatInput
Signed airfoil parameter values in ``[-1, 1]``.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(d2x/du2, d2y/du2)`` of ``float64`` arrays.
Notes
-----
At listed breakpoints, this method returns the minus-side second
derivative so array-valued evaluations remain single-valued.
For CST class exponents with singular ``dy/dx`` or ``d2y/dx2``
behavior at the leading edge, this native derivative can remain
singular.
"""
u_array = as_float_array(u)
x_uu = np.zeros_like(u_array)
y_uu = np.empty_like(u_array)
analytic_mask = np.ones_like(u_array, dtype=bool)
breakpoints = self.breakpoints()
flat_u = u_array.ravel()
flat_x_uu = x_uu.ravel()
flat_y_uu = y_uu.ravel()
flat_mask = analytic_mask.ravel()
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):
analytic_u = u_array[analytic_mask]
analytic_x = np.abs(analytic_u)
upper_mask = analytic_u >= 0.0
lower_mask = ~upper_mask
analytic_y_uu = np.empty_like(analytic_u)
if np.any(upper_mask):
upper_x = analytic_x[upper_mask]
upper_y_xx = self.upper.y_xx(upper_x)
analytic_y_uu[upper_mask] = upper_y_xx
if np.any(lower_mask):
lower_x = analytic_x[lower_mask]
lower_y_xx = self.lower.y_xx(lower_x)
analytic_y_uu[lower_mask] = lower_y_xx
y_uu[analytic_mask] = analytic_y_uu
return as_float_array(x_uu), as_float_array(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 CST airfoil parameters matching ``xi`` on the
selected surface.
Notes
-----
General CST geometry 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 CST airfoil parameters in ``[-1, 1]``.
Returns
-------
SurfaceMappedValues
Surface-local ``xi`` values and upper-surface membership flags.
Notes
-----
General CST geometry uses the linear mapping ``xi = |u|``.
"""
return self._xi_from_u_signed_linear(u)
[docs]
@override
def breakpoints(self) -> list[FloatScalar]:
"""
Return the boundary and leading-edge breakpoints.
Returns
-------
list[float]
Ordered parameter locations where surface branches meet or
derivative one-sided limits may differ.
"""
return [-1.0, 0.0, 1.0]
@staticmethod
def _breakpoint_index(
value: FloatScalar,
*,
breakpoints: list[FloatScalar],
) -> int | None:
"""
Return the matching native breakpoint index within tolerance.
Parameters
----------
value : FloatScalar
Native CST 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 index.
Parameters
----------
index : int
Index into :meth:`breakpoints`.
Returns
-------
tuple[tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]]
``((x_u_minus, y_u_minus), (x_u_plus, y_u_plus))``.
Notes
-----
Endpoint breakpoints return the same exact boundary value for both
entries. The interior leading-edge breakpoint returns the lower and
upper side values explicitly.
"""
u_breakpoint = self.breakpoints()[index]
if u_breakpoint <= -1.0:
boundary_values = (
-1.0,
-as_float_scalar(self.lower.y_x(1.0)),
)
return boundary_values, boundary_values
if u_breakpoint >= 1.0:
boundary_values = (
1.0,
as_float_scalar(self.upper.y_x(1.0)),
)
return boundary_values, boundary_values
minus_values = (
-1.0,
as_float_scalar(self.lower.y_x(np.abs(u_breakpoint))),
)
plus_values = (
1.0,
as_float_scalar(self.upper.y_x(np.abs(u_breakpoint))),
)
return minus_values, plus_values
[docs]
def xy_uu_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided second derivatives at one breakpoint index.
Parameters
----------
index : int
Index into :meth:`breakpoints`.
Returns
-------
tuple[tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]]
``((x_uu_minus, y_uu_minus), (x_uu_plus, y_uu_plus))``.
Notes
-----
Endpoint breakpoints return the same exact boundary value for both
entries. The interior leading-edge breakpoint returns the lower and
upper side values explicitly.
"""
u_breakpoint = self.breakpoints()[index]
if u_breakpoint <= -1.0:
boundary_values = (
0.0,
as_float_scalar(self.lower.y_xx(1.0)),
)
return boundary_values, boundary_values
if u_breakpoint >= 1.0:
boundary_values = (
0.0,
as_float_scalar(self.upper.y_xx(1.0)),
)
return boundary_values, boundary_values
minus_values = (
0.0,
as_float_scalar(self.lower.y_xx(np.abs(u_breakpoint))),
)
plus_values = (
0.0,
as_float_scalar(self.upper.y_xx(np.abs(u_breakpoint))),
)
return minus_values, plus_values
[docs]
def xy_s_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided arc-length derivatives at one breakpoint index.
Notes
-----
This method composes the exact arc-length tangent values from the
exact native breakpoint derivatives returned by
:meth:`xy_u_breakpoint`.
"""
minus_u, plus_u = self.xy_u_breakpoint(index=index)
return (
self._xy_s_from_native_breakpoint(minus_u),
self._xy_s_from_native_breakpoint(plus_u),
)
[docs]
def xy_ss_breakpoint(
self,
*,
index: int,
) -> tuple[
tuple[FloatScalar, FloatScalar], tuple[FloatScalar, FloatScalar]
]:
"""
Return one-sided arc-length second derivatives at one 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`.
"""
minus_u, plus_u = self.xy_u_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),
)