FEM Core

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

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

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

For development, from the repo root:

uv sync
source .venv/bin/activate

Building a model directly

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:

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:

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:

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.

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:

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.

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:

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):

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:

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.

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:

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 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) 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().

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:

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)
  • 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.

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:

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.