FEM/BEM Fredkin–Koehler demagnetization¶

Last changes: 23:23 29.08.2026Documentation changelog

The Fredkin–Koehler realization keeps the magnetic volume in FEM and represents the unbounded exterior by a boundary integral. It therefore differs fundamentally from the finite airbox realization: there is no exterior air volume in the body-only mesh, and no Dirichlet or Robin condition is imposed on an artificial outer airbox boundary.

The current implementation is a dense CPU reference path. The page documents the full algorithm, including the two potential solves, surface extraction, solid-angle diagonal, dense operator, linear-solver policy, potential combination, field recovery, energy reduction, and current qualification boundary. It does not promote the generic bem or fmm vocabulary to an executable implementation when the planner has not resolved it to Fredkin–Koehler.

1. Physical problem and domain¶

Let \(\Omega_m\) be a closed tetrahedral magnetic body with boundary \(\Gamma_m\). The exterior is magnetostatic free space. The magnetic source is the volume magnetization \(\mathbf M=M_s\mathbf m\); the exterior contribution is represented through the trace of the potential on \(\Gamma_m\).

The body-only requirement is strict. A mesh with an airbox is a different FEM realization and must use the Poisson-airbox chapter. A non-closed surface cannot define the exterior boundary integral used here and must fail during workspace construction rather than being silently filled with air.

2. Governing equations¶

The scalar potential satisfies

(1)¶\[\mathbf H_{\mathrm d}=-\nabla u, \qquad \Delta u=\nabla\cdot\mathbf M\quad\text{in }\Omega_m.\]

The implementation decomposes the total potential into an interior Poisson part and a harmonic boundary correction:

(2)¶\[u=u_1+u_2, \qquad \int_{\Omega_m}\nabla u_1\cdot\nabla v\,\mathrm dV =\int_{\Omega_m}\mathbf M\cdot\nabla v\,\mathrm dV.\]

The first solve is a Neumann problem and has a constant-potential nullspace. The implementation copies the source right-hand side and pins the first true degree of freedom to zero before the first solve:

(3)¶\[b_{1,0}=0, \qquad u_1\leftarrow u_1+C\ \text{is fixed by the gauge choice }C=0.\]

Let \(\gamma u_1\) be the FEM trace on \(\Gamma_m\). The dense boundary operator maps this trace to the Dirichlet data for the harmonic correction:

(4)¶\[g_2=\mathcal B\,\gamma u_1, \qquad \Delta u_2=0\ \text{in }\Omega_m, \qquad \gamma u_2=g_2\ \text{on }\Gamma_m.\]

The dense diagonal contains the solid-angle contribution. For a boundary node \(x_i\) and adjacent surface triangles \(T\), the diagonal contribution has the implemented form

(5)¶\[\mathcal B_{ii}=\frac{\sum_{T\ni x_i}\omega_T(x_i)}{4\pi}-1.\]

Off-diagonal entries use the linear-triangle boundary weights, including triangle area, oriented unit normal, edge logarithms, and solid angle. Vertex-coincident terms are excluded by the source coincidence guard; this is a singular/self-term treatment, not a general-purpose quadrature approximation that may be removed without changing the operator.

The final field and energy use the combined potential:

(6)¶\[\mathbf H_{\mathrm d}=-\nabla(u_1+u_2), \qquad E_{\mathrm d}=-\frac{\mu_0}{2} \int_{\Omega_m}\mathbf M\cdot\mathbf H_{\mathrm d}\,\mathrm dV.\]

No separate airbox Robin surface energy is added because this realization has no Robin airbox operator. The reported energy is the same magnetic-field functional used by the Poisson path, evaluated with the Fredkin–Koehler recovered field.

3. Symbols and SI units¶

Symbol	Definition	SI unit
\(\mathbf M\)	magnetization field	\(\mathrm{A\,m^{-1}}\)
\(M_s\)	saturation magnetization	\(\mathrm{A\,m^{-1}}\)
\(\mathbf m\)	reduced magnetization	\(1\)
\(\mathbf H_{\mathrm d}\)	demagnetizing field	\(\mathrm{A\,m^{-1}}\)
\(u\)	total scalar potential	\(\mathrm{A}\)
\(u_1\)	interior Poisson potential	\(\mathrm{A}\)
\(u_2\)	harmonic boundary correction	\(\mathrm{A}\)
\(v\)	FEM test function	\(1\)
\(\Omega_m\)	magnetic volume	\(\mathrm{m^3}\)
\(\Gamma_m\)	closed magnetic surface	\(\mathrm{m^2}\)
\(\gamma\)	FEM trace operator on \(\Gamma_m\)	\(1\)
\(\mathcal B\)	dense boundary-integral map	\(1\)
\(g_2\)	boundary values for \(u_2\)	\(\mathrm{A}\)
\(\omega_T\)	oriented solid angle of surface triangle \(T\)	\(1\)
\(T\)	oriented boundary triangle	\(\mathrm{m^2}\)
\(x_i\)	boundary node location	\(\mathrm{m}\)
\(E_{\mathrm d}\)	demagnetization energy	\(\mathrm{J}\)
\(\mu_0\)	vacuum permeability	\(\mathrm{N\,A^{-2}}\)
\(b_1\)	Neumann right-hand-side vector for \(u_1\)	\(\mathrm{A\,m}\)
\(C\)	additive potential gauge constant	\(\mathrm{A}\)
\(i\)	boundary-node index	\(1\)

