Geometry API

flowdesign.geometry

Geometry helpers built around Kulfan/CST airfoil representations.

AirfoilCSTDataset

Bases: Dataset

Dataset for airfoils parameterized by a Kulfan/CST vector.

Design vector length = 2 * n_weights_per_side.

It reads .dat files, fits CST parameters, and caches them to a .pt file.

Source code in flowdesign/geometry/airfoil.py

class AirfoilCSTDataset(Dataset):
    """
    Dataset for airfoils parameterized by a Kulfan/CST vector.

    Design vector length = 2 * n_weights_per_side.

    It reads .dat files, fits CST parameters, and caches them to a .pt file.
    """

    def __init__(
        self,
        folder_path: str,
        cache_path: str = None,
        n_weights_per_side: int = 8,
        device: str = "cpu",
    ):
        self.folder_path = folder_path
        self.files = sorted([
            f for f in os.listdir(folder_path)
            if f.lower().endswith((".dat", ".csv"))
        ])
        self.cache_path = cache_path or os.path.join(folder_path, "cst_cache.pt")
        self.n_weights_per_side = n_weights_per_side
        self.device = device

        if os.path.exists(self.cache_path):
            print(f"[AirfoilCSTDataset] Loading cached CST data from {self.cache_path}")
            cached = torch.load(self.cache_path, map_location="cpu")

            if isinstance(cached, torch.Tensor):
                # If stored as a single stacked tensor
                self.data = [v for v in cached]
            else:
                # If stored as a list of tensors
                self.data = cached

        else:
            print("[AirfoilCSTDataset] Precomputing CST parameters...")
            self.data = []
            self._precompute_data()
            torch.save(self.data, self.cache_path)
            print(f"[AirfoilCSTDataset] Saved cache to {self.cache_path}")

    def _precompute_data(self):
        """
        Loop over all .dat files, fit CST parameters,
        and store the design vectors as tensors.
        """
        expected_len = 2 * self.n_weights_per_side

        for file in self.files:
            path = os.path.join(self.folder_path, file)
            try:
                geom = CSTGeometry.from_dat_file(
                    path,
                    n_weights_per_side=self.n_weights_per_side,
                )
                dv_tensor = geom.to_tensor(device="cpu")

                if dv_tensor.numel() != expected_len:
                    raise RuntimeError(
                        f"{file}: wrong DV length {dv_tensor.numel()}, "
                        f"expected {expected_len}"
                    )

                self.data.append(dv_tensor)

            except Exception as e:
                print(f"[AirfoilCSTDataset] Error processing {file}: {e}")
                continue

    def __len__(self):
        return len(self.data)

    def __getitem__(self, idx):
        """
        Return one design vector tensor on the target device.
        """
        return self.data[idx].to(self.device)

__getitem__(idx)

Return one design vector tensor on the target device.

Source code in flowdesign/geometry/airfoil.py

def __getitem__(self, idx):
    """
    Return one design vector tensor on the target device.
    """
    return self.data[idx].to(self.device)

CSTAirfoilGeometry

Problem-level CST/Kulfan airfoil parameterization.

Source code in flowdesign/geometry/airfoil.py

class CSTAirfoilGeometry:
    """Problem-level CST/Kulfan airfoil parameterization."""

    def __init__(
        self,
        n_upper: int = 8,
        n_lower: int = 8,
        prior_scale: float = 0.12,
        prior_mean: float = 0.0,
    ):
        self.n_upper = int(n_upper)
        self.n_lower = int(n_lower)
        self.design_dim = self.n_upper + self.n_lower
        self.prior_scale = float(prior_scale)
        self.prior_mean = prior_mean

    def sample_prior(self, n: int, device: str | torch.device | None = None) -> torch.Tensor:
        mean = torch.as_tensor(self.prior_mean, dtype=torch.float32, device=device)
        if mean.ndim == 0:
            mean = mean.expand(self.design_dim)
        if mean.numel() != self.design_dim:
            raise ValueError(f"prior_mean must be scalar or length {self.design_dim}, got {mean.numel()}")
        return mean.reshape(1, -1) + self.prior_scale * torch.randn(
            int(n), self.design_dim, dtype=torch.float32, device=device
        )

    def decode(self, params: torch.Tensor):
        params = self._as_batch(params)
        if get_kulfan_coordinates is None:
            return params
        coords = []
        for row in params.detach().cpu():
            geom = CSTGeometry(
                upper_weights=row[: self.n_upper].numpy(),
                lower_weights=row[self.n_upper :].numpy(),
            )
            coords.append(torch.tensor(geom.to_coordinates(), dtype=torch.float32))
        return torch.stack(coords)

    def encode(self, raw_geometry) -> torch.Tensor:
        if isinstance(raw_geometry, torch.Tensor):
            return raw_geometry.float()
        geom = CSTGeometry.from_coordinates(np.asarray(raw_geometry), n_weights_per_side=self.n_upper)
        return torch.tensor(geom.to_param_vector(), dtype=torch.float32)

    def validate(self, params: torch.Tensor) -> torch.Tensor:
        params = self._as_batch(params)
        if params.shape[1] != self.design_dim:
            raise ValueError(f"Expected CST params [B,{self.design_dim}], got {tuple(params.shape)}")
        return torch.isfinite(params).all(dim=1)

    def plot(self, params: torch.Tensor):
        import matplotlib.pyplot as plt

        decoded = self.decode(params)
        if isinstance(decoded, torch.Tensor) and decoded.ndim == 2:
            raise ImportError("Plotting CST coordinates requires aerosandbox.")
        coords = decoded[0].detach().cpu().numpy()
        fig, ax = plt.subplots()
        ax.plot(coords[:, 0], coords[:, 1])
        ax.axis("equal")
        return fig

    @staticmethod
    def _as_batch(params: torch.Tensor) -> torch.Tensor:
        params = torch.as_tensor(params, dtype=torch.float32)
        if params.ndim == 1:
            params = params.reshape(1, -1)
        if params.ndim == 3 and params.shape[1] == 1:
            params = params.reshape(params.shape[0], -1)
        return params

