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, (complex-valued) Helmholtz, Reissner-Mindlin plate bending and Naghdi shells out of the box, plus a CUSTOM formulation driven by symbolic weak-form expressions
• 🌐 Surface Models and Shells — FemDimension.SHELL for surfaces in 3D (faces as the domain, edges as boundaries), with a locking-free linear Naghdi shell (ShellModel)
• 📐 Stress — the stress of elasticity, plate and shell solutions (face stresses, membrane forces, moments, shear forces) as CoefficientFunctions, ready to plot with fem-post (stress_fields)
• 🧮 Multi-field Spaces — PRODUCT function spaces (e.g. a plate's deflection and rotations), with Dirichlet conditions that fix chosen components only
• 🎯 Domain Facades — expressive, discoverable builder APIs (e.g. AcousticModel, PlateModel) 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)
• 📊 Run Statistics Log — opt-in plain-text log of DOFs, matrix nonzeros, assembly/solve times and memory per stage (fem_core.stats_log), on Linux, macOS and Windows

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. Add --log RUN.log to write the run statistics (see below) to a text file.

Run statistics log

Each pipeline stage can log one plain-text line with its size, wall time and process memory to the fem_core.stats logger. The library adds no handler, so the lines are off (and nothing is measured) until you turn them on:

from fem_core.stats_log import enable_stats_log

enable_stats_log("run.log")   # or enable_stats_log() for stderr

You can also configure logging.getLogger("fem_core.stats") yourself. A static run on a 2D beam (unloaded, hence the zero residual) logs:

2026-10-03 01:50:54 | mesh        | 0.003 s | ne=78 nv=62 | rss 204.4 MB (peak 204.9 MB)
2026-10-03 01:50:54 | fespace     | 0.001 s | ndof=402 free=392 complex=False | rss 205.8 MB (peak 205.8 MB)
2026-10-03 01:50:54 | forms       | 0.001 s | rss 209.5 MB (peak 209.5 MB)
2026-10-03 01:50:54 | assemble    | 0.001 s | ndof=402 free=392 complex=False nnz=7884 | rss 211.4 MB (peak 211.4 MB)
2026-10-03 01:50:54 | solve       | 0.002 s | inverse=sparsecholesky | rss 212.9 MB (peak 212.9 MB)
2026-10-03 01:50:54 | residual    | 0.000 s | norm=0
2026-10-03 01:50:54 | study       | 0.003 s | name=Static type=static | rss 212.9 MB (peak 212.9 MB)

The sweep's per-point combine/solve/residual lines are logged at DEBUG, so pass enable_stats_log(level=logging.DEBUG) to see them. A failed stage is marked status=FAILED. Memory is measured with the standard library only. Linux and Windows report current and peak resident memory; macOS reports the peak only.

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, keyed by Study.name
studies.save_results("results.h5", [sweep.study_id], append=True)  # add to an existing file (one model)

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.

Several models on one mesh. When the studies belong to several models that share one mesh (say an elasticity model on one solid and a Helmholtz model on another, each restricted with FemFormulation.domain), save_results writes them to one file, each model under its name (model names must then be unique; append is single-model only). Load them back on one mesh with studies.load_all_results(path) (or ResultStorageService.load_models), which returns the mesh once and a StoredResults per model, all sharing that mesh object. Every reader above takes an optional model_key; list_models(path) lists the keys. Studies of models on different meshes are still refused: save them separately. ResultStorageService.save_models( path, mesh, {key: StoredModel(fem_model, fes, studies)}) is the underlying call.

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. That is format version 1, still what a single-model save writes. Version 2 (several models) keeps /mesh once and moves each model and its studies to /models/<model> and /models/<model>/studies/<name>; a version-1 file reads back as one model under the key "model".

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.

2D models from planar faces

A 2D model (FemDimension.DIM_2D) can be meshed from planar netgen.occ faces, glued like solids are in 3D. GeometryImportService.index_shape then names the faces f1..fN (the 2D materials, what FemFormulation.domain refers to) and the edges e1..eN (the 2D boundaries, what a BoundaryCondition refers to), and FemModelManager.build_mesh(model_id, geo=shape) meshes the shape as a 2D OCCGeometry:

from netgen.occ import Glue, WorkPlane

shape = Glue([WorkPlane().Rectangle(2, 1).Face(), WorkPlane().MoveTo(2, 0).Rectangle(1, 1).Face()])
GeometryImportService.index_shape(shape)  # faces f1, f2; edges e1..e7

