Metadata-Version: 2.5
Name: flapjax-model-gen
Version: 0.1.29
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: matplotlib>=3.8.4
Requires-Dist: pynastran>=1.3.4
Description-Content-Type: text/markdown

# flapjax-model-gen

Generate a [flapjax](https://github.com/ben-l-p/flapjax) aeroelastic model — structure, mass, aerodynamic mesh,
control surfaces, half-models — directly from NASTRAN input files, using
[pyNastran](https://github.com/SteveDoyle2/pyNastran) to parse bulk data. Also ships a companion workflow for
matching strip-wise aerodynamic corrections to a NASTRAN-ready twist distribution.

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

## Building a flapjax model from NASTRAN inputs

`flapjax_model_gen.nastran` converts a NASTRAN aircraft model into a fully populated flapjax
`CoupledAeroelastic`, ready for `.reference_configuration()`/`.static_solve()`. It expects three input files,
matching how these models are actually organised:

- **Stiffness model** (`.bdf`, `GRID` + `CBEAM` + `PBEAM`, referencing `MAT1`) — `parse_stiffness_bdf` builds
  node coordinates, connectivity, per-element orientation (`y_vector`) and per-element 6x6 stiffness.
- **Mass model** (one or more `.nsb` files, `CONM2` entries only — e.g. one structural, one fuel) —
  `parse_mass_files` sums them into one lumped 6x6 mass matrix per node.
- **Aerodynamic mesh** (`.bdf`, `CAERO1` panels + `SPLINE1`/`SET1`) — `read_aero_bdf` reads it, and
  `parse_spline_surfaces` groups CAERO1 panels into continuous surfaces using the model's own `SPLINE1`/`SET1`
  entries, instead of requiring you to hand-list which CAERO eids and structural GRID ids belong together.

Everything derivable from these files is handled automatically: node topology, stiffness, lumped mass, which
CAERO1 panels form one continuous surface (including auto-combining a chordwise split, e.g. a main surface +
control surface, and a spanwise split, e.g. inboard/outboard discretisation), each surface's aerodynamic grid
rediscretized to sit exactly at the structural node positions, the aero-to-structure `dof_mapping`, and (for a
named control surface) the `grid_func` that deflects it. Everything else — flight condition, wake truncation,
boundary constraints, gravity, mirror symmetry — has no NASTRAN-derivable value and is a plain keyword argument
to `build_coupled_aeroelastic`, mirroring the underlying flapjax constructor it's forwarded to; see that
function's own docstring for the full parameter list.

```python
import jax.numpy as jnp
from flapjax.aero import ConstantFlowField
from flapjax_model_gen.nastran import (
    parse_stiffness_bdf, parse_mass_files, read_aero_bdf, parse_spline_surfaces, reorder_by_surfaces,
    SurfaceGridSpec, build_coupled_aeroelastic,
)

beam_model = parse_stiffness_bdf("stiffness.bdf")
aero_model = read_aero_bdf("aero.bdf")
surfaces = parse_spline_surfaces(aero_model)  # one SplineSurface per SPLINE1/SET1 grouping

beam_model = reorder_by_surfaces(beam_model, [s.node_ids for s in surfaces], span_axis=1)
m_lumped, _ = parse_mass_files(
    ["structure.nsb", "fuel.nsb"], node_ids=beam_model.node_ids, node_coords=beam_model.coords
)

wing = SurfaceGridSpec(name="wing", surface=surfaces[0], span_axis=1, m=8, m_star=20)
case = build_coupled_aeroelastic(
    beam_model,
    [wing],
    aero_model,
    m_lumped=m_lumped,
    dt=0.01,
    flowfield=ConstantFlowField(u_inf=jnp.array([30.0, 0.0, 0.0]), rho=1.225, relative_motion=True),
)
```

### Multiple aerodynamic surfaces

A model with more than one surface (wing, tail, fin, ...) just needs more `SurfaceGridSpec`s — no manual
`grid_shapes`/`x0_aero`/`dof_mapping` bookkeeping. Each spec names its surface (an arbitrary label, just used to
cross-reference control surfaces below) and its own `span_axis` — e.g. a wing/horizontal tail spans in y, a
vertical tail in z; a single shared axis can't sort both correctly:

```python
surfaces = parse_spline_surfaces(aero_model)
# parse_spline_surfaces's own order is ascending SET1-id, not semantically meaningful -- identify each
# surface by its caero_ids/node_ids (from your model), not by position.
wing_surface, htail_surface, vtail_surface = surfaces

beam_model = reorder_by_surfaces(
    beam_model,
    [wing_surface.node_ids, htail_surface.node_ids, vtail_surface.node_ids],
    span_axis=[1, 1, 2],
)

specs = [
    SurfaceGridSpec(name="wing", surface=wing_surface, span_axis=1, m=8, m_star=20),
    SurfaceGridSpec(name="htail", surface=htail_surface, span_axis=1, m=6, m_star=20),
    SurfaceGridSpec(name="vtail", surface=vtail_surface, span_axis=2, m=6, m_star=20),
]
case = build_coupled_aeroelastic(beam_model, specs, aero_model, m_lumped=m_lumped, dt=..., flowfield=...)
```

`reorder_by_surfaces` only groups/orders the *structural* side; it's a plain `reorder_nodes(beam_model,
node_order)` under the hood, callable directly if you want to reorder by something other than a spline
grouping. Only `SPLINE1` is supported (raises `NotImplementedError` for `SPLINE2`-`5`, which pyNastran does
parse into an object but this package doesn't handle). `SPLINE6`/`SPLINE7` aren't implemented by pyNastran
itself at all -- they land in `model.reject_cards`, not `model.splines` -- and a `SPLINE7`'s surface can't be
recovered reliably from the raw card either, so any surface using one of these is skipped entirely, with a
`UserWarning`, rather than guessed at or left to raise.

### Control surfaces

Name a control surface with its spanwise extent and chordwise hinge fraction; `build_coupled_aeroelastic` locates
it on the already-built grid and wires it into a `grid_func` that exposes the name as a deflection-angle keyword:

```python
from flapjax_model_gen.nastran import ControlSurfaceRequest, nudged_x_over_c

# a hinge fraction rarely lands on a nice multiple of 1/m -- nudged_x_over_c keeps the rest of the chordwise
# discretisation uniform and only adjusts the one station nearest the hinge onto it exactly.
wing = SurfaceGridSpec(name="wing", surface=wing_surface, span_axis=1, x_over_c=nudged_x_over_c(20, 0.73))

case = build_coupled_aeroelastic(
    beam_model,
    [wing],
    aero_model,
    control_surfaces={
        "aileron": ControlSurfaceRequest(surface_name="wing", span_lo=8.0, span_hi=10.0, hinge_frac=0.73),
    },
    m_lumped=m_lumped,
    dt=...,
    flowfield=...,
)

deflected = case.aero.grid_func(case.aero.zeta_b0, aileron=jnp.deg2rad(5.0))
```

### Half (symmetric) models

For a half-model mirrored about a plane (e.g. the aircraft's x-z plane, normal to the global y axis), reduce
the full-model pieces first: `half_model` drops the negative-side structure (keeping the mirror-plane nodes)
and halves the stiffness of any element lying entirely in the plane; `half_model_lumped_mass` does the nodal
equivalent for lumped mass; `restrict_surface_nodes` pares each surface's node list to match. Each surface's
own CAERO1 panels are then filtered automatically (dropping the negative/mirror-plane side) via
`SurfaceGridSpec.mirror_axis`:

```python
from flapjax_model_gen.nastran import half_model, half_model_lumped_mass, reindex_lumped_mass, restrict_surface_nodes

full_node_ids = beam_model.node_ids
beam_model = half_model(beam_model, mirror_axis=1)
m_lumped, _ = reindex_lumped_mass(m_lumped, full_node_ids, beam_model.node_ids)
m_lumped = half_model_lumped_mass(m_lumped, beam_model.coords, mirror_axis=1)
surfaces = [restrict_surface_nodes(s, beam_model.node_ids) for s in surfaces]

wing = SurfaceGridSpec(name="wing", surface=surfaces[0], span_axis=1, m=8, mirror_axis=1)
case = build_coupled_aeroelastic(
    beam_model,
    [wing],
    aero_model,
    m_lumped=m_lumped,
    dt=...,
    flowfield=...,
    mirror_point=jnp.zeros(3),
    mirror_normal=jnp.array([0.0, 1.0, 0.0]),  # tells flapjax's UVLM about the symmetry plane
)
```

### Auxiliary spline nodes

Some models don't splice a `SPLINE1` directly to the beam's own structural nodes -- instead the `SET1` lists a
separate set of auxiliary leading/trailing-edge `GRID` points (defined in their own file, e.g.
`spline_nodes.bdf`), each rigidly tied to one real structural node via an `RBE2` (any `PLOTEL` entries, used
only for plotting, are ignored). `read_spline_nodes_bdf`/`resolve_spline_surfaces` translate a
`parse_spline_surfaces` result's node ids through those `RBE2`s into real structural node ids before calling
`reorder_by_surfaces` -- an LE/TE pair rigidly tied to the same station collapses to that one node id, and an
id not covered by any `RBE2` (i.e. the `SET1` already lists a real structural node) passes through unchanged,
so it's always safe to call:

```python
from flapjax_model_gen.nastran import read_spline_nodes_bdf, resolve_spline_surfaces

surfaces = resolve_spline_surfaces(surfaces, read_spline_nodes_bdf("spline_nodes.bdf"))
```

### Checking panel/surface groupings visually

`plot_caero1_panels` plots each `CAERO1`'s quadrilateral outline in 3D (matplotlib), labelled by its element
id, and `plot_beam_model` does the equivalent for the structural side -- both useful for sanity-checking
`span_axis`/`mirror_axis` choices and `parse_spline_surfaces` groupings before wiring up the full model:

```python
from flapjax_model_gen.nastran import parse_caero1_panels, plot_beam_model, plot_caero1_panels

ax = plot_caero1_panels(parse_caero1_panels(aero_model), surfaces=surfaces)  # colours each panel by its surface
ax.figure.show()  # or ax.figure.savefig("panels.png")

ax = plot_beam_model(beam_model)
ax.figure.show()
```

### Lower-level building blocks

`build_coupled_aeroelastic`/`SurfaceGridSpec` are a convenience layer over functions you can call directly for
finer control: `build_surface_grid` (one surface's aero grid, given its `SplineSurface` and node coordinates)
is itself `combine_surface_panels` (auto-detect/merge a chordwise split into one full-chord panel, order
root-to-tip) + `build_local_grid` (rediscretize). `control_surface_indices` computes the `m_slice`/`n_slice`/
`hinge_axis` flapjax's own `add_control_surface` needs, if you're assembling a `grid_func` by hand.
`dof_mapping_for_surfaces` resolves a surface's node ids to structural indices, if you're assembling a `UVLM`
by hand rather than through `build_coupled_aeroelastic`.

### 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`), `MAT1`, and `SPLINE1` 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`); mass is taken entirely from `CONM2` mass
  files, not `PBEAM`'s `NSM` or any distributed cross-section mass.
- 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) — only the chordwise discretization (`m`/`x_over_c`) is a free choice.
- A control surface's hinge must land exactly on a chordwise station (`nudged_x_over_c`/`hinged_x_over_c`
  build one that does for an arbitrary hinge fraction), and its spanwise extent must land exactly on existing
  structural nodes — there's no clipping at an arbitrary location on either axis.
- There's no NASTRAN `SPC`/`SPC1` parsing, so boundary constraints (a clamped root, a multibody hinge) are
  built as flapjax constraint objects by hand and passed to `build_coupled_aeroelastic(constraints=...)`.
- 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.

## Matching aerodynamic corrections to a twist distribution

A secondary workflow: given 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), find the NASTRAN-
ready twist distribution that best reproduces the corrected forces without needing the correction itself at
solve time. Most of the underlying machinery (`apply_polar_correction`/`strip_alpha`/`project_forcing_to_beam`,
wired into `UVLM` via `polar_data`/`polar_function`) already lives in `flapjax` itself — see its
`models/cantilever_wing/polar_correction.ipynb` tutorial. This package adds:

- `linear_polar.py` — a `PolarFunction` evaluating `cl = cl0 + cla * alpha`, `cm = cm0 + cma * alpha`, `cd = 0`,
  so `cl0`/`cla`/`cm0`/`cma` data drops straight into `UVLM(polar_data=..., polar_function=...)`.
- `rigid.py` — `rotate_hg`, rigidly rotating a structure's reference SE(3) frames, for representing a change in
  angle of attack as a rotation of the geometry rather than the freestream.
- `twist_grid.py` — `make_twisted_grid`, a per-spanwise-station generalisation of flapjax's
  `make_rectangular_grid(..., twist=...)` — 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.
- `twist_match.py` — `match_twist`, a JAX-native Levenberg-Marquardt solver (autodiff Jacobian,
  `jax.lax.while_loop`) finding the twist distribution minimising the residual between a candidate model's
  forces and a target force field.

```python
import jax.numpy as jnp
from flapjax_model_gen import LinearPolar, linear_polar_function, make_twisted_grid, match_twist, rigid_corrected_forces, rotate_hg

# Target forces from the correction data on the rigid flight shape (uvlm_corrected built with
# polar_data=[LinearPolar(cl0, cla, cm0, cma)], polar_function=[linear_polar_function]), at two angles of attack.
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)))

# Match an *uncorrected* model's per-station twist (the free variable) against those targets.
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(structure.hg0, jnp.deg2rad(1.0)))
    return jnp.stack([f0, f1])

result = match_twist(forces_fn, jnp.stack([target_0, target_1]), 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.

## Development

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

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