Metadata-Version: 2.5
Name: flapjax-model-gen
Version: 0.1.4
Summary: Match strip-wise aerodynamic corrections (cl0/cla/cm0/cma) to a NASTRAN-ready twist distribution using flapjax.
Author-email: Ben Preston <b.preston23@imperial.ac.uk>
License-Expression: MIT
License-File: LICENSE
Keywords: aeroelasticity,flapjax,nastran,vortex lattice method
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.12
Requires-Dist: flapjax>=1.3.6
Requires-Dist: jax>=0.8.1
Requires-Dist: pynastran>=1.3.4
Description-Content-Type: text/markdown

# flapjax-model-gen

Match strip-wise aerodynamic corrections (`cl0`, `cla`, `cm0`, `cma`) — defined about an aircraft's flight
shape, with no knowledge of local angle of attack or elastic deformation — to a NASTRAN-ready aerodynamic twist
distribution, using [flapjax](https://github.com/ben-l-p/flapjax) as the differentiable solver in the loop.

This is a thin, fast-moving companion package to `flapjax`, kept separate so it can be iterated on and released
independently.

## What's here

Most of the workflow this package supports is already implemented in `flapjax` itself
(`flapjax.aero.utils.apply_polar_correction` / `strip_alpha` / `project_forcing_to_beam`, wired into `UVLM` via
its `polar_data` / `polar_function` arguments — see `flapjax`'s `models/cantilever_wing/polar_correction.ipynb`
tutorial). This package adds the pieces that were missing:

- `linear_polar.py` — a `PolarFunction` (`flapjax`'s `(alpha, data) -> (cl, cd, cm)` interface) that evaluates
  an affine correction, `cl = cl0 + cla * alpha`, `cm = cm0 + cma * alpha`, `cd = 0`, so `cl0`/`cla`/`cm0`/`cma`
  data can be dropped straight into `UVLM(polar_data=..., polar_function=...)` exactly like a tabulated polar.
- `rigid.py` — rigidly rotates a structure's reference SE(3) frames (`rotate_hg`), for representing a change in
  angle of attack as a *rotation of the geometry*, matching the correction data's own convention, rather than
  rotating the freestream direction.
- `twist_grid.py` — `make_twisted_grid`, a per-spanwise-station generalisation of `flapjax`'s
  `make_rectangular_grid(..., twist=...)` (which applies one scalar twist to the whole grid). This is the search
  space the optimizer below works in.
- `strip_forces.py` — thin helpers around `AeroCase.project_forcing_to_beam` for pulling global-frame,
  per-node force/moment vectors out of a solved case, on either the rigid reference shape or a deformed one.
- `twist_match.py` — `match_twist`, a JAX-native Levenberg-Marquardt solver (autodiff Jacobian,
  `jax.lax.while_loop`) that finds the twist distribution minimising the residual between a candidate model's
  forces and a target force field. This is the genuinely new piece: no optimizer for this existed in `flapjax`.

## Suggested workflow

This mirrors the four-step process it was built for; wire the pieces together to fit your model rather than
treating this as a locked pipeline — in particular, **check the assumption flagged in step 3** before trusting
the result.

```python
import jax.numpy as jnp
from flapjax.aero import ConstantFlowField
from flapjax_model_gen import (
    LinearPolar, linear_polar_function, rotate_hg, rigid_corrected_forces,
    make_twisted_grid, match_twist,
)

# 1. Target forces directly from the correction data, on the rigid, undeformed flight shape.
#    `uvlm_corrected` is your UVLM/CoupledAeroelastic built with
#    polar_data=[LinearPolar(cl0, cla, cm0, cma)], polar_function=[linear_polar_function].
target_0 = rigid_corrected_forces(uvlm_corrected, hg=structure.hg0)
target_1 = rigid_corrected_forces(uvlm_corrected, hg=rotate_hg(structure.hg0, jnp.deg2rad(1.0)))

# 2. Static deformation under the corrected aero at AoA=0 (full FSI solve).
sol_0 = coupled_corrected.static_solve(prescribed_dofs=range(6), horseshoe=True)

# 3. Rigidly rotate the *converged, deformed* shape and take a single (no-FSI) aero pass at AoA=1.
#    NOTE: apply_polar_correction evaluates cl0/cla/cm0/cma at each strip's *locally computed* AoA, which in a
#    coupled solve includes any elastic twist -- extrapolating a correction that was only derived from rigid
#    sweeps. That's usually the intended reading of "sectional" data, but confirm it before relying on it.
hg_1 = rotate_hg(sol_0.structure.hg, jnp.deg2rad(1.0))
model_1 = rigid_corrected_forces(uvlm_corrected, hg=hg_1)

# 4. Build an *uncorrected* aero model with per-station twist as the free variable, and match it (on the same,
#    fixed deformed geometry from step 2/3) against whichever of the above you're treating as ground truth.
def forces_fn(twist):
    x0_aero = [make_twisted_grid(m, n, chord, ea, twist=twist)]
    uvlm_plain.set_design_variables(..., x0_aero=x0_aero)
    f0 = rigid_corrected_forces(uvlm_plain, hg=structure.hg0)
    f1 = rigid_corrected_forces(uvlm_plain, hg=rotate_hg(sol_0.structure.hg, jnp.deg2rad(1.0)))
    return jnp.stack([f0, f1])

target = jnp.stack([target_0, target_1])
result = match_twist(forces_fn, target, twist0=jnp.zeros(n + 1))
```

A single per-strip twist angle is a pure shift of local incidence: it can reproduce a `cl0`-like offset against
whatever lift-curve slope the plain panel aerodynamics already has, but it cannot independently fix a
mismatched `cla`, and it can't inject an independent camber-driven `cm0` (a flat panel has no camber). Worth
checking that the plain UVLM's native `cla` per strip is already close to the target before trusting the fit —
if it isn't, geometry/discretisation needs adjusting, not just twist.

## Building a flapjax model from NASTRAN inputs

`flapjax_model_gen.nastran` converts a NASTRAN aircraft model into flapjax `BeamStructure`/`UVLM` inputs, using
[pyNastran](https://github.com/SteveDoyle2/pyNastran) to parse bulk data. It expects three input files, matching
how these models are actually organised:

- **Aerodynamic mesh** (`.bdf`, `CAERO1` panels) — `nastran.parse_caero1_panels` extracts each panel's leading/
  trailing-edge corner points (already resolved to the global frame, including its own `CP` coordinate system);
  `nastran.build_local_grid` then rediscretizes them into a flapjax local grid. Pass `include=[...]` or
  `exclude=[...]` (CAERO element ids, mutually exclusive) to convert only a subset of surfaces from a
  multi-surface aero model, or to drop an unsupported (non-`CAERO1`) entry you don't need.
- **Mass model** (one or more `.nsb` files, `CONM2` entries only — e.g. one structural, one fuel) —
  `nastran.parse_mass_files` sums them into one lumped 6x6 mass matrix per node.
- **Stiffness model** (`.bdf`, `GRID` + `CBEAM` + `PBEAM`, referencing `MAT1`) — `nastran.parse_stiffness_bdf`
  builds node coordinates, connectivity, per-element orientation (`y_vector`) and per-element 6x6 stiffness.

```python
import jax.numpy as jnp
from flapjax.aero import GridDiscretisation, UVLM
from flapjax.structure import BeamStructure
from flapjax.coupled import CoupledAeroelastic
from flapjax_model_gen.nastran import (
    parse_stiffness_bdf, parse_mass_files, read_aero_bdf, parse_caero1_panels, build_local_grid,
)

beam_model = parse_stiffness_bdf("stiffness.bdf")
m_lumped, m_lumped_index = parse_mass_files(["structure.nsb", "fuel.nsb"], node_ids=beam_model.node_ids)

structure = BeamStructure(
    num_nodes=len(beam_model.node_ids),
    connectivity=beam_model.connectivity,
    y_vector=beam_model.y_vector,
    k_cs_index=beam_model.k_cs_index,
    m_lumped_index=m_lumped_index,
)

aero_model = read_aero_bdf("aero.bdf")
panels = parse_caero1_panels(aero_model)  # group/order multi-panel surfaces yourself, e.g. sorted by span
local_grid = build_local_grid(panels, beam_model.coords, m=8, span_axis=1)

uvlm = UVLM(
    grid_shapes=[GridDiscretisation(m=8, n=len(beam_model.node_ids) - 1, m_star=20)],
    dof_mapping=jnp.arange(len(beam_model.node_ids)),
)
wing = CoupledAeroelastic(structure, uvlm)
wing.set_design_variables(
    coords=beam_model.coords, k_cs=beam_model.k_cs, m_cs=None, m_lumped=m_lumped,
    dt=..., flowfield=..., x0_aero=[local_grid],
)
```

### Multiple aerodynamic surfaces

A model with more than one aerodynamic surface (wing, tail, fin, ...) needs `grid_shapes`/`x0_aero`/
`dof_mapping` as parallel *lists*, one entry per surface, each pointing at its own subset of the single,
shared `BeamStructure`'s nodes (`flapjax`'s `UVLM.dof_mapping` is `(n_surf, )(zeta_n, )`). Since a NASTRAN
aero model's own `SPLINE`/`SET1` cards already record exactly this "which CAERO panels + which structural
GRID ids form one continuous surface" grouping, `nastran.parse_spline_surfaces` reads it directly instead of
requiring it to be hand-listed, and `nastran.reorder_by_surfaces`/`nastran.reorder_nodes` reindex a
`NastranBeamModel`'s nodes (and elements) to match:

```python
from flapjax_model_gen.nastran import parse_spline_surfaces, reorder_by_surfaces

surfaces = parse_spline_surfaces(aero_model)  # groups CAERO ids sharing a SET1 into one SplineSurface each
beam_model = reorder_by_surfaces(
    beam_model, [s.node_ids for s in surfaces], span_axis=1
)  # nodes grouped surface-by-surface, root-to-tip within each; unspliced nodes appended at the end

grid_shapes, x0_aero, dof_mapping = [], [], []
index_of = {nid: i for i, nid in enumerate(beam_model.node_ids)}
for surface in surfaces:
    panels = sorted(
        parse_caero1_panels(aero_model, include=surface.caero_ids), key=lambda p: p.le_root[1]
    )
    node_ids = sorted(surface.node_ids, key=lambda nid: beam_model.coords[index_of[nid], 1])
    mapping = jnp.array([index_of[nid] for nid in node_ids])
    grid_shapes.append(GridDiscretisation(m=8, n=len(node_ids) - 1, m_star=20))
    x0_aero.append(build_local_grid(panels, beam_model.coords[mapping], m=8, span_axis=1))
    dof_mapping.append(mapping)

uvlm = UVLM(grid_shapes=grid_shapes, dof_mapping=dof_mapping)
```

`reorder_by_surfaces` only groups/orders the *structural* side; it's a plain `nastran.reorder_nodes(beam_model,
node_order)` under the hood, callable directly if you want to reorder by something other than a spline
grouping. `SPLINE1` and `SPLINE7` are supported (raises `NotImplementedError` for `SPLINE2`-`5`). `SPLINE6`
isn't implemented by pyNastran itself, and `SPLINE7` isn't either -- pyNastran drops it into
`model.reject_cards` rather than parsing it into an object -- so `parse_spline_surfaces` hand-parses
`SPLINE7`'s raw fields itself (per the MSC Nastran QRG's `SPLINE7 EID CAERO AELIST SETG ...` layout) and
explicitly detects+rejects `SPLINE6` rather than silently dropping that surface.

**What's deliberately out of scope / assumed**, each backed by a clear error or warning rather than a silent
wrong answer if violated:

- Only `CAERO1`, `CBEAM`/`PBEAM` (not `PBEAML`/`PBCOMP`), and `MAT1` are supported.
- `CBEAM`'s `OFFT` must be the default `'GGG'`, and end offsets (`WA`/`WB`) must be zero — no offset beam ends.
- `CONM2`'s `CID` must be the basic coordinate system (`0` or `-1`).
- A tapered `PBEAM` (multiple stations) is collapsed to one constant cross-section per element by averaging its
  station values — flapjax elements are constant-property, so this is an approximation for a genuinely tapered
  element, not a limitation you can configure around.
- `PBEAM`'s `I12` (bend-bend coupling) is not incorporated into `k_cs` — a warning fires if it's nonzero, since
  results will be approximate for that element.
- The aerodynamic grid's spanwise stations are placed exactly at the given structural node positions (flapjax's
  `dof_mapping` ties each aero column rigidly to one structural node, with no separate spline layer like
  NASTRAN's beam splines) — `build_local_grid` doesn't take an independent spanwise panel count.
- The axis convention (`CBEAM`'s orientation vector = flapjax's `y_vector` directly; PBEAM `I1`→bending about
  local z, `I2`→about local y) was verified against the NASTRAN QRG and cross-checked against flapjax's own
  beam code, not assumed — see the docstring in `nastran/stiffness.py` for the derivation.

## Development

```bash
uv sync --dev
uv run pytest
uv run ruff check src
```

Tests are deliberately lightweight (geometry/algebra invariants and a synthetic least-squares problem for the
optimizer) rather than full aeroelastic solves, so the suite stays fast to iterate against.