4. Assumptions and validity limits¶

The implementation assumes a closed, consistently oriented triangular surface extracted from a tetrahedral magnetic body. It uses a dense reference boundary operator; no FMM, hierarchical matrix, fast multipole, or periodic acceleration is implied by the API vocabulary.

The two sparse solves are not interchangeable:

Stage	Problem	Boundary data	Numerical consequence
\(u_1\)	Interior Poisson/Neumann problem	Magnetic RHS plus gauge pin	Constant nullspace must be removed.
\(u_2\)	Interior harmonic/Dirichlet problem	Dense BEM trace \(g_2\)	Boundary values are injected into the FEM RHS before solve.

The dense operator is sensitive to surface orientation, degenerate triangles, node coincidence, solid-angle convention, and boundary-node indexing. A body-only mesh with an open or inconsistent surface is not a lower-accuracy input; it changes or invalidates the mathematical problem.

5. Python API and complete solver parameters¶

The Fredkin–Koehler path uses a closed body-only mesh. The following complete CPU scenario does not declare an airbox: the exterior problem is represented by the boundary operator built from the sphere surface.

# %% Imports and body-only FEM/BEM study
import fullmag as fm

nm = 1.0e-9
study = fm.study("demag_fredkin_koehler_sphere")
study.engine("fem")
study.device("cpu", precision="double")
study.mode("strict")

# %% Closed magnetic body, material, initial state, and volume mesh
body = study.geometry(
    fm.Sphere(radius=30 * nm, name="sphere"),
    name="sphere",
)
body.Ms = 8.0e5
body.Aex = 1.3e-11
body.alpha = 0.1
body.m = fm.init.UniformMagnetization((0.0, 0.0, 1.0))
body.mesh(
    maximum_element_size=12 * nm,
    order=1,
    algorithm_2d=1,
    algorithm_3d=1,
    size_factor=1,
    size_from_curvature=0,
    smoothing_steps=1,
    optimize_iterations=1,
    narrow_regions=0,
    compute_quality=False,
    per_element_quality=False,
)

# %% Fredkin–Koehler interaction and the two-solve FEM policy
study.exchange(enabled=False)
study.demag(model="fredkin_koehler")
study.fem_demag_solver(
    solver="CG",
    preconditioner="AMG",
    rtol=1.0e-10,
    max_iterations=500,
)

# %% Ordered stage and scientific outputs
study.stages.add_relax(
    stage_id="relax",
    algorithm="nonlinear_cg",
    max_steps=2_000,
    tolT=1.0e-6,
).autosave(
    fm.StageAutosave(
        table=fm.TableAutosave(
            every_steps=10,
            quantities=["step", "e_demag", "e_total", "max_torque_T"],
        ),
        fields=[fm.FieldAutosave("H_demag", every_steps=100)],
    )
)