CSTGeometry

Generic NeuralFoil-compatible Kulfan (CST) geometry representation.

Design vector (exposed to models): [upper_weights(n), lower_weights(n)] -> length = 2 * n_weights_per_side

Internally, we also keep: leading_edge_weight, TE_thickness, N1, N2 for geometric reconstruction, but they are not required to be part of the DV.

Source code in flowdesign/geometry/airfoil.py

class CSTGeometry:
    """
    Generic NeuralFoil-compatible Kulfan (CST) geometry representation.

    Design vector (exposed to models):
        [upper_weights(n), lower_weights(n)]  -> length = 2 * n_weights_per_side

    Internally, we also keep:
        leading_edge_weight, TE_thickness, N1, N2
    for geometric reconstruction, but they are not required to be part of the DV.
    """

    def __init__(
        self,
        upper_weights: np.ndarray,
        lower_weights: np.ndarray,
        leading_edge_weight: float = 0.0,
        TE_thickness: float = 0.0,
        N1: float = 0.5,
        N2: float = 1.0,
    ):
        upper_weights = np.asarray(upper_weights, dtype=np.float32)
        lower_weights = np.asarray(lower_weights, dtype=np.float32)

        if upper_weights.shape != lower_weights.shape:
            raise ValueError("Upper and lower Kulfan weights must have the same dimension.")

        self.upper_weights = upper_weights
        self.lower_weights = lower_weights

        # Internal reconstruction parameters
        self.leading_edge_weight = float(leading_edge_weight)
        self.TE_thickness = float(TE_thickness)
        self.N1 = float(N1)
        self.N2 = float(N2)

    # -------------------------------------------------------
    # Constructors: from coordinates / from .dat file
    # -------------------------------------------------------

    @classmethod
    def from_coordinates(
        cls,
        coordinates: np.ndarray,
        n_weights_per_side: int = 8,
        N1: float = 0.5,
        N2: float = 1.0,
    ):
        """
        Fit Kulfan (CST) parameters from raw airfoil coordinates (N, 2).
        """
        if get_kulfan_parameters is None:
            raise ImportError("CST fitting requires aerosandbox.")
        coordinates = np.asarray(coordinates, dtype=float)

        params = get_kulfan_parameters(
            coordinates=coordinates,
            n_weights_per_side=n_weights_per_side,
            N1=N1,
            N2=N2,
            normalize_coordinates=True,
            use_leading_edge_modification=True,
            method="least_squares",
        )

        return cls(
            upper_weights=params["upper_weights"],
            lower_weights=params["lower_weights"],
            leading_edge_weight=params.get("leading_edge_weight", 0.0),
            TE_thickness=params.get("TE_thickness", 0.0),
            N1=params.get("N1", N1),
            N2=params.get("N2", N2),
        )

    @classmethod
    def from_dat_file(
        cls,
        dat_path: str,
        n_weights_per_side: int = 8,
        N1: float = 0.5,
        N2: float = 1.0,
    ):
        """
        Build CSTGeometry from a .dat/.csv file with airfoil coordinates.
        """
        try:
            coords = np.loadtxt(dat_path)
        except ValueError:
            # Handle possible header line
            coords = np.loadtxt(dat_path, skiprows=1)

        coords = np.asarray(coords, dtype=float)
        if coords.ndim != 2 or coords.shape[1] != 2:
            raise ValueError(f"{dat_path} is not a valid (N, 2) airfoil coordinate file.")

        return cls.from_coordinates(
            coordinates=coords,
            n_weights_per_side=n_weights_per_side,
            N1=N1,
            N2=N2,
        )

    # -------------------------------------------------------
    # Export: coordinates / param vector / tensor
    # -------------------------------------------------------

    def to_coordinates(self, n_points_per_side: int = 200) -> np.ndarray:
        """
        Reconstruct airfoil surface coordinates from stored Kulfan parameters.
        """
        if get_kulfan_coordinates is None:
            raise ImportError("CST coordinate reconstruction requires aerosandbox.")
        coords = get_kulfan_coordinates(
            lower_weights=self.lower_weights,
            upper_weights=self.upper_weights,
            leading_edge_weight=self.leading_edge_weight,
            TE_thickness=self.TE_thickness,
            n_points_per_side=n_points_per_side,
            N1=self.N1,
            N2=self.N2,
        )
        return np.asarray(coords, dtype=np.float32)

    def to_param_vector(self) -> np.ndarray:
        """
        Return design vector: [upper_weights, lower_weights].
        Length = 2 * n_weights_per_side.
        """
        return np.concatenate(
            [self.upper_weights, self.lower_weights],
            axis=0
        ).astype(np.float32)

    def to_tensor(self, device: str = "cpu") -> torch.Tensor:
        """
        Convert design vector to a PyTorch tensor on the given device.
        """
        return torch.tensor(self.to_param_vector(), dtype=torch.float32, device=device)

    def __repr__(self):
        n = self.upper_weights.shape[0]
        return f"CSTGeometry(n_weights_per_side={n})"