The faces may lie in any plane: MeshService.build_mesh_from_geo2d first rotates the shape into the xy plane (into_xy_plane; netgen's 2D mesher drops z), keeping the region names, so the mesh coordinates are in-plane coordinates. planar_frame(shape) gives the plane's origin and normal, or None for a shape that is not planar. A 2D model can therefore use one face of a solid: index the glued 3D geometry once and mesh the chosen faces, which keep their f#/e# names. Faces that belong to no solid (a plate glued next to a body) are indexed after the faces of the solids.

Linear elasticity in 2D is plane strain by default (a long body: the 3D Lamé constants). LinearElasticityParameters(plane_mode=PlaneMode.PLANE_STRESS) models a thin plate instead, with lambda* = 2 mu lambda / (lambda + 2 mu). No thickness is needed: for a uniformly thick body the loads are per unit thickness and the displacement does not depend on it. Plane stress in a 3D model and a nonzero body_force_z in a 2D one are validation errors. tests/test_two_dimensional_real_ngsolve.py checks both modes against beam theory on a cantilever.

Surface (shell) models

FemDimension.SHELL models a surface in 3D space: build it from a shape's faces (an OCC face, or several glued together; meshing a solid instead is reported as an error), and fem-core meshes the surface with surface elements only. Its faces are the model's material regions (a formulation.domain index means a face, f2) and its edges are its boundaries (a boundary condition index means an edge, e3); points are available with region_kind=RegionKind.POINT. Each formulation lists the dimensions it supports (FormulationType.supported_dimensions): HELMHOLTZ (the surface, Laplace-Beltrami, Helmholtz equation) and CUSTOM work on shells, where CUSTOM's dx integrates over the faces and its grad is the surface gradient. SHELL (below) is the shell-structure formulation.

from netgen.occ import Box, Pnt, Sphere
from fem_core.enums import FemDimension

hemisphere = (Sphere(Pnt(0, 0, 0), 1.0).faces[0] * Box(Pnt(-2, -2, 0), Pnt(2, 2, 2))).faces[0]
model = FemModel(dimension=FemDimension.SHELL, ...)  # HELMHOLTZ or CUSTOM
manager.build_model(model, geo=hemisphere)

tests/test_shell_real_ngsolve.py checks the surface Laplacian's spectrum on a hemisphere (l(l+1), exact to 1e-4) and a surface Poisson problem against its closed form. fem-post renders surface meshes directly (see below).

Shell structures (Naghdi)

FormulationType.SHELL is a linear Naghdi (Reissner-Mindlin) shell on a SHELL model: the unknowns are the displacement u and the rotation vector theta, both in global x, y, z, as six H1 fields of one order (at least 2, default 3). Membrane, bending, transverse shear and a small drilling term make up the energy; the membrane strain is interpolated into Regge (HCurlCurl) and the shear strain into H(curl), which keeps thin shells free of membrane and shear locking, and the 3D rotation vector carries moments across folds between faces. Static, modal (translational and rotary inertia) and frequency-sweep studies all apply.

from fem_core.facades.shell_model import ShellModel

shell = ShellModel(E=4.32e8, nu=0.0, thickness=0.25, density=11.19)
shell.set_self_weight((0, 0, -32.174))               # or set_pressure / set_surface_force
shell.rigid_diaphragm("diaphragm", "x")              # uy = uz = 0
shell.symmetry_plane("symmetry_x", "x")              # ux = 0, no rotation about y or z
shell.fix("support", uz=True, rx=True)               # any global axes of u and theta
shell.clamp("root").pin("rim")                       # all six / the three displacements
shell.add_line_force("tip", (0, 0, -5.0))            # force per length on an edge
shell.add_point_force("corner", (1.0, 0, 0))         # a named or indexed geometry point
model = shell.to_fem_model()                         # edges with no condition are free

Each support fixes chosen components of the six-field space (BoundaryCondition.components, indices SHELL_UX ... SHELL_RZ); a named point needs region_kind=RegionKind.POINT, which add_point_force sets. The shell's curvature comes from the curved mesh, so mesh_config.curvature_order must be at least the field order (ShellModel defaults it to order + 1). Results are evaluated on the surface: gridfunction.components[2](mesh(x, y, z, ngsolve.BND)) is uz. ROBIN supports and edge moments are not supported yet.

validations/shell/ validates the formulation against the shell obstacle course with convergence studies (Scordelis-Lo roof, pinched cylinder, pinched hemisphere), exact membrane states (a pressurized cylinder, a sphere's breathing mode), the flat-plate limit, folds (a square tube in torsion against Bredt) and a locking sweep to R/t = 1e5, and explains the benchmarks whose references differ from Naghdi theory. tests/test_shell_structure_real_ngsolve.py runs the core checks, and examples/shell_scordelis_lo_example.py runs the roof.

Conditions on edges and points

A BoundaryCondition usually targets the model's boundary kind (faces in 3D, edges in 2D). region_kind=RegionKind.EDGE or RegionKind.POINT targets a lower-dimensional region of the same index instead, e.g. BoundaryCondition(bc_type=NEUMANN, region=3, region_kind=RegionKind.EDGE, vector_value=(0, 0, -5)) for edge e3:

• Dirichlet fixes the field there (dirichlet_bbnd / dirichlet_bbbnd of the space), e.g. pinning corners against rigid-body motion.
• Neumann integrates the load over the region at its own dimension: per length on a 3D edge (a line load) and at the point itself (a concentrated force). A point force, and to a lesser degree a 3D line load, is singular in continuum theory: the field right there depends on the mesh, while away from it the solution converges (Saint-Venant).
• Robin stays on the boundary kind.

tests/test_lower_dimensional_conditions_real_ngsolve.py checks a line load and a point force against a distributed load of the same total. Point regions are numbered by position like the other kinds (IndexedRegion.centroid is the vertex position).

Plate bending (Reissner-Mindlin)

FormulationType.PLATE models a plate by its mid-surface, a 2D model whose unknowns are the transverse deflection w and the rotation vector theta. They live on a PRODUCT function space with components [H1 (w), VECTOR_H1 (theta)], so a solution's gridfunction.components[0] is the deflection. Both fields use the same polynomial order (at least 2, default 3), and the transverse shear strain grad(w) - theta is interpolated into H(curl) of one order lower before it enters the shear energy (the MITC idea). That keeps thin plates free of shear locking on any mesh, curved elements included. Static, modal (translational and rotary inertia) and frequency-sweep studies all apply.

from fem_core.facades.plate_model import PlateModel

plate = PlateModel(E=210e9, nu=0.3, thickness=0.01, density=7850.0).set_pressure(1000.0)
plate.clamp("left")              # w = 0 and theta = 0
plate.simply_support("right")    # w = 0, rotation free (soft)
plate.hard_simply_support("top") # w = 0, no twist along the edge (hard)
plate.add_elastic_support("bottom", stiffness=1e6)  # transverse springs per length
plate.add_line_force("bottom", 50.0)             # transverse force per length
model = plate.to_fem_model()     # edges with no condition are free

Without the facade, a DIRICHLET condition's components field picks which PRODUCT space components it fixes (components=(0,) is a simple support; None fixes all of them), and plate_function_space(order) builds the matching space.

There are two simple supports, which differ only for thick plates. simply_support() is the soft one (w = 0): the edge may twist, so a thick plate is more flexible, by about 0.8 t/a in centre deflection on a square plate. hard_simply_support() also fixes the rotation along the edge; it is the support of closed-form Mindlin solutions. It works on any edge. On an edge parallel to the x or y axis the rotation along it is a single component, which fem-core fixes exactly (NGSolve's per-direction Dirichlet conditions, with the direction read from the mesh). On an inclined or curved edge it is imposed weakly by Nitsche's method, which converges with the mesh: on a circular plate's rim the remaining tangential rotation is about 1e-5 of the rotation at maxh = R/10.

Named supports

The facades' supports are data in fem_core.facades.supports, so a UI can offer them and recognise them in an existing condition. supports(FormulationType.PLATE | SHELL) lists the Support presets (name, components, tangential_components), and component_names() names the components they fix. Support.condition(region, region_kind) builds the DIRICHLET condition. support_of(bc, formulation_type) returns the preset a condition applies, or None; fixing every component counts as clamped, whether listed or None.

from fem_core.enums import FormulationType, RegionKind
from fem_core.facades.supports import SHELL_PINNED, support_of, supports

[s.name for s in supports(FormulationType.SHELL)]
# ['Clamped', 'Pinned', 'Symmetry plane ⊥ x', ..., 'Rigid diaphragm ⊥ z']
shell.support("p1", SHELL_PINNED, RegionKind.POINT)
support_of(shell.boundary_conditions[-1], FormulationType.SHELL).name  # 'Pinned'

validations/plate_bending/ checks the formulation against Kirchhoff and exact Mindlin theory (static, modal and harmonic, t/a from 1e-4 to 0.2), including a shear-locking test. For curved surfaces in 3D, see Shell structures.

Stress

fem_core.services.derived_fields.stress_fields(model) returns the stress of a solved LINEAR_ELASTICITY, PLATE or SHELL model as CoefficientFunctions. Nothing is solved again, and complex (harmonic) solutions and mode shapes work too. It returns a callable that takes a solution (a GridFunction or a result) and gives {name: field}:

from fem_core.services.derived_fields import stress_fields

stress = stress_fields(managed.fem_model)(static.result)
moment = stress["moment"](managed.mesh(0.5, 0.5))  # a plate's M at its centre, row-major 2 x 2

Sign and axis conventions:

• Faces. For a plate, top is the face at z = +t/2 along +w (the pressure direction). For a shell, it is the face at +t/2 along the surface normal.
• Moments. M = ∫ z σ dz per length, so σ(z) = N/t + 12 M z/t³. A positive moment puts the top face in tension, so a simply supported plate under pressure has a positive centre moment, as in Timoshenko.
• Axes. Shell forces and moments are 3 × 3 tensors tangential to the surface, in global axes.

The plate's and shell's strains are the same locking-free ones the weak form uses:

• the shell's membrane strain is interpolated into Regge elements;
• both models' shear strain is interpolated into H(curl).

These interpolated strains are recovered by a projection, which costs one small sparse solve per field. tests/test_stress_fields_real_ngsolve.py checks:

• the exact stress of linear displacement fields;
• a cantilever's M c/I;
• a hard simply supported plate's centre moment (0.04789 q a², Kirchhoff exact 0.04789) and edge shear force (0.338 q a);
• a tilted clamped shell's centre moment (0.0231 q a²), with its shear force matching the plate's;
• a pressurized cylinder's hoop force q R;
• a square tube's torsional shear flow T/(2A).

The result has the shape fem-post's views take as fields= (see below).

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 an int, a str, a tuple of those, or None, resolved the same way BoundaryCondition.region is. A tuple (a list is converted to one) restricts the formulation to several material regions at once, e.g. domain=(1, 3) for solids 1 and 3; it resolves to the alternation "s1|s3", which ngsolve's mesh.Materials(...) matches as a full-match regular expression, so exactly those regions are selected. An empty tuple is invalid (use None for the whole mesh). Each element is resolved on its own: 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.

• A SHELL model's surface mesh is rendered from its surface elements, colored by face.

• A PLATE or SHELL solution bundles several fields, so pass fields=plate_fields or fields=shell_fields. The view then colors by the deflection (plate) or the displacement (shell), warps by the displacement, and also stores the rotation:

  from fem_post import PostProcessor, shell_fields
  
  post = PostProcessor(managed.mesh, subdivision=3)
  post.add_field_view("roof", static.result, fields=shell_fields, deformed=True).export_png("roof.png")
  

• A stress is plotted with fields=stress_fields(model) (see Stress). A tensor is colored by its von Mises value by default; component= picks a component ("xx", …) or a principal stress, entry= another entry (e.g. a shell's "moment"), and smooth=True averages it across elements within each material region:

  stress = post.add_field_view(
      "stress", static.result, fields=stress_fields(managed.fem_model), smooth=True, deformed=True
  )
  stress.export_png("von_mises_top.png")
  stress.configure(entry="moment", component="xx").export_png("moment_xx.png")
  

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;
• plate_bending_example.py (plate deflection, modes, bending moment and face stress);
• shell_scordelis_lo_example.py (the deformed Scordelis-Lo roof, its modes, the von Mises stress on both faces and the longitudinal membrane force).

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, and PRODUCT spaces of those) from FunctionSpaceConfig
  • FormulationService — assembles unassembled BilinearForm/LinearForm pairs per FormulationType (linear elasticity, Helmholtz, plate, 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 — uses the formulation's registered mass form (build_mass, as PLATE registers), or else builds one 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, LINEAR_ELASTICITY, PLATE and SHELL today), then re-solves the model at each point, reusing the built mesh/FESpace and rebuilding only the forms (via FormulationService) per point
  • stress_fields (derived_fields.py) — a solution's stress, from the formulation's registered build_stress (LINEAR_ELASTICITY, PLATE and SHELL), as CoefficientFunctions in fem-post's field-extractor shape
  • 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
  • PlateModel — Reissner-Mindlin plate bending (thickness, pressure, clamped/simply supported/elastic edges, line and point forces) on a 2D mid-surface model
  • ShellModel — linear Naghdi shells on a SHELL model (thickness, pressure, traction and self-weight, line and point forces, clamped/pinned/symmetry/diaphragm and per-axis supports)

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); a BoundaryCondition targets a region of codim >= 1. 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 (as is any codim-0 kind, e.g. FACE in 2D); ROBIN is restricted to the model's own boundary kind, while DIRICHLET (via dirichlet/dirichlet_bbnd/dirichlet_bbbnd, all wired into FunctionSpaceConfig/FunctionSpaceService) and NEUMANN (integrated at the region's codimension, BoundaryConditionService.region_selector) are valid on FACE, EDGE or POINT (see "Conditions on edges and points").

• 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, PLATE, SHELL, CUSTOM	A single linear solve of whatever forms were built; valid for anything.
  MODAL	LINEAR_ELASTICITY, HELMHOLTZ, PLATE, SHELL	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, PLATE, SHELL	Repeated solves at a varying driving frequency, driven by omega on each formulation's parameters — HELMHOLTZ derives its spatial wave number k = omega/wave_speed from it, LINEAR_ELASTICITY, PLATE and SHELL use 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 format-check  # ruff formatter in check mode (what CI runs)
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.