cufem in Jupyter
Drive the GPU finite-element engine directly from a notebook — the same calls the production worker makes.
cufem is the compiled C++/CUDA engine that everything on the platform runs on: 2-D magnetostatic FEM in float64, assembled and solved on an NVIDIA GeForce RTX 5090 (sm_120). The standard analyses (pmsm_operating_torque, cogging_sweep_2d, …) are thin dispatchers over this engine. In a notebook with engine access you skip the dispatcher and build the problem yourself: a Mesh, one or more material MagnetostaticRegions, Dirichlet boundary nodes, then solve().
Check your device
First confirm the engine imports and sees the GPU. There is exactly one device in the container.
import cufem print(cufem.device_count()) # -> 1 info = cufem.device_info() print(info.name) # -> NVIDIA GeForce RTX 5090 print(info.compute_capability_major, info.compute_capability_minor) # -> 12 0 (sm_120)
If you also have PyTorch in the same container (torch 2.12.0+cu130), torch.cuda.is_available() returns True on the same RTX 5090 — useful for surrogate / ML post-processing alongside the FEM.
The core objects
The whole magnetostatic path is a handful of classes. These are the real compiled types — nothing here is a wrapper or a future plan.
| Object | What it is |
|---|---|
| cufem.Mesh(et, n_nodes, n_elems) | The 2-D mesh. Fill it with set_node(nid, x, y) and set_element(eid, [n0,n1,n2], region_tag). Read-backs (all methods): .n_nodes(), .n_elems(), .bounding_box(), .min_element_area(), .element_type(). |
| cufem.ElementType | Element family: TRIANGLE3 (linear triangles) or QUAD4. |
| cufem.MagnetostaticProblem(mesh) | The solver. set_region_material(tag, region), add_dirichlet_nodes(nodes, value), then solve() → MagnetostaticResult, or solve_torque_vs_angle(spec) → TorqueSweepResult. |
| cufem.MagnetostaticRegion() | Per-region material. Set .mu_r (relative permeability, constant only), and for magnets/coils .M_x, .M_y (magnetization, A/m) and .J_z (out-of-plane current density, A/m²). |
| MagnetostaticResult | Solution fields: .Az (nodal potential), cell fields .Bx, .By, .Bmag, scalar methods .bmax(), .airgap_b_median(); solver health .converged, .cg_iterations, .cg_residual; timings .assembly_time_ms, .solve_time_ms, .total_time_ms. |
| cufem.TorqueSweepSpec() | Defines a torque-vs-angle sweep: .angles_rad, airgap radii .airgap_r1/.airgap_r2, .axial_length, rotor magnet tags/strength .rotor_region_tags/.rotor_M_mag/.rotor_M_angle0_rad, plus optional slot-gate fields. |
| TorqueSweepResult | Sweep output: .torque_Nm (list, one per angle), .net_work_J, .converged (list), .cg_iterations, .total_gpu_ms. |
Geometry / CSG helpers (Circle, Rectangle, Annulus, Union, Difference, build_delaunay_2d, rasterize_sdf_gpu) and I/O (read_gmsh_v2, write_vtk_legacy) round out the module; the magnetostatic objects above are what you need for a motor solve.
A minimal magnetostatic solve
A complete solve has four moves: build the mesh, attach a material to every region tag, pin the outer boundary (Az=0), and call solve(). Below, a small structured grid of TRIANGLE3 elements carries three regions — air (mu_r=1), a soft-iron block (mu_r=1000) and a permanent magnet (mu_r=1.05) whose magnetization comes from its remanence.
/mu_r form is the validated magnet model (it is what production uses) — the helper cufem.br_theta_to_mxy(Br, theta) applies it for you.import cufem, math MU0 = 4e-7 * math.pi MU_R = 1.05 # magnet recoil permeability Br = 1.2 # remanence, T Mmag = Br / (MU0 * MU_R) # validated magnet model: M = Br/(mu0*mu_r) # --- build a small structured TRIANGLE3 grid on [-L, L]^2 -------------------- N, L = 40, 0.05 # 40x40 cells, 50 mm half-box nnx = N + 1 mesh = cufem.Mesh(cufem.ElementType.TRIANGLE3, nnx*nnx, 2*N*N) xs = [-L + 2*L*i/N for i in range(nnx)] bnd = [] for j in range(nnx): for i in range(nnx): nid = j*nnx + i mesh.set_node(nid, xs[i], xs[j]) if i in (0, N) or j in (0, N): bnd.append(nid) # outer boundary nodes e = 0 for j in range(N): for i in range(N): bl = j*nnx + i for nodes in ((bl, bl+1, bl+nnx+1), (bl, bl+nnx+1, bl+nnx)): xc = (xs[nodes[0]%nnx] + xs[nodes[1]%nnx] + xs[nodes[2]%nnx]) / 3 yc = (xs[nodes[0]//nnx] + xs[nodes[1]//nnx] + xs[nodes[2]//nnx]) / 3 r = math.hypot(xc, yc) tag = 4 if r < 0.012 else (3 if r < 0.020 else 1) # 4=magnet 3=iron 1=air mesh.set_element(e, list(nodes), tag); e += 1 # --- materials: one MagnetostaticRegion per tag ----------------------------- prob = cufem.MagnetostaticProblem(mesh) air = cufem.MagnetostaticRegion(); air.mu_r = 1.0 prob.set_region_material(1, air) iron = cufem.MagnetostaticRegion(); iron.mu_r = 1000.0 prob.set_region_material(3, iron) magnet = cufem.MagnetostaticRegion(); magnet.mu_r = MU_R magnet.M_x, magnet.M_y = cufem.br_theta_to_mxy(Br, 0.0) # = (Mmag, 0): +x magnetized prob.set_region_material(4, magnet) # --- boundary + solve ------------------------------------------------------- prob.add_dirichlet_nodes(bnd, 0.0) # Az = 0 on outer box res = prob.solve() # GPU CG solve, float64 print("converged :", res.converged, res.cg_iterations, "it") print("|B| max :", round(res.bmax(), 4), "T") print("airgap Bmed :", round(res.airgap_b_median(), 4), "T") print("solve time :", round(res.total_time_ms, 1), "ms") # res.Bx / res.By are the per-cell flux-density components used for # Arkkio airgap torque (sum r*Br*Bt*area over an airgap annulus).
res.Bx and res.By are the per-element flux-density components — these are exactly the cell fields the operating-torque dispatcher integrates over an airgap annulus to get Arkkio torque. res.Az is the nodal vector potential if you want to plot contours.
Cogging / torque-vs-angle
For a rotor sweep you do not re-build and re-solve in a Python loop — you hand the engine a TorqueSweepSpec and it steps the magnet angle on the GPU, returning one torque per position. Point the rotor magnets at a region tag, give the remanence-derived M magnitude, and the airgap annulus radii for the torque integral.
import cufem, math MU0, MU_R, Br = 4e-7*math.pi, 1.05, 1.2 Mmag = Br / (MU0 * MU_R) # prob is a MagnetostaticProblem built on a motor mesh whose rotor # magnets carry region tag 4 (air=1, stator iron=3, rotor iron=5). spec = cufem.TorqueSweepSpec() spec.angles_rad = [2*math.pi*k/48 for k in range(48)] # 48 positions / rev spec.airgap_r1 = 0.0255 # inner airgap radius, m spec.airgap_r2 = 0.0265 # outer airgap radius, m spec.axial_length = 0.050 # stack length, m spec.rotor_region_tags = [4] spec.rotor_M_mag = [Mmag] # Br/(mu0*mu_r), not Br/mu0 spec.rotor_M_angle0_rad = [0.0] sweep = prob.solve_torque_vs_angle(spec) # one GPU dispatch, all angles T = sweep.torque_Nm # list, one value per angle amp = (max(T) - min(T)) / 2 print("all converged :", all(sweep.converged)) print("cogging amp :", round(amp, 4), "Nm") # small cogging amplitude; see /#validation for the cross-checked number print("GPU time :", round(sweep.total_gpu_ms, 1), "ms")
For the reference 12-slot surface-PM geometry this is the validated 2-D no-load cogging/torque case — independently cross-checked against FEMM and GetDP inside pre-declared thresholds. The full validation fact and numbers live on the homepage (see Validation) so there is a single source of truth. Adding stator current (a .J_z per slot region) moves from no-load cogging to on-load operating torque; that q-axis-current operating-torque path (pmsm_operating_torque, documented on the Motor workflow page) is roadmap / under development — not part of the validated reference, which is the 2-D no-load cogging/torque above.
Import CAD / mesh from gmsh
You are not limited to hand-built grids. The container ships gmsh 4.15.2 with the OpenCASCADE kernel, so you can import a STEP/IGES part, mesh it, write a legacy v2 mesh, and read it straight into a cufem.Mesh with read_gmsh_v2.
import gmsh, cufem # Option A: build the geometry directly in gmsh (self-contained, no external file) gmsh.initialize() gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) # cufem.read_gmsh_v2 needs legacy v2 ASCII gmsh.model.add("rotor") outer = gmsh.model.occ.addDisk(0, 0, 0, 0.05, 0.05) inner = gmsh.model.occ.addDisk(0, 0, 0, 0.02, 0.02) ring, _ = gmsh.model.occ.cut([(2, outer)], [(2, inner)]) # annulus (outer disk - inner disk) gmsh.model.occ.synchronize() gmsh.model.addPhysicalGroup(2, [t for (d, t) in ring], 1) # region tag -> carried into the .msh gmsh.option.setNumber("Mesh.MeshSizeMax", 0.004) gmsh.model.mesh.generate(2) # 2-D mesh (3-D meshing also available) gmsh.write("part.msh") gmsh.finalize() mesh = cufem.read_gmsh_v2("part.msh") # -> cufem.Mesh print(mesh.n_nodes(), mesh.n_elems(), mesh.bounding_box()) # -> 570 1029 (-0.04996, 0.05, -0.04999, 0.04999) (tested output) # Option B: import a real STEP/IGES CAD file instead of building geometry -- # gmsh.model.occ.importShapes("rotor.step") # same flow; keep MshFileVersion = 2.2 # Then assign materials by the physical-group tags, and # prob = cufem.MagnetostaticProblem(mesh); prob.solve() as in the examples above.
To export a solution for ParaView, use cufem.write_vtk_legacy("out.vtk", mesh, nodal_fields=[("Az", res.Az)]).
mu_r, with no nonlinear B–H saturation curve yet (that lands in v0.4). Pick a representative iron permeability (e.g. mu_r=1000) for the operating flux level. (2) VTK precision. write_vtk_legacy truncates coordinates to ~6 significant figures; if you need an exact roundtrip of node positions, keep your authoritative geometry elsewhere and treat the VTK as a visualization export, not a data store.Where to go next
API reference
Every real cufem class and method, plus the analysis catalog, input schemas and the locked result envelope.
Motor workflow
The motor analysis workflow: inputs, 2-mesh convergence, ripple gating and the independent 2-D no-load cogging cross-check (on-load operating torque is roadmap, not validated).
Hard limitations
What the engine does not do today — stated plainly so a number is never mistaken for out-of-scope physics.