from_coordinates(coordinates, n_weights_per_side=8, N1=0.5, N2=1.0) classmethod

Fit Kulfan (CST) parameters from raw airfoil coordinates (N, 2).

Source code in flowdesign/geometry/airfoil.py

@classmethod
def from_coordinates(
    cls,
    coordinates: np.ndarray,
    n_weights_per_side: int = 8,
    N1: float = 0.5,
    N2: float = 1.0,
):
    """
    Fit Kulfan (CST) parameters from raw airfoil coordinates (N, 2).
    """
    if get_kulfan_parameters is None:
        raise ImportError("CST fitting requires aerosandbox.")
    coordinates = np.asarray(coordinates, dtype=float)

    params = get_kulfan_parameters(
        coordinates=coordinates,
        n_weights_per_side=n_weights_per_side,
        N1=N1,
        N2=N2,
        normalize_coordinates=True,
        use_leading_edge_modification=True,
        method="least_squares",
    )

    return cls(
        upper_weights=params["upper_weights"],
        lower_weights=params["lower_weights"],
        leading_edge_weight=params.get("leading_edge_weight", 0.0),
        TE_thickness=params.get("TE_thickness", 0.0),
        N1=params.get("N1", N1),
        N2=params.get("N2", N2),
    )

from_dat_file(dat_path, n_weights_per_side=8, N1=0.5, N2=1.0) classmethod

Build CSTGeometry from a .dat/.csv file with airfoil coordinates.

Source code in flowdesign/geometry/airfoil.py

@classmethod
def from_dat_file(
    cls,
    dat_path: str,
    n_weights_per_side: int = 8,
    N1: float = 0.5,
    N2: float = 1.0,
):
    """
    Build CSTGeometry from a .dat/.csv file with airfoil coordinates.
    """
    try:
        coords = np.loadtxt(dat_path)
    except ValueError:
        # Handle possible header line
        coords = np.loadtxt(dat_path, skiprows=1)

    coords = np.asarray(coords, dtype=float)
    if coords.ndim != 2 or coords.shape[1] != 2:
        raise ValueError(f"{dat_path} is not a valid (N, 2) airfoil coordinate file.")

    return cls.from_coordinates(
        coordinates=coords,
        n_weights_per_side=n_weights_per_side,
        N1=N1,
        N2=N2,
    )

to_coordinates(n_points_per_side=200)

Reconstruct airfoil surface coordinates from stored Kulfan parameters.

Source code in flowdesign/geometry/airfoil.py

def to_coordinates(self, n_points_per_side: int = 200) -> np.ndarray:
    """
    Reconstruct airfoil surface coordinates from stored Kulfan parameters.
    """
    if get_kulfan_coordinates is None:
        raise ImportError("CST coordinate reconstruction requires aerosandbox.")
    coords = get_kulfan_coordinates(
        lower_weights=self.lower_weights,
        upper_weights=self.upper_weights,
        leading_edge_weight=self.leading_edge_weight,
        TE_thickness=self.TE_thickness,
        n_points_per_side=n_points_per_side,
        N1=self.N1,
        N2=self.N2,
    )
    return np.asarray(coords, dtype=np.float32)

to_param_vector()