Python parameter	Type	Default	SI unit	Validation	Meaning	Backend support	ProblemIR
`Demag.model`	`optional str`	`None`	\(1\)	`fredkin_koehler` selects the executable body-only FEM/BEM path; `bem` and `fmm` remain separate planner requests.	Selects the FEM/BEM realization family.	FEM CPU reference path; GPU is unsupported.	`energy[].realization`
`Demag.variant`	`optional str`	`None`	\(1\)	Must be omitted for non-airbox models; otherwise constructor validation fails.	Airbox-only boundary selector, not a BEM parameter.	Not applicable to Fredkin–Koehler.	`energy[].realization`
`FEM.order`	`int`	required	\(1\)	Integer greater than or equal to \(1\).	Polynomial order of the interior FEM space.	FEM CPU body-only path subject to mesh support.	`backend_policy.discretization_hints.fem.order`
`FEM.maximum_element_size`	`float`	required	\(\mathrm{m}\)	Finite positive value; `hmax` is an alias.	Target maximum element size for the magnetic body.	FEM CPU.	`backend_policy.discretization_hints.fem.hmax`
`FEM.demag_solver_policy`	`FemLinearSolverPolicy \| None`	`None`	\(1\)	Policy fields are validated individually.	Linear-solver policy for the \(u_1\) and \(u_2\) sparse solves.	FEM CPU; no GPU BEM lane.	`backend_policy.discretization_hints.fem.demag_solver_policy`
`FemLinearSolverPolicy.solver`	`Literal["CG", "GMRES"]`	`"CG"`	\(1\)	Must be `CG` or `GMRES`.	Krylov solver used for each sparse FEM solve.	FEM CPU; both solves use the selected policy.	`demag_solver_policy.solver`
`FemLinearSolverPolicy.preconditioner`	`Literal["AMG", "JACOBI", "NONE"]`	`"AMG"`	\(1\)	Must be `AMG`, `JACOBI`, or `NONE`.	Preconditioner for each sparse solve.	FEM CPU; Hypre path depends on runtime.	`demag_solver_policy.preconditioner`
`FemLinearSolverPolicy.rtol`	`float`	`1e-8`	\(1\)	Finite and strictly positive.	Relative tolerance applied to both \(u_1\) and \(u_2\) solves; it is algebraic, not a physical discretization tolerance.	FEM CPU.	`demag_solver_policy.rtol`
`FemLinearSolverPolicy.atol`	`float \| None`	`None`	\(1\)	When supplied, finite and strictly positive.	Optional absolute stopping tolerance for each sparse solve.	FEM CPU.	`demag_solver_policy.atol`
`FemLinearSolverPolicy.max_iterations`	`int`	`500`	\(1\)	Integer greater than or equal to \(1\).	Maximum iterations for each sparse solve; telemetry reports \(u_1\) and \(u_2\) separately.	FEM CPU.	`demag_solver_policy.max_iterations`
`FemLinearSolverPolicy.print_level`	`int`	`0`	\(1\)	Integer greater than or equal to \(0\).	Solver diagnostic verbosity.	FEM CPU.	`demag_solver_policy.print_level`

6. ProblemIR and normalization¶

The example lowers the interaction to

{
  "kind": "demag",
  "realization": "fredkin_koehler"
}

The FEM hint carries the solver policy separately. Demag.model does not carry rtol or max_iterations; those values belong to FEM.demag_solver_policy because they govern the sparse FEM systems, not the physical demagnetization term. The planner additionally resolves body-only mesh legality, surface closure, available dense operator support, and the requested CPU/GPU lane.

7. Round-trip and failure semantics¶

Requested intent contains fredkin_koehler, the body-only FEM hint, and the complete solver policy. Resolved execution contains the actual surface topology, dense matrix size, solver kind, preconditioner, tolerances, \(u_1\) and \(u_2\) iterations/residuals, total field, energy, and runtime identity.

Invalid model/variant combinations, non-positive tolerances, invalid iteration limits, missing closed surface, degenerate surface triangles, and a failed gauge preparation are validation errors. An unavailable GPU BEM implementation, unsupported generic bem, and unresolved fmm are unsupported combinations. The planner must never silently substitute a Poisson airbox, FDM convolution, or a CPU run for a GPU request.

8. Discrete realization and CPU/GPU separation¶

FEM CPU¶

The CPU sequence is:

assemble the magnetic Poisson RHS;
copy it to the Neumann RHS and pin its first true degree of freedom;
solve for \(u_1\) using the selected sparse solver policy;
extract the boundary trace of \(u_1\);
apply the dense boundary operator to obtain \(g_2\);
prepare the Dirichlet RHS and boundary values for \(u_2\);
solve the harmonic correction problem;
add \(u_1\) and \(u_2\) pointwise;
recover \(\mathbf H_{\mathrm d}\) and reduce energy.

The dense operator build uses oriented triangles, solid-angle diagonal terms, linear-triangle weights, and an explicit node-to-boundary map. The sparse solver may use a cached Hypre operator, preconditioner, and solver. Telemetry records the two solve residuals and iteration counts rather than collapsing them into one number.

FEM GPU¶

No production GPU FEM/BEM lane is claimed. The CUDA Poisson airbox path is a different realization and cannot be described as GPU FEM/BEM. A future GPU BEM implementation would require separate device ownership, dense-operator storage, surface-kernel execution, two-solve residual telemetry, and executed-device validation.

9. Implementation mapping¶

Stage	Stable source owner	Exact responsibility
Source RHS	`assemble_demag_poisson_rhs`	Volume magnetic RHS used by \(u_1\).
Gauge	`prepare_demag_fem_bem_neumann_rhs`	Copies Neumann RHS and pins the gauge DOF.
Dense operator	`DenseDemagBemOperator::build`	Surface matrix, solid-angle diagonal and triangle weights.
Sparse solves	`solve_demag_fem_bem_sparse_system`	Policy-controlled \(u_1\) and \(u_2\) solves.
Boundary values	`prepare_demag_fem_bem_dirichlet_rhs`	Injects dense BEM boundary values into the second solve.
Potential	`combine_demag_fem_bem_total_potential`	Forms \(u=u_1+u_2\).
Orchestration	`context_compute_demag_fem_bem`	Orders all stages and publishes solve telemetry.
Field and energy	`recover_demag_poisson_field` and `demag_poisson_energy_from_field`	Reuses the common field recovery and energy contract.

