Metadata-Version: 2.5
Name: fem-core
Version: 0.1.2
Summary: NGSolve/Netgen FEM model-building core — pure Python, no PySide6
Project-URL: Homepage, https://github.com/gcmartins/fem-core
Project-URL: Repository, https://github.com/gcmartins/fem-core
Author-email: Gustavo Martins <gucmartins@gmail.com>
License-Expression: MIT
License-File: LICENSE
Keywords: fea,fem,finite-element,netgen,ngsolve,simulation
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Developers
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Software Development :: Libraries :: Python Modules
Requires-Python: >=3.10
Requires-Dist: h5py>=3.10
Requires-Dist: ngsolve>=6.2.2404
Requires-Dist: numpy>=1.24.0
Requires-Dist: scipy>=1.11.0
Provides-Extra: viz
Requires-Dist: fem-post; extra == 'viz'
Description-Content-Type: text/markdown

# FEM Core

NGSolve/Netgen finite-element model-building core — pure Python, no PySide6 or Qt
dependency. Runs headlessly (no display required).

![Python Version](https://img.shields.io/badge/python-3.10+-blue.svg)
![License](https://img.shields.io/badge/license-MIT-green.svg)

## Features

- 🧱 **Layered FEM Pipeline** — mesh → function space → weak-form assembly, each stage a
  swappable service
- 🧩 **Typed Model Config** — serializable dataclasses (`FemModel`, `FemFormulation`,
  `BoundaryCondition`, `MeshConfig`, `FunctionSpaceConfig`) with `to_dict()`/`from_dict()`
- 🌊 **Physics Formulations** — linear elasticity and (complex-valued) Helmholtz out of the
  box, plus a `CUSTOM` formulation driven by symbolic weak-form expressions
- 🎯 **Domain Facades** — expressive, discoverable builder APIs (e.g. `AcousticModel`) that
  translate physics vocabulary into the generic `FemModel`
- 🔒 **Safe Expression Evaluation** — boundary condition and custom-form expressions are
  parsed with a restricted AST walker, never `eval()`/`exec()`
- 🗂️ **Geometry Import** — load STEP/IGES/BREP files via `GeometryImportService` and index
  their whole containment hierarchy (solids, faces, edges, points), deterministically by
  geometry, so a `BoundaryCondition` can target an imported region directly by index
- 🧵 **Boundary Conditions** — Dirichlet, Neumann, Robin, and periodic, with symbolic
  expressions or vector-valued tractions
- ✅ **Model Validation** — cross-cutting consistency checks across formulation, function
  space, and boundary conditions before a model is registered
- 🧪 **Dependency-Injectable Services** — every service exposes a `Protocol` so tests can
  inject mocks instead of driving real NGSolve/Netgen calls
- 🖥️ **Headless** — no Qt/PySide6 dependency; runs anywhere NGSolve runs
- 🧮 **Solver/Study Workflow** — run a **static** (steady-state), **modal** (eigensolution),
  or **frequency-sweep** (harmonic) analysis against a built model via one `StudyManager`
- 💾 **Result Storage** — save study results with their model and mesh to one HDF5 file and
  load them back as the same GridFunctions (`StudyManager.save_results`/`load_results`)

## Scope

`fem-core` builds the mesh, function space, and `BilinearForm`/`LinearForm` pair, *and* solves
them: a `Study` (static/modal/frequency-sweep) can be run against a built `FemModel` via
`StudyManager`. What's still out of scope is **transient/time-stepping analysis** — that
needs a time-integration layer this package doesn't build toward yet (see Key Design
Decisions below for why, and which study types apply to which formulations).

## Quick Start

### Installation

```bash
pip install fem-core            # from PyPI
pip install "fem-core[viz]"     # plus fem-post (VTK post-processing)
```

For development, from the repo root:

```bash
uv sync
source .venv/bin/activate
```

### Building a model directly

```python
from fem_core.enums import BoundaryConditionType, FemDimension, FormulationType, FunctionSpaceType
from fem_core.managers.fem_model_manager import FemModelManager
from fem_core.models.boundary_condition import BoundaryCondition
from fem_core.models.fem_formulation import FemFormulation
from fem_core.models.fem_model import FemModel
from fem_core.models.formulation_parameters import LinearElasticityParameters
from fem_core.models.function_space_config import FunctionSpaceConfig
from fem_core.models.mesh_config import MeshConfig
from netgen.occ import OCCGeometry, Box, Pnt

geo = OCCGeometry(Box(Pnt(0, 0, 0), Pnt(100, 20, 10)))

model = FemModel(
    name="Bracket",
    dimension=FemDimension.DIM_3D,
    formulation=FemFormulation(
        formulation_type=FormulationType.LINEAR_ELASTICITY,
        parameters=LinearElasticityParameters(E=210e9, nu=0.3, body_force_y=-9.81),
    ),
    function_space=FunctionSpaceConfig(space_type=FunctionSpaceType.VECTOR_H1, order=2),
    mesh_config=MeshConfig(maxh=5.0),
    boundary_conditions=[
        BoundaryCondition(bc_type=BoundaryConditionType.DIRICHLET, region="fixed_face"),
    ],
)

manager = FemModelManager()
managed = manager.build_model(model, geo=geo.shape)  # registers + mesh -> function space -> forms

managed.mesh, managed.fes, managed.assembled_forms  # ngsolve.Mesh, FESpace, AssembledForms
```

### Building a model with a physics facade

Facades expose domain vocabulary instead of requiring callers to construct
`FemFormulation`/`BoundaryCondition`/`FunctionSpaceConfig` by hand:

```python
from fem_core.facades.acoustic_model import AcousticModel

acoustic = (
    AcousticModel(name="Cabin Acoustics", speed_of_sound=343.0, density=1.225, frequency=250.0)
    .add_sound_hard_wall("rigid_walls")
    .add_impedance_boundary("absorber_panel", impedance=415.0 + 50.0j)
)

fem_model = acoustic.to_fem_model()  # -> FemModel, ready for FemModelManager
```

### Running a study

`StudyManager` wraps a `FemModelManager` and runs a `Study` against one of its built models.
Study types are only valid against certain formulations — see
[Key Design Decisions](#key-design-decisions):

```python
from fem_core.enums import StudyType
from fem_core.managers.study_manager import StudyManager
from fem_core.models.study import Study
from fem_core.models.study_parameters import ModalStudyParameters, StaticStudyParameters

manager = FemModelManager()
managed = manager.build_model(model, geo=geo.shape)  # the LINEAR_ELASTICITY model from above

studies = StudyManager(manager)

static_study = studies.create_study(managed.model_id, Study(parameters=StaticStudyParameters()))
studies.run_study(static_study.study_id)
static_study.result.gridfunction, static_study.result.residual_norm

modal_study = studies.create_study(
    managed.model_id, Study(study_type=StudyType.MODAL, parameters=ModalStudyParameters(num_modes=6))
)
studies.run_study(modal_study.study_id)
modal_study.result.eigenvalues  # omega^2, ascending
```

#### Running studies off the GUI thread

`run_study` is `start_run` → `solve` → `finish_run`. A GUI can split it so the
long solve runs on a worker while the manager stays on its own thread:

```python
run = studies.start_run(study_id)       # owning thread: checks, RUNNING, study_started
result = studies.solve(run)             # worker: pure -- no manager state, no notifications
studies.finish_run(study_id, result)    # owning thread: COMPLETED/FAILED + notifications
```

`start_run` returns a `StudyRun` snapshot (copies of the study and model, plus
the built function space and forms). A study removed while solving is ignored by
`finish_run`.

### Undo and redo

`FemModelManager` and `StudyManager` keep an undo history (50 steps):

- **Recorded:** `create_model`, `update_model` (only when the model actually
  changes), `remove_model`, `create_study`, `remove_study` and `clear_all`.
- **Not recorded:** builds, runs and failed calls.
- **Undoing a model change** drops that model's mesh, function space and
  forms, as `update_model` does.
- **A removed study** comes back with its status and result.

```python
from fem_core.services.undo_redo_service import UndoCommand, UndoRedoService

history = UndoRedoService()                      # one stack for both managers
manager = FemModelManager(undo_redo_service=history)
studies = StudyManager(manager, undo_redo_service=history)

managed = manager.create_model(model)
manager.update_model(managed.model_id, edited_model)
manager.undo()                                   # back to `model`
manager.can_undo, manager.can_redo
manager.clear_undo_history()                     # e.g. after loading a project

# Record your own steps on the same stack; apply(state, is_undo) restores
# `before` on undo and `after` on redo.
history.push(UndoCommand("Rename", {"name": "old"}, {"name": "new"}, apply=my_restore))
```

Pass `UndoRedoService(max_steps=0)` to keep no history. Subclasses get
`_notify_undo_redo_changed(can_undo, can_redo)` after every change.

### Solving in a separate process

A direct factorization is one C call that no thread can interrupt. To make a
long solve stoppable, run it in a child process and kill that instead:

```bash
python -m fem_core.study_runner INPUT.h5 STUDY.json OUTPUT.h5
```

`INPUT.h5` holds the model and its mesh as `ResultStorageService.save(path,
model, mesh, fes, {})` writes them; `STUDY.json` is `Study.to_dict()`. The
model is rebuilt on the mesh, the study is run and its result is saved to
`OUTPUT.h5` under the study's name (load it with `ResultStorageService().load`).
Exit status 0 on success, 1 with the reason on stderr.

### Saving and loading results (HDF5)

`StudyManager.save_results()` writes completed studies to one HDF5 file with their model and
mesh, and `load_results()` rebuilds them exactly: same mesh (curved elements included), same
FESpace, same DOF vectors, so the loaded GridFunctions are the ones that were saved. What is
stored is the raw solution, not a sampled picture: for the 3D cantilever example, the model,
mesh and 6 modes take 0.65 MB, against about 16 MB for its `.vtu` files.

`examples/result_storage_example.py` shows the whole workflow. It solves a static, a modal
and a frequency-sweep study on one beam and saves all three to one file. It then works from
the file alone: it lists the studies, reads one mode's raw values, loads everything, plots
the *loaded* results with fem-post, and checks that they match the saved ones exactly.

```python
studies.save_results("results.h5")  # every completed study (all of one model), keyed by Study.name
studies.save_results("results.h5", [sweep.study_id], append=True)  # add to an existing file

stored = studies.load_results("results.h5")  # None on failure, like the other manager calls
stored.fem_model, stored.mesh, stored.fes
modes = stored.studies["Beam Natural Modes"].result  # a ModalResult with real GridFunctions
```

`ResultStorageService` (`src/fem_core/services/result_storage_service.py`) is the layer
underneath. It adds partial reads: `list_studies(path)` returns each study's definition, type
and step count, and `read_vectors(path, key, steps=[2])` returns raw DOF vectors as numpy.
Neither builds a mesh. `load_model(path)` plus `load_study(path, key, fes, steps=...)` rebuild
only the modes or frequencies you ask for.

File layout: `/model` stores the `FemModel` as JSON plus ndof and complexity. `/mesh` stores
netgen's serialized mesh. Each study is in `/studies/<name>`, with its `Study` as JSON,
`vectors[n_steps, ndof]` (float64 or complex128, one compressed chunk per step), and
`eigenvalues`, `frequencies` or `residual_norms` as they apply. The root attributes record the
format version and the fem-core, ngsolve and netgen versions.

On load, the FESpace is rebuilt from the stored model and checked against the stored ndof and
complexity, so a file that no longer matches fails with an error instead of loading wrong values.

The mesh is netgen's own serialization, a pickle stream. It is loaded with an allow-list of the
netgen mesh and geometry classes only, and anything else is refused before it runs. Still,
only load result files from sources you trust.

### Importing external geometry

`GeometryImportService` loads a STEP/IGES/BREP file into a `netgen.occ` shape and indexes its
whole containment hierarchy — solids, then each solid's faces, each face's edges, each edge's
points. A shape imported from a file has no named regions at all, so without this step every
boundary collapses into netgen's single `"default"` region; `index_shape()` numbers every
region 1..N *within its own kind* (the 3rd face is `(FACE, 3)` regardless of how many solids/
edges/points the shape also has), deterministically by geometry, so a `BoundaryCondition` can
target one directly by index — no name to look up or construct by hand:

```python
from fem_core.services.geometry_import_service import GeometryImportService

loaded = GeometryImportService.load_indexed_geometry("part.step")
shape, boundary_index = loaded

print(boundary_index.describe())  # inspect the whole hierarchy to find "which region is which"

model = FemModel(
    ...,
    boundary_conditions=[
        BoundaryCondition(bc_type=BoundaryConditionType.DIRICHLET, region=3),  # face 3
    ],
)

manager = FemModelManager()
managed = manager.build_model(model, geo=shape)  # shape plugs straight into the normal build pipeline
```

`region=3` alone means "face 3" for a 3D model (or "edge 3" for a 2D one) — `region_kind`
only needs setting to target an `EDGE`/`POINT` index instead (Dirichlet only; see
[Key Design Decisions](#key-design-decisions)):

```python
from fem_core.enums import RegionKind

BoundaryCondition(bc_type=BoundaryConditionType.DIRICHLET, region=7, region_kind=RegionKind.EDGE)
```

One condition can cover several regions: pass a tuple (or list) as `region`, and optionally
a `name` for display. The regions resolve to one NGSolve pattern, `"f1|f3|f7"`, which NGSolve
matches as a full regular expression, so exactly those regions are selected:

```python
BoundaryCondition(bc_type=BoundaryConditionType.DIRICHLET, region=(1, 3, 7), name="Fixed support")
BoundaryCondition(
    bc_type=BoundaryConditionType.NEUMANN, region=(2, 4), name="Tip load", vector_value=(0.0, 0.0, -1e3)
)
```

For a vector-valued formulation (linear elasticity), a Neumann load is a traction and must be
given as `vector_value` (one component per dimension); scalar formulations use `value` or
`expression`. Dirichlet conditions are homogeneous (the field is fixed to zero on the region).

Two solids that merely touch (e.g. two boxes side by side in an assembly) share no topology at
all, so every region still gets exactly one owner and one index with no extra step needed —
see [Key Design Decisions](#key-design-decisions).

### Restricting a formulation to a subdomain

`FemFormulation.domain` restricts where a formulation is solved, leaving the rest of the mesh
with no dofs and no equations at all — not just a weakly-suppressed contribution. Left unset
(the default), a formulation applies mesh-wide exactly as before this field existed. This is
what makes it possible to, say, solve a Helmholtz model in just one box of a two-box assembly
and have the other box be genuinely inactive:

```python
from fem_core.enums import RegionKind

model = FemModel(
    dimension=FemDimension.DIM_3D,
    formulation=FemFormulation(
        formulation_type=FormulationType.HELMHOLTZ,
        parameters=HelmholtzParameters(omega=5.0, source=1.0, wave_speed=340.0),
        domain=1,  # solid 1 only -- see GeometryImportService above for where "1" comes from
    ),
    ...,
)
```

`domain` is a `Union[int, str, None]`, resolved the same way `BoundaryCondition.region` is:
an `int` is an index within `domain_kind` (defaulting to the model's own *material* kind —
`SOLID` in 3D, `FACE` in 2D — the only kind a domain restriction is valid for, since a
material region is always codim 0); a `str` is a raw netgen material name, for hand-built
geometry. `FemModelManager.build_function_space` passes the resolved name straight to
`FunctionSpaceService` as `definedon=`, so the `FESpace` itself gets fewer dofs; every
`FormulationService` weak-form term (including the tabular-Robin sweep-form path, and a
`CUSTOM` formulation's own `dx` symbol) is restricted the same way, so it doesn't matter
whether the `FESpace` passed in was pre-restricted or not.

### Visualizing meshes and results with VTK (fem-post)

Post-processing lives in a separate package, **fem-post** (`packages/fem-post/`, import name
`fem_post`). It is kept out of `fem_core` and does not import it. For now it is a uv workspace
member of this repository, and it is self-contained so it can move to its own repository unchanged.
See [`packages/fem-post/README.md`](packages/fem-post/README.md) for the full documentation.

fem-post is an API built around a `PostProcessor`. It wraps one mesh, sampled once, and creates
named **views** on it: the mesh, a single field, a frequency sweep or a set of mode shapes. Each
view can be updated afterwards and exported to `.vtu` (for [ParaView](https://www.paraview.org/))
or `.png`:

- **Update:** new field data with `update()`, display settings with `configure()`, and the
  current step (mode or frequency) with `select_step()`.
- **Export:** `.vtu` with `export_vtu()`, one step to `.png` with `export_png()`, every step with
  `export_pngs()`.

```python
from fem_post import PostProcessor

post = PostProcessor(built.mesh)
post.add_mesh_view("mesh").export_png("vtk_output/my_mesh.png")  # colored by material region

sweep = post.add_frequency_sweep_view("pressure", sweep_study.result)
sweep.export_vtu("vtk_output/my_sweep")  # one field per frequency, for ParaView
sweep.export_pngs("vtk_output")  # one PNG per frequency

modes = post.add_modal_view("modes", modal.result, deformed=True)
modes.configure(warp_scale=50.0).select_step(1).export_png("vtk_output/mode2.png")
```

fem-post samples geometry and fields through `Mesh.MapToAllElements` on a refined reference
lattice, so curved elements and higher-order fields render correctly.

How fields are shown:

- A complex field (any `HELMHOLTZ` result) is stored as `<name>_re`, `<name>_im` and `<name>_abs`,
  and it is colored by `<name>_abs` unless you set `part="real"` or `part="imag"`.
- A vector field is colored by magnitude.
- `deformed=True` also draws a vector field as the deformed shape it describes. The deformation is
  auto-scaled unless you pass `warp_scale`, and the factor used is printed on the image.

The following examples render their mode shapes and harmonic responses this way:

- `linear_elasticity_modal_example.py` (2D);
- `linear_elasticity_modal_3d_example.py` (3D). It is validated against beam theory in
  `validations/linear_elasticity_modal_3d/linear_elasticity_modal_3d_validation.md`, with an
  executable companion notebook `validations/linear_elasticity_modal_3d/linear_elasticity_modal_3d_validation.ipynb`;
- `linear_elasticity_frequency_sweep_example.py`.

fem-post is installed by fem-core's `viz` extra. The example scripts need it:

```bash
uv sync --extra viz
xvfb-run -a uv run python examples/frequency_sweep_example.py
```

Off-screen rendering needs an OpenGL context. On headless Linux, `xvfb-run -a` provides one. Every
example exports its mesh/result `.vtu` files and renders PNGs into
`examples/vtk_output/<example name>/`. That directory is gitignored, and each example gets its own
subdirectory so same-named fields from different examples cannot collide.

Like `linear_elasticity_modal_3d_example.py`, the acoustic `helmholtz_modal_3d_tube_example.py`
(the modes of a rigid, closed cylindrical tube) is validated against the exact Bessel-function solution in
`validations/helmholtz_modal_3d/helmholtz_modal_3d_validation.md`, with an executable companion
notebook `validations/helmholtz_modal_3d/helmholtz_modal_3d_validation.ipynb`.

#### Interactive viewing (needs a real display, run locally)

`view.show()` opens an interactive window with mouse-driven rotate, pan and zoom. For a view with
several steps (modes, frequencies), `n`/`Right` and `p`/`Left` flip between them. The window needs a
**real display attached to the machine it runs on**. It cannot be shown from a remote Claude Code
session. Do not run it under `xvfb-run`, which would give it an invisible display.

## Architecture

FEM Core follows the same layered, signal-free style as `cad-core`/`cad-widgets`, minus the
Qt widget layer:

```
Facade Layer → Manager Layer → Service Layer → Models Layer
```

- **Models** (`src/fem_core/models/`) — Dataclasses describing a FEM problem:
  - `FemModel` — top-level aggregate (dimension, formulation, function space, mesh config,
    boundary conditions)
  - `FemFormulation` / `formulation_parameters.py` — physics type plus its typed parameters
    (`LinearElasticityParameters` — including `density`, a material property used by
    `EigenSolveService` — `HelmholtzParameters` — including `wave_speed`, used by
    `FrequencySweepService` — `CustomParameters`); `domain`/`domain_kind` optionally restrict
    the formulation to one material region of the mesh (see
    [Restricting a formulation to a subdomain](#restricting-a-formulation-to-a-subdomain))
  - `FunctionSpaceConfig` — FESpace type, order, complex-valued flag, and a
    `dirichlet_boundaries` field derived from boundary conditions
  - `MeshConfig` — `maxh`, per-region `maxh`, curvature order, grading, optimization steps
  - `BoundaryCondition` — type (Dirichlet/Neumann/Robin/Periodic), optional display `name`,
    region (an int index within `region_kind`, a raw netgen region name, or a tuple of them),
    value or symbolic expression, Robin coefficient, vector value (traction), or periodic
    partner region
  - `boundary_index.py` — `IndexedRegion` / `BoundaryIndex`, produced by
    `GeometryImportService.index_shape()`: every region of an imported shape's containment
    hierarchy (solid/face/edge/point), each numbered 1..N within its own kind with a
    deterministic name/centroid/measure, inspectable via `BoundaryIndex.describe()`
  - `Study` / `study_parameters.py` — analysis-run type plus its typed parameters
    (`StaticStudyParameters`, `ModalStudyParameters`, `FrequencySweepStudyParameters`)
  - `study_result.py` — `StaticResult`/`ModalResult`/`FrequencySweepResult`, runtime
    containers holding live `ngsolve.GridFunction` objects (not serializable — see Key Design
    Decisions)

  All models except the study results support `to_dict()`/`from_dict()` for serialization.

- **Services** (`src/fem_core/services/`) — Business logic with no manager/UI dependencies,
  each behind a `Protocol` for dependency injection in tests:
  - `GeometryImportService` — loads an external CAD file (STEP/IGES/BREP) into a
    `netgen.occ` shape via `import_shape()`, and walks its whole containment hierarchy
    (solids, faces, edges, points) via `index_shape()`, assigning each region a deterministic,
    geometry-sorted index and name within its own kind (`load_indexed_geometry()` does both
    in one call)
  - `MeshService` — builds NGSolve meshes from 2D/3D netgen geometry (or wraps a pre-built
    mesh), applies curving and per-region `maxh`
  - `FunctionSpaceService` — builds NGSolve FESpace objects (H1, L2, VectorH1, HCurl, HDiv,
    NumberSpace) from `FunctionSpaceConfig`
  - `FormulationService` — assembles unassembled `BilinearForm`/`LinearForm` pairs per
    `FormulationType` (linear elasticity, Helmholtz, custom); folds a
    `-omega**2*density*InnerProduct(u, v)` dynamic-stiffness term into the linear elasticity
    form whenever `LinearElasticityParameters.omega` is nonzero
  - `BoundaryConditionService` — resolves symbolic BC expressions and builds Neumann/Robin
    integrator terms
  - `ExpressionEvaluator` (`expression_evaluator.py`) — restricted AST-based `safe_eval()`
    used by both `FormulationService` and `BoundaryConditionService`
  - `FemModelValidationService` — cross-cutting consistency checks (formulation ↔ function
    space, boundary condition completeness, mesh config sanity), syncs
    `dirichlet_boundaries`/`dirichlet_bboundaries`/`dirichlet_bbboundaries` from boundary
    conditions, and (once a mesh is built) checks every resolved region actually exists in it
    via `validate_regions_against_mesh()`
  - `LinearSolveService` — solves a single assembled system for a `STATIC` study
    (`a.Assemble(); f.Assemble(); u = a.mat.Inverse(...) * f.vec`)
  - `EigenSolveService` — builds a mass matrix from the model's
    `LinearElasticityParameters.density` or `HelmholtzParameters.wave_speed`
    (`1/wave_speed**2`, matching `FrequencySweepService`'s mass-term scaling) and
    solves the generalized eigenproblem `K*phi = omega^2*M*phi` on the free-dof
    submatrices (via `scipy.sparse.linalg.eigsh`, shift-invert) for a `MODAL` study
  - `FrequencySweepService` — converts each `FREQUENCY_SWEEP` study frequency (Hz) into
    `omega` (`2*pi*f`), the angular driving frequency every compatible formulation's
    parameters carry (`HELMHOLTZ` and `LINEAR_ELASTICITY` today), then re-solves the model
    at each point, reusing the built mesh/FESpace and rebuilding only the forms (via
    `FormulationService`) per point
  - `StudyValidationService` — checks a `Study` is compatible with its target model's
    formulation type (see Key Design Decisions) and that its parameters are sane
  - `ResultStorageService` — saves study results with their `FemModel` and mesh to HDF5 and
    loads them back as GridFunctions on a rebuilt, verified FESpace; also reads raw DOF
    vectors of selected steps without building anything

- **Managers** (`src/fem_core/managers/`) — Orchestration layer.
  - `FemModelManager` maintains a registry of `ManagedFemModel` instances and drives the
    build pipeline (`build_mesh` → `build_function_space` → `build_forms`, or `build_all` in
    one call). `build_model(fem_model, geo=...)` is the recommended entry point when geometry
    is already at hand — it registers the model and runs `build_all` in a single call.
  - `StudyManager` composes over a `FemModelManager` (not a subclass of it) and drives the
    solve pipeline: `create_study(model_id, study)` validates and registers a `Study` against
    a built model, `run_study(study_id)` dispatches to the matching solve service and stores
    the result on the returned `ManagedStudy`. One model can have multiple studies
    (`get_studies_for_model`). `save_results(path)`/`load_results(path)` persist completed
    studies through `ResultStorageService`.

  Both managers are pure Python with no Qt dependency — subclasses can override their
  `_notify_*` hook methods to react to state changes (e.g. to emit Qt signals in a UI layer),
  but the base classes themselves have none.

- **Facades** (`src/fem_core/facades/`) — Small stateful builders that give each physics
  domain an expressive setup API instead of forcing callers through raw
  `FemModel`/`FemFormulation`/`BoundaryCondition` construction:
  - `PhysicsModelBuilder` (`base.py`) — base class; subclasses implement `_function_space()`
    and `_formulation()` and inherit `to_fem_model()`
  - `AcousticModel` — frequency-domain Helmholtz acoustics (sound-hard/sound-soft walls,
    impedance boundaries, frequency/speed of sound/density), translated into the generic
    Helmholtz `omega`; `frequency=0.0` builds a model for a `MODAL` (natural-frequency) study
    instead

## Key Design Decisions

- **No Qt dependency.** `fem-core` is pure Python and runs headlessly; any UI integration
  layer subclasses `FemModelManager`/`StudyManager` and overrides their `_notify_*` hooks
  rather than the reverse.
- **Facades are builders, not models.** A facade (e.g. `AcousticModel`) is a small stateful
  object exposing domain vocabulary; `FemModel` (via `to_fem_model()`) remains the single
  serialization boundary.
- **Formulations stay domain-agnostic.** `FormulationType.HELMHOLTZ` only knows about
  `omega`/`source`/`wave_speed` — acoustics-specific inputs (speed of sound, density,
  frequency) are translated into `omega` by the `AcousticModel` facade, so the same
  Helmholtz formulation is reusable by other physics domains.
- **Expressions are never `eval()`'d.** `BoundaryCondition.expression` and
  `FemFormulation.custom_bilinear_form`/`custom_linear_form` are parsed by a restricted AST
  walker (`expression_evaluator.safe_eval`) that only permits numeric literals, `+ - * / **`,
  and calls to a caller-supplied namespace — no attribute access, imports, or arbitrary code
  execution.
- **`dirichlet_boundaries` is derived, not hand-maintained.** It is recomputed from the
  model's `DIRICHLET` boundary conditions by
  `FemModelValidationService.sync_dirichlet_boundaries()` whenever a model is created or
  updated through `FemModelManager`.
- **Every service is DI-friendly.** `MeshServiceProtocol`, `FunctionSpaceServiceProtocol`,
  `FormulationServiceProtocol`, `BoundaryConditionServiceProtocol`, and
  `GeometryImportServiceProtocol` let tests inject mocks into `FemModelManager` without a real
  NGSolve/Netgen call.
- **`BoundaryCondition.region` identifies a mesh region either way: by index or by name.**
  An `int` is an index within `region_kind` (defaulting to the model's own boundary kind —
  `FACE` in 3D, `EDGE` in 2D), resolved to its actual netgen name via `region_name(mesh_dim)`;
  a `str` is a raw netgen region name used directly, unaffected by `region_kind` — the same
  mechanism hand-built geometry (`SplineGeometry.AddRectangle(bcs=[...])`) already used.
  `region_name()` is the single place that resolution happens; every ngsolve-facing service
  (`BoundaryConditionService`, `FormulationService`, `FunctionSpaceService`) calls it rather
  than touching `region` directly.
- **Boundary indices are geometry-sorted and hierarchical, not file-order.** Raw netgen/OCC
  traversal order is only reproducible for one already-loaded shape; regenerating the source
  STEP file with the originating CAD tool can shift it and silently reassign a
  `BoundaryCondition` to the wrong region. `GeometryImportService.index_shape()` instead walks
  the shape's containment hierarchy — solids, then each solid's faces, each face's edges, each
  edge's points — ordering each level by a rounded `(centroid, measure)` key and numbering
  every region 1..N *within its own kind*, so an index keeps meaning the same physical region
  across re-imports of the same geometry.
- **`RegionKind` covers the whole hierarchy (SOLID/FACE/EDGE/POINT), but a `BoundaryCondition`
  can only target a boundary.** `resolved_kind()`/`enums.region_codim()` map each kind to the
  ngsolve selector it reaches at the model's dimension (`Materials`/`Boundaries`/
  `BBoundaries`/`BBBoundaries`) — `SOLID` is always codim 0, a material/volume region, and is
  rejected for any `BoundaryCondition`; `NEUMANN`/`ROBIN` (surface integrals) are further
  restricted to the model's own boundary kind, while `DIRICHLET` is valid on `FACE`, `EDGE`,
  *or* `POINT` (via `dirichlet`/`dirichlet_bbnd`/`dirichlet_bbbnd`, all wired into
  `FunctionSpaceConfig`/`FunctionSpaceService`).
- **Coincident boundaries between touching solids are resolved by construction, not by
  gluing.** Two solids that merely touch (e.g. two adjacent boxes in an assembly) share no
  topology at all — verified against the installed netgen/ngsolve: 0 shared edges, 0 shared
  vertices, and their coincident interface faces compare unequal. `index_shape()` walking the
  hierarchy per-solid therefore gives every region exactly one owner and one index with no
  ambiguity; there is no glue step, and none is needed.
- **A resolved region that doesn't exist in the mesh is a build error, not a silent no-op.**
  `mesh.Boundaries(name)` (and the other selectors) simply select nothing for a name the mesh
  doesn't have — e.g. an out-of-range index — so `FemModelManager.build_mesh` calls
  `FemModelValidationService.validate_regions_against_mesh()` right after meshing, failing the
  build with a message naming what regions of that kind actually do exist.
- **The build pipeline is staged and stateful.** `FemModelManager.build_mesh` →
  `build_function_space` → `build_forms` (or `build_all`) short-circuit on failure and store
  intermediate artifacts on `ManagedFemModel`; updating a model's config resets those
  artifacts to `None` since they no longer match the new configuration. `build_model` is the
  recommended combined entry point for the common case where geometry is already available;
  the staged calls stay available directly for workflows where geometry isn't picked until
  after the model's config is registered (e.g. a UI form filled out before geometry import).
- **A study type is only valid for the formulations it's physically meaningful for:**

  | Study type | Compatible formulations | Why |
  | --- | --- | --- |
  | `STATIC` | `LINEAR_ELASTICITY`, `HELMHOLTZ`, `CUSTOM` | A single linear solve of whatever forms were built; valid for anything. |
  | `MODAL` | `LINEAR_ELASTICITY`, `HELMHOLTZ` | Eigenproblem `K*phi = omega^2*M*phi`, valid when `omega == 0` (no baked-in driving frequency) so the stiffness form is a genuine `K` rather than a dynamic stiffness `K - omega^2*M`. |
  | `FREQUENCY_SWEEP` | `LINEAR_ELASTICITY`, `HELMHOLTZ` | Repeated solves at a varying driving frequency, driven by `omega` on either formulation's parameters — `HELMHOLTZ` derives its spatial wave number `k = omega/wave_speed` from it, `LINEAR_ELASTICITY` uses it directly in a dynamic-stiffness `K - omega^2*M` term, for a harmonic forced-vibration response. |

  `StudyValidationService` enforces this against `StudyType.compatible_formulations`
  (`fem_core.enums`) before a study is registered.
- **Solvers stay generic and read material/medium properties from the model, not from study
  parameters.** A study's parameters describe how to run the analysis (`num_modes`,
  `shift`, a list of `frequencies`, ...); the model describes what's being analyzed,
  material and medium properties included:
  - The modal mass matrix (`density * ∫u·v dx` for `LINEAR_ELASTICITY`, `(1/wave_speed**2) *
    ∫u*v dx` for `HELMHOLTZ`) is assembled in `EigenSolveService` rather than
    `FormulationService`, but `density`/`wave_speed` themselves live on each formulation's own
    parameters — not on `ModalStudyParameters`.
  - `FrequencySweepStudyParameters.frequencies` is just the list of driving frequencies
    (Hz) to solve at. `FrequencySweepService` converts each one into `omega` (`2*pi*f`),
    the angular driving frequency every compatible `FormulationParameters` carries —
    `HelmholtzParameters.omega` and `LinearElasticityParameters.omega` are the same field
    with the same units and conversion, so the sweep drives both formulations identically,
    with no per-formulation branching. Each formulation then derives whatever it actually
    needs from its own generic parameters when building its weak form:
    `FormulationService._build_helmholtz` divides by `HelmholtzParameters.wave_speed` to
    get the spatial wave number `k = omega/wave_speed`; `_build_linear_elasticity`
    multiplies by `LinearElasticityParameters.density` for the `omega**2*density` mass
    term. `wave_speed` and `density` are generic parameters of their own formulation
    (relating `omega` to the formulation's own driven quantity for *any* problem of that
    kind, not just acoustics or a particular structure), unlike domain vocabulary such as
    an acoustic medium's speed of sound, which stays out of the generic formulation layer
    entirely and is only ever resolved inside `AcousticModel`.
- **Study results are runtime containers, like `AssembledForms`.** `StaticResult`,
  `ModalResult`, and `FrequencySweepResult` carry live `ngsolve.GridFunction` objects that
  can't round-trip through JSON; only their scalar/list fields (`residual_norm`,
  `eigenvalues`, `frequencies`) are serializable, and it's on the caller to extract them.
- **A subdomain restriction excludes dofs, not just contributions.** `FemFormulation.domain`
  resolves to an ngsolve `Region` (`mesh.Materials(name)`) that is passed straight to
  `ngs.dx(definedon=...)` at every `FormulationService` weak-form site and, when the `FESpace`
  itself is built through `FemModelManager`, to `definedon=` there too — the unrestricted
  region outside `domain` gets no dofs at all, not a small or zero-valued contribution.
  `ngs.dx(definedon=None)` is a verified no-op passthrough, so every call site uses the same
  `domain` variable unconditionally with no branching for the (default) unrestricted case.

## Testing

Most service tests run against real NGSolve/Netgen (no mocking of the FEM engine itself);
manager/model tests use the injectable protocols to run without it.

```bash
uv sync --group dev --extra viz
invoke check       # lint -> typecheck -> test -> test_post
invoke lint        # ruff linter (fem-core and fem-post)
invoke format      # ruff formatter
invoke typecheck   # ty
invoke test        # pytest (fem-core)
invoke test_post   # pytest (fem-post; run under xvfb-run -a for the rendering tests)
invoke test_cov    # pytest with coverage report

# Or run pytest directly
pytest tests -v
xvfb-run -a pytest packages/fem-post/tests -v
```

## Releasing

`fem-core` and `fem-post` are published to PyPI together, as pure-Python
wheels and sdists, by `.github/workflows/release.yml`. To cut a release,
push `main` to the `release` branch:

```bash
git push origin main:release
```

The workflow then:

1. bumps the patch version of both packages in lockstep and commits it back
   to `release`;
2. builds both packages (`uv build --all-packages`) and runs
   `twine check`;
3. publishes them through PyPI trusted publishing and creates a GitHub
   release `v<version>`.

One-time setup before the first release: on PyPI, add a **pending trusted
publisher** for each project, **fem-core** and **fem-post**, with:

- owner `gcmartins`;
- repository `fem-core`;
- workflow `release.yml`;
- environment `pypi`.

Also create a `pypi` environment in this repository's GitHub settings.