Return design vector: [upper_weights, lower_weights]. Length = 2 * n_weights_per_side.

Source code in flowdesign/geometry/airfoil.py

def to_param_vector(self) -> np.ndarray:
    """
    Return design vector: [upper_weights, lower_weights].
    Length = 2 * n_weights_per_side.
    """
    return np.concatenate(
        [self.upper_weights, self.lower_weights],
        axis=0
    ).astype(np.float32)

to_tensor(device='cpu')

Convert design vector to a PyTorch tensor on the given device.

Source code in flowdesign/geometry/airfoil.py

def to_tensor(self, device: str = "cpu") -> torch.Tensor:
    """
    Convert design vector to a PyTorch tensor on the given device.
    """
    return torch.tensor(self.to_param_vector(), dtype=torch.float32, device=device)

SuperWingConfigDataset

Bases: Dataset

Dataset for SuperWing geometry parameters from configs.dat.

SuperWing's current configs.dat contains one integer shape index followed by numeric columns. For geometry-only pretraining, the default behavior keeps the first 38 numeric columns after the index:

7 global planform + 10 spanwise airfoil/twist + 21 baseline airfoil terms

The remaining columns encode operating conditions and are intentionally not part of the generated design vector in this MVP.

Source code in flowdesign/geometry/superwing.py

class SuperWingConfigDataset(Dataset):
    """
    Dataset for SuperWing geometry parameters from configs.dat.

    SuperWing's current configs.dat contains one integer shape index followed by
    numeric columns. For geometry-only pretraining, the default behavior keeps
    the first 38 numeric columns after the index:

        7 global planform + 10 spanwise airfoil/twist + 21 baseline airfoil terms

    The remaining columns encode operating conditions and are intentionally not
    part of the generated design vector in this MVP.
    """

    def __init__(
        self,
        configs_path: str,
        design_dim: int = 38,
        device: str | torch.device = "cpu",
        download: bool = False,
    ):
        if download:
            download_superwing_configs(configs_path)
        if not os.path.exists(configs_path):
            raise FileNotFoundError(
                f"SuperWing configs.dat not found at {configs_path}. "
                "Pass --download-configs or provide --configs-path."
            )

        self.configs_path = configs_path
        self.design_dim = int(design_dim)
        self.device = torch.device(device)

        rows: list[list[float]] = []
        with open(configs_path, "r", encoding="utf-8") as handle:
            for line in handle:
                stripped = line.strip()
                if not stripped or stripped.startswith("VARIABLES"):
                    continue
                values = [float(token) for token in stripped.split()]
                # Drop the shape index in column 0.
                numeric = values[1:]
                if len(numeric) < self.design_dim:
                    raise ValueError(
                        f"Requested design_dim={self.design_dim}, but row has only "
                        f"{len(numeric)} numeric columns after the shape index."
                    )
                rows.append(numeric[: self.design_dim])

        if not rows:
            raise ValueError(f"No SuperWing rows parsed from {configs_path}")

        self.data = torch.tensor(rows, dtype=torch.float32)

    def __len__(self) -> int:
        return self.data.shape[0]

    def __getitem__(self, idx: int) -> torch.Tensor:
        return self.data[idx].to(self.device)

SuperWingParameterGeometry

Problem-level SuperWing parameter geometry.

Source code in flowdesign/geometry/superwing.py

class SuperWingParameterGeometry:
    """Problem-level SuperWing parameter geometry."""

    def __init__(
        self,
        design_dim: int = 38,
        configs_path: str | None = None,
        prior_scale: float = 0.05,
        device: str | torch.device = "cpu",
    ):
        self.design_dim = int(design_dim)
        self.configs_path = configs_path
        self.prior_scale = float(prior_scale)
        self.device = torch.device(device)
        self.builder = TorchSuperWingGeometryBuilder()
        self._dataset: SuperWingConfigDataset | None = None
        if configs_path and os.path.exists(configs_path):
            self._dataset = SuperWingConfigDataset(configs_path, design_dim=self.design_dim, device=self.device)

    def sample_prior(self, n: int, device: str | torch.device | None = None) -> torch.Tensor:
        device = torch.device(device) if device is not None else self.device
        if self._dataset is not None:
            idx = torch.randint(0, len(self._dataset), (int(n),))
            base = torch.stack([self._dataset[int(i)].cpu() for i in idx]).to(device)
            return base + self.prior_scale * torch.randn_like(base)
        return self._default_prior(int(n), device)

    def decode(self, params: torch.Tensor):
        return self.builder.to(params.device)(params)

    def encode(self, raw_geometry) -> torch.Tensor:
        return torch.as_tensor(raw_geometry, dtype=torch.float32, device=self.device)

    def validate(self, params: torch.Tensor) -> torch.Tensor:
        params = self._as_batch(params)
        return torch.isfinite(params).all(dim=1) & (params.shape[1] >= self.design_dim)

    def plot(self, params: torch.Tensor):
        import matplotlib.pyplot as plt

        geom = self.decode(self._as_batch(params[:1])).origingeom[0].detach().cpu()
        fig = plt.figure()
        ax = fig.add_subplot(111, projection="3d")
        surface = geom.permute(1, 2, 0)
        ax.plot_surface(surface[:, :, 0], surface[:, :, 1], surface[:, :, 2], linewidth=0, antialiased=False)
        return fig

    def _default_prior(self, n: int, device: torch.device) -> torch.Tensor:
        base = torch.zeros(n, self.design_dim, dtype=torch.float32, device=device)
        defaults = [
            25.0, 3.0, 1.5, 9.0, 0.35, 0.45, 0.0,
            0.12, 0.10, 0.08, 0.0, 0.0, 0.0, 2.0, 1.0, 0.0, -1.0,
        ]
        base[:, : len(defaults)] = torch.tensor(defaults, dtype=base.dtype, device=device)
        if self.design_dim >= 38:
            base[:, 17:27] = torch.linspace(0.15, 0.05, 10, device=device)
            base[:, 27:37] = -torch.linspace(0.12, 0.03, 10, device=device)
            base[:, 37] = 0.12
        return base + self.prior_scale * torch.randn_like(base)

    @staticmethod
    def _as_batch(params: torch.Tensor) -> torch.Tensor:
        params = torch.as_tensor(params, dtype=torch.float32)
        if params.ndim == 1:
            params = params.reshape(1, -1)
        if params.ndim == 3 and params.shape[1] == 1:
            params = params.reshape(params.shape[0], -1)
        return params

SuperWingGeometry dataclass

Reference meshes reconstructed from SuperWing geometry parameters.

Source code in flowdesign/geometry/superwing_builder.py

@dataclass
class SuperWingGeometry:
    """Reference meshes reconstructed from SuperWing geometry parameters."""

    origingeom: torch.Tensor
    geom0: torch.Tensor
    ref_area: torch.Tensor
    half_span: torch.Tensor
    ref_chord: torch.Tensor
    xref: torch.Tensor
    yref: torch.Tensor

TorchSuperWingGeometryBuilder

Bases: Module

Differentiable SuperWing reference-mesh builder.

Input rows follow configs.dat after removing the leading shape index: SA, DA, DA_kink, AR, TR, kink, rootadj, rtcs(3), cambers(3), twists(4), cst_u(10), cst_l(10), tcroot.

The generated meshes follow the AeroTransformer/SuperWing reference format: origingeom is a vertex mesh [B, 3, 129, 257]; geom0 is the cell-center mesh [B, 3, 128, 256] used by ATsurf models. This implementation ports the public cfdpost.wing.multi_section_wing.MultiSecWing meshing equations into torch, with a differentiable cubic interpolation approximation for spanwise distributions.

Source code in flowdesign/geometry/superwing_builder.py