10. Validation and qualification¶

Test	What it proves	What it does not prove
Closed-body surface extraction	Boundary-node and triangle topology are valid.	Physical convergence of the dense operator.
Solid-angle/self-term contract	Diagonal convention and singular guard are stable.	GPU parity.
Uniform sphere or ellipsoid	Demag factors and energy against an independent open-boundary result.	Arbitrary geometry accuracy.
\(u_1/u_2\) residual telemetry	Both sparse solves satisfy their policies.	Airbox or surface-discretization error.
Energy parity	Combined-potential field and negative half-\(\mu_0\) reduction agree.	Runtime qualification on an untested device.
CPU qualification	Executed CPU result with mesh and policy provenance.	Any GPU implementation.

11. Limitations¶

The current reference path is dense in the number of boundary nodes and is therefore not a large-scale FMM implementation. rtol and max_iterations affect only the sparse interior solves; they do not control dense boundary quadrature error, surface orientation error, or mesh error. A converged sparse solve on an open surface is still an invalid physical result.

The current page deliberately reports GPU FEM/BEM as unsupported. Source files for GPU Poisson demag do not constitute a GPU BEM implementation or executed-device proof.

12. Scientific bibliography¶

D. R. Fredkin and T. R. Koehler, “A boundary-element method for computing the magnetostatic energy of domain walls,” IEEE Transactions on Magnetics, 26, 4154–4156 (1990).
W. F. Brown, Micromagnetics, Wiley, 1963.
FullMag internal reference: docs/physics/0870-fem-bem-demag-open-boundary.md.

13. Source-code index¶

Repository path	Stable symbol	Responsibility	Lane
`packages/fullmag-py/src/fullmag/model/energy.py`	`class Demag`	BEM model validation and IR normalization.	Public authoring
`packages/fullmag-py/src/fullmag/model/discretization.py`	`class FEM`	Body-only mesh and solver-policy hint.	Public authoring
`packages/fullmag-py/src/fullmag/model/discretization.py`	`class FemLinearSolverPolicy`	`rtol`, `atol`, iteration and solver controls.	Public authoring
`backends/fem/cpu/mfem/interactions/demag_fem_bem_rhs.cpp`	`prepare_demag_fem_bem_neumann_rhs`	Gauge-pinned Neumann RHS.	FEM CPU
`backends/fem/cpu/mfem/interactions/demag_fem_bem_operator.cpp`	`DenseDemagBemOperator::build`	Dense surface operator and solid-angle terms.	FEM CPU
`backends/fem/cpu/mfem/interactions/demag_fem_bem_linear_solve.cpp`	`solve_demag_fem_bem_sparse_system`	Policy-controlled sparse solves.	FEM CPU
`backends/fem/cpu/mfem/interactions/demag_fem_bem_boundary_values.cpp`	`prepare_demag_fem_bem_dirichlet_rhs`	Boundary-value injection for \(u_2\).	FEM CPU
`backends/fem/cpu/mfem/interactions/demag_fem_bem_solve.cpp`	`context_compute_demag_fem_bem`	Full Fredkin–Koehler sequence and telemetry.	FEM CPU
`backends/fem/cpu/mfem/interactions/demag_fem_bem_potential.cpp`	`combine_demag_fem_bem_total_potential`	Potential sum.	FEM CPU
`backends/fem/cpu/mfem/interactions/demag_poisson_recovery.cpp`	`recover_demag_poisson_field`	Common field recovery and energy handoff.	FEM CPU
`backends/fem/cpu/mfem/interactions/demag_poisson_energy.cpp`	`demag_poisson_energy_from_field`	Common demagnetization energy reduction.	FEM CPU

Control Room crosswalk¶

Use Model Explorer -> Objects -> <object> -> Physics when PhysicsInteractionPanel exposes the interaction. Status: partial. frontend support is not implemented applies to physical parameters without a matching control. See {doc}/frontend/capability-register; do not infer UI support from backend or Python availability.

Python/API crosswalk¶

The linked Python API page is authoritative for exact functions, arguments, units, and failure semantics. If this page is a foundation or category overview, runnable Python is ot applicable here and must be taken from the terminal API page.

Bibliography and source scope¶

Use the scientific bibliography and source-code index on the linked terminal page. This block adds no new equation or unverified implementation claim.

Source-code index¶

Public Python and lowering sources are linked by the applicable terminal API page. Runtime realization is in the relevant backends/fdm or backends/fem lane; frontend ownership is apps/control-room/src/modules/inspector/panels/PhysicsInteractionPanel.tsx where a live control exists.