class TorchSuperWingGeometryBuilder(nn.Module):
    """
    Differentiable SuperWing reference-mesh builder.

    Input rows follow `configs.dat` after removing the leading shape index:
    `SA, DA, DA_kink, AR, TR, kink, rootadj, rtcs(3), cambers(3),
    twists(4), cst_u(10), cst_l(10), tcroot`.

    The generated meshes follow the AeroTransformer/SuperWing reference format:
    `origingeom` is a vertex mesh `[B, 3, 129, 257]`; `geom0` is the cell-center
    mesh `[B, 3, 128, 256]` used by ATsurf models. This implementation ports the
    public `cfdpost.wing.multi_section_wing.MultiSecWing` meshing equations into
    torch, with a differentiable cubic interpolation approximation for spanwise
    distributions.
    """

    def __init__(
        self,
        nx: int = 129,
        nz: int = 129,
        tail: float = 0.008,
        cabin_ratio: float = 0.1,
        xref: float = 0.25,
        yref: float = 0.0,
        clamp_inputs: bool = True,
    ) -> None:
        super().__init__()
        self.nx = int(nx)
        self.nz = int(nz)
        self.tail = float(tail)
        self.cabin_ratio = float(cabin_ratio)
        self.clamp_inputs = bool(clamp_inputs)

        x_dist = self._dist_clustcos(self.nx)
        eta = torch.linspace(self.cabin_ratio, 1.0, self.nz)
        basis = self._cst_basis(x_dist, n_coef=10)

        self.register_buffer("x_dist", x_dist)
        self.register_buffer("eta", eta)
        self.register_buffer("cst_basis", basis)
        self.register_buffer("ref_chord_value", torch.tensor(1.0))
        self.register_buffer("xref_value", torch.tensor(float(xref)))
        self.register_buffer("yref_value", torch.tensor(float(yref)))

    @staticmethod
    def _dist_clustcos(nn: int, a0: float = 0.0079, a1: float = 0.96) -> torch.Tensor:
        yt = torch.linspace(0.0, 1.0, nn)
        aa = (1.0 - math.cos(a0 * math.pi)) / 2.0
        dd = (1.0 - math.cos(a1 * math.pi)) / 2.0 - aa
        angle = math.pi * (a0 * (1.0 - yt) + a1 * yt)
        return ((1.0 - torch.cos(angle)) / 2.0 - aa) / dd

    @staticmethod
    def _cst_basis(x: torch.Tensor, n_coef: int) -> torch.Tensor:
        rows = []
        class_fn = torch.sqrt(torch.clamp(x, min=0.0)) * (1.0 - x)
        for i in range(n_coef):
            coeff = math.comb(n_coef - 1, i)
            rows.append(coeff * x.pow(i) * (1.0 - x).pow(n_coef - 1 - i) * class_fn)
        basis = torch.stack(rows, dim=0)
        basis[:, 0] = 0.0
        basis[:, -1] = 0.0
        return basis

    def _sanitize(self, x: torch.Tensor) -> torch.Tensor:
        if not self.clamp_inputs:
            return x
        cols = [x[:, i] for i in range(x.shape[1])]
        cols[3] = cols[3].clamp_min(1.0)       # AR
        cols[4] = cols[4].clamp_min(0.02)      # TR
        cols[5] = cols[5].clamp(0.101, 0.98)   # kink, avoid singular inner span
        for idx in range(7, 10):
            cols[idx] = cols[idx].clamp_min(0.05)
        cols[37] = cols[37].clamp_min(0.02)    # tcroot
        return torch.stack(cols, dim=1)

    def _cst_foil(self, cst_u: torch.Tensor, cst_l: torch.Tensor, tcroot: torch.Tensor) -> tuple[torch.Tensor, torch.Tensor]:
        yu = cst_u @ self.cst_basis.to(dtype=cst_u.dtype, device=cst_u.device)
        yl = cst_l @ self.cst_basis.to(dtype=cst_l.dtype, device=cst_l.device)

        thick = yu - yl
        max_thick, idx = thick.max(dim=1)
        x = self.x_dist.to(dtype=cst_u.dtype, device=cst_u.device)
        x_at_max = x.gather(0, idx)
        scale = (tcroot - self.tail * x_at_max) / (max_thick + 1e-8)
        yu = yu * scale[:, None]
        yl = yl * scale[:, None]

        tail_term = 0.5 * self.tail * x[None, :]
        return yu + tail_term, yl - tail_term

    @staticmethod
    def _natural_cubic_eval(xk: torch.Tensor, yk: torch.Tensor, xq: torch.Tensor) -> torch.Tensor:
        """Evaluate a batched natural cubic spline at batched query points."""
        batch, n_knots = xk.shape
        if n_knots == 2:
            denom = (xk[:, 1:2] - xk[:, 0:1]).clamp_min(1e-8)
            weight = (xq - xk[:, 0:1]) / denom
            return yk[:, 0:1] + weight * (yk[:, 1:2] - yk[:, 0:1])

        h = (xk[:, 1:] - xk[:, :-1]).clamp_min(1e-8)
        rhs_inner = 6.0 * (
            (yk[:, 2:] - yk[:, 1:-1]) / h[:, 1:]
            - (yk[:, 1:-1] - yk[:, :-2]) / h[:, :-1]
        )
        n_inner = n_knots - 2
        lower = h[:, :n_inner]
        diag = 2.0 * (h[:, :n_inner] + h[:, 1 : n_inner + 1])
        upper = h[:, 1 : n_inner + 1]

        c_prime: list[torch.Tensor] = []
        d_prime: list[torch.Tensor] = []
        denom = diag[:, 0].clamp_min(1e-8)
        c_prime.append(upper[:, 0] / denom if n_inner > 1 else torch.zeros_like(denom))
        d_prime.append(rhs_inner[:, 0] / denom)
        for i in range(1, n_inner):
            denom = (diag[:, i] - lower[:, i] * c_prime[i - 1]).clamp_min(1e-8)
            c_prime.append(upper[:, i] / denom if i < n_inner - 1 else torch.zeros_like(denom))
            d_prime.append((rhs_inner[:, i] - lower[:, i] * d_prime[i - 1]) / denom)

        inner: list[torch.Tensor] = [d_prime[-1]]
        for i in range(n_inner - 2, -1, -1):
            inner.append(d_prime[i] - c_prime[i] * inner[-1])
        inner = list(reversed(inner))
        second = torch.cat(
            [
                torch.zeros(batch, 1, dtype=xk.dtype, device=xk.device),
                torch.stack(inner, dim=1),
                torch.zeros(batch, 1, dtype=xk.dtype, device=xk.device),
            ],
            dim=1,
        )

        out = torch.zeros(xq.shape, dtype=xq.dtype, device=xq.device)
        assigned = torch.zeros(xq.shape, dtype=torch.bool, device=xq.device)
        for i in range(n_knots - 1):
            left = xk[:, i : i + 1]
            right = xk[:, i + 1 : i + 2]
            hi = (right - left).clamp_min(1e-8)
            a = (right - xq) / hi
            b = (xq - left) / hi
            val = (
                a * yk[:, i : i + 1]
                + b * yk[:, i + 1 : i + 2]
                + ((a.pow(3) - a) * second[:, i : i + 1] + (b.pow(3) - b) * second[:, i + 1 : i + 2])
                * hi.pow(2)
                / 6.0
            )
            if i == 0:
                mask = xq <= right
            elif i == n_knots - 2:
                mask = xq >= left
            else:
                mask = (xq >= left) & (xq <= right)
            mask = mask & (~assigned)
            out = torch.where(mask, val, out)
            assigned = assigned | mask
        return out

    @staticmethod
    def _outer_dihedral(y: torch.Tensor, kink_y: torch.Tensor, tip_y: torch.Tensor, z_kink: torch.Tensor, z_tip: torch.Tensor, slope_kink: torch.Tensor) -> torch.Tensor:
        span = (tip_y - kink_y).clamp_min(1e-8)
        s = (y - kink_y) / span
        s = s.clamp(0.0, 1.0)
        slope_tip = (z_tip - z_kink) / span
        h00 = 2.0 * s.pow(3) - 3.0 * s.pow(2) + 1.0
        h10 = s.pow(3) - 2.0 * s.pow(2) + s
        h01 = -2.0 * s.pow(3) + 3.0 * s.pow(2)
        h11 = s.pow(3) - s.pow(2)
        return h00 * z_kink + h10 * span * slope_kink + h01 * z_tip + h11 * span * slope_tip

    def forward(self, design: torch.Tensor) -> SuperWingGeometry:
        if design.dim() == 3 and design.shape[1] == 1:
            design = design.squeeze(1)
        if design.dim() == 1:
            design = design.reshape(1, -1)
        if design.dim() != 2 or design.shape[1] < 38:
            raise ValueError(f"Expected SuperWing geometry design [B,>=38], got {tuple(design.shape)}")

        x = self._sanitize(design[:, :38].float())
        dtype, device = x.dtype, x.device
        batch = x.shape[0]

        sa, da, da_kink = x[:, 0], x[:, 1], x[:, 2]
        ar, tr, kink, rootadj = x[:, 3], x[:, 4], x[:, 5], x[:, 6]
        rtcs = x[:, 7:10]
        cambers = x[:, 10:13]
        twists = x[:, 13:17]
        cst_u = x[:, 17:27]
        cst_l = x[:, 27:37]
        tcroot = x[:, 37]

        eta = self.eta.to(dtype=dtype, device=device)[None, :]
        foil_x = self.x_dist.to(dtype=dtype, device=device)
        half_span = 0.25 * ar * (1.0 + tr)
        y = half_span[:, None] * eta
        cabin_y = half_span[:, None] * self.cabin_ratio
        kink_y = half_span[:, None] * kink[:, None]
        tip_y = half_span[:, None] + 1e-9

        deg = math.pi / 180.0
        sa_tan = torch.tan(sa[:, None] * deg)
        da_tan = torch.tan(da[:, None] * deg)
        da_kink_tan = torch.tan(da_kink[:, None] * deg)

        x_le = sa_tan * y
        x_le_kink = sa_tan * kink_y
        z_inner = da_kink_tan * y
        z_kink = da_kink_tan * kink_y
        z_tip = da_tan * half_span[:, None]
        z_outer = self._outer_dihedral(y, kink_y, tip_y, z_kink, z_tip, da_kink_tan)
        z_le = torch.where(eta <= kink[:, None], z_inner, z_outer)

        chord_tip = tr[:, None]
        chord_kink = tr[:, None] * kink[:, None] + (1.0 - kink[:, None])
        root_adj_l = rootadj[:, None] * (x_le_kink + chord_kink - 1.0)
        chord_cabin = (
            tr[:, None] * self.cabin_ratio
            + (1.0 - self.cabin_ratio)
            + root_adj_l * (kink[:, None] - self.cabin_ratio) / kink[:, None].clamp_min(1e-8)
        )
        chord_inner = chord_cabin + (chord_kink - chord_cabin) * (y - cabin_y) / (kink_y - cabin_y).clamp_min(1e-8)
        chord_outer = chord_kink + (chord_tip - chord_kink) * (y - kink_y) / (half_span[:, None] - kink_y).clamp_min(1e-8)
        chords = torch.where(eta < kink[:, None], chord_inner, chord_outer)
        ref_area = (
            (chord_cabin[:, 0] + chord_kink[:, 0]) * (kink_y[:, 0] - cabin_y[:, 0]) * 0.5
            + (chord_kink[:, 0] + chord_tip[:, 0]) * (half_span - kink_y[:, 0]) * 0.5
        ).clamp_min(1e-6)

        yu, yl = self._cst_foil(cst_u, cst_l, tcroot)
        camber_root = 0.5 * (yu + yl)
        thickness_root = yu - yl

        mid_outer = 0.5 * (kink_y + half_span[:, None])
        mid_inner = 0.5 * kink_y
        rtcs_knots = torch.stack(
            [
                torch.zeros_like(half_span),
                kink_y[:, 0],
                mid_outer[:, 0],
                tip_y[:, 0],
            ],
            dim=1,
        )
        rtcs_values = torch.stack(
            [
                torch.ones_like(half_span),
                rtcs[:, 0],
                rtcs[:, 0] * rtcs[:, 1],
                rtcs[:, 0] * rtcs[:, 1] * rtcs[:, 2],
            ],
            dim=1,
        )
        thickness_scale = self._natural_cubic_eval(rtcs_knots, rtcs_values, y)

        camber_knots = torch.stack(
            [
                torch.zeros_like(half_span),
                mid_inner[:, 0],
                kink_y[:, 0],
                mid_outer[:, 0],
                tip_y[:, 0],
            ],
            dim=1,
        )
        camber_values = torch.stack(
            [
                torch.zeros_like(half_span),
                cambers[:, 0] * cambers[:, 1],
                cambers[:, 1],
                torch.ones_like(half_span),
                cambers[:, 2],
            ],
            dim=1,
        )
        camber_scale = self._natural_cubic_eval(camber_knots, camber_values, y)

        y_upper = camber_root[:, :, None] * camber_scale[:, None, :] + 0.5 * thickness_root[:, :, None] * thickness_scale[:, None, :]
        y_lower = camber_root[:, :, None] * camber_scale[:, None, :] - 0.5 * thickness_root[:, :, None] * thickness_scale[:, None, :]
        section_z = torch.cat([torch.flip(y_lower[:, 1:, :], dims=[1]), y_upper], dim=1)
        section_x = torch.cat([torch.flip(foil_x[1:], dims=[0]), foil_x], dim=0)

        root_twist = torch.zeros(batch, 1, dtype=dtype, device=device)
        twist_angles = torch.cat([root_twist, twists], dim=1).cumsum(dim=1)
        twist_scale = self._natural_cubic_eval(camber_knots, twist_angles, y)
        twist_rad = twist_scale * deg

        z_le = z_le + 0.25 * chords * torch.sin(twist_rad)
        x_section = section_x[None, :, None]
        cos_twist = torch.cos(twist_rad)[:, None, :]
        sin_twist = torch.sin(twist_rad)[:, None, :]
        x_rot = cos_twist * x_section + sin_twist * section_z
        z_rot = -sin_twist * x_section + cos_twist * section_z

        x_coord = x_rot * chords[:, None, :] + x_le[:, None, :]
        y_coord = z_rot * chords[:, None, :] + z_le[:, None, :]
        z_coord = y[:, None, :].expand(batch, 2 * self.nx - 1, self.nz)

        block = torch.stack([x_coord, y_coord, z_coord], dim=1).transpose(2, 3)
        geom0 = 0.25 * (
            block[:, :, 1:, 1:]
            + block[:, :, 1:, :-1]
            + block[:, :, :-1, 1:]
            + block[:, :, :-1, :-1]
        )

        return SuperWingGeometry(
            origingeom=block,
            geom0=geom0,
            ref_area=ref_area,
            half_span=half_span,
            ref_chord=self.ref_chord_value.to(dtype=dtype, device=device).expand(batch),
            xref=self.xref_value.to(dtype=dtype, device=device).expand(batch),
            yref=self.yref_value.to(dtype=dtype, device=device).expand(batch),
        )

download_superwing_configs(path, url=SUPERWING_CONFIGS_URL)

Download SuperWing configs.dat if it is not already present.

The file is small enough for an MVP workflow and avoids downloading the multi-TB field data.

Source code in flowdesign/geometry/superwing.py

def download_superwing_configs(path: str, url: str = SUPERWING_CONFIGS_URL) -> str:
    """
    Download SuperWing configs.dat if it is not already present.

    The file is small enough for an MVP workflow and avoids downloading the
    multi-TB field data.
    """
    if os.path.exists(path):
        return path

    os.makedirs(os.path.dirname(os.path.abspath(path)), exist_ok=True)
    urllib.request.urlretrieve(url, path)
    return path