 |
BELFEM 0.9.0
Berkeley Lab Finite Element Framework
|
|
BELFEM 0.9.0
Berkeley Lab Finite Element Framework
|
Date: 2026-08-14 Purpose: The weak-form derivations behind BELFEM's electromagnetic formulations: the fundamental lemma of variational calculus, the two divergence-theorem corollaries, the least squares projection, a thermal warm-up example, and the b-conform and h-conform weak forms of the quasi-magnetostatic Maxwell equations. Module: src/fem/maxwell Origin: Distilled from Christian Messe's quasi-magnetostatic theory notes (LaTeX, written in the Lagrange-multiplier era of the code). Transcription errata in the notes were fixed during extraction. Each equation block names its TeX source for traceability. Implementation status: BELFEM solves the h-phi formulation only (en_Maxwell_Formulations.hpp). The b-conform derivations (a and a-v) are included for completeness of the theory; they are not implemented as solving formulations. The least squares projection is live in the L2* postprocessing kernels (matrices/mt_maxwell_l2_h.cpp, mt_maxwell_l2_b.cpp, mt_maxwell_l2_phi.cpp).
(source: weakform/fundamental_lemma.tex)
When calculus of variations was discovered, people investigated energy levels of potential functions. Let phi(x) be a potential field that describes a physics problem, and let the functional Pi be a measure for the compulsion the field experiences within a domain Omega. From observation we know that such systems follow Le Chatelier's principle of assuming a state where this compulsion is minimized over the full domain:
This is called the strong form, and G is the governing equation of the problem in an integral sense. Applying the chain rule gives the field that fulfills it:
Some problems, like linear elasticity, can be formulated directly this way. For many other problems, however, the functional Pi is not known or does not exist at all. Galerkin's great discovery is that the statement generalizes to a fundamental lemma:
where g(phi) = 0 is a partial differential equation describing the conservation law of the problem. In this context, delta phi loses its physical meaning as a variation of phi and becomes an abstract test function. But that does not matter: as long as delta phi is smooth and delta Pi = 0 is fulfilled, the conservation law is automatically fulfilled within the domain. The formulation of (flcv) that can be discretized in an appropriate way is called the weak form of the governing equation.
Regardless of the physical problem, the divergence theorem or one of its corollaries is always required to bring the weak form into a suitable shape. Two corollaries are of significant importance. Let c be a scalar field, f and g vector fields in Omega, and n the outward normal on the boundary dOmega. Combining Gauss's theorem with the chain rule gives, for the expression div(c g):
while a combination of Stokes' theorem (the Stokes-Cartan theorem) with the chain rule yields, for div(f x g):
(source: weakform/leastsquares.tex; live in `matrices/mt_maxwell_l2_.cpp`)*
Since the edge degrees of freedom h_tilde have an abstract nature, they cannot be directly visualized. Moreover, the vast majority of visualization tools such as ParaView require a node or element based dataset. The goal is therefore to project the edge based field onto the nodes while minimizing the projection error. One very efficient way is a least-squares mapping, which is one of the easiest finite element problems to formulate. We apply the fundamental lemma to the trivial statement that the magnetic field h equals the magnetic field h:
The test function delta_h introduces a set of virtual degrees of freedom. They are not physical; their sole purpose is to span a linear system that can be solved. For the nodal field we write h ≈ N h_hat; following the idea of Boris Galerkin, the same interpolation is used for the virtual degrees of freedom, delta_h^T ≈ delta_h_hat^T N^T. Since the edge dofs are known, the field can also be expressed as h ≈ E h_tilde. Inserting both yields the system A h_hat = b that projects the edge dofs onto the nodes:
Derivative field properties such as the electric current are computed the same way with the curl operator C in place of E. The least squares mapping is not only used to project fields and their derivatives onto different bases and even meshes; in the Lagrange-multiplier era of these notes it also underlay the interface conditions for mixed formulations. In today's code the interface coupling is hanging-edge condensation, and the L2 kernels serve postprocessing.
(source: weakform/heat.tex; the production thermal solver lives in src/fem/thermal)
The heat conduction problem is one of the simplest problems in physics that can be solved with finite element theory, which makes it the ideal warm-up for constructing a weak form. The governing equation reads:
with density rho, specific heat capacity c, temperature T, heat flux q_dot and a volumetric heat load q_dot_v (imposed, for example, by electromagnetic induction). Having identified T as the degree of freedom, the fundamental lemma gives:
The Gauss corollary (gauss_theorem) moves the divergence onto the test function:
and applying Fourier's law q_dot = -k grad T yields the weak form:
The expression -n^T q_dot is the scalar heat load into the structure. With the interpolation rules
(the notes carried a spurious transpose on the gradient rule; B is (dim x n), so grad T ≈ B T_hat), the discretized problem becomes:
with mass matrix M, stiffness matrix K and load vector f. Iterative solution methods even allow c and k to depend on the temperature. A special thought should be invested into the temperature dependence of the density. Thermal expansion does change rho with temperature; but if the problem is solved on the undeformed mesh Omega_0, mass conservation requires that the reference density rho_0 = rho(T_0) is used at all times:
(source: weakform/maxwell.tex)
The Maxwell equations are Gauss's electric law, Gauss's magnetic law, Faraday's law of induction, and the Ampère-Maxwell equation:
Here d is the displacement field, b the magnetic flux density, h the magnetic field, j the electric current density and e the electric field. The model is closed with Ohm's law and the constitutive equations, which in the absence of polarization and magnetization simplify to:
with electric conductivity sigma (resistivity rho = sigma^-1), permittivity epsilon and reluctivity nu. In general these are spatially dependent tensors; many metals and the vacuum behave isotropically, in which case they collapse to a scalar times the identity.
In the magneto-quasistatic simplification used throughout BELFEM, the change of the displacement field d_,t and the volumetric charges rho_v are neglected.
For thermal conduction it is obvious that the temperature is the field to solve for and that the heat flux is the flux term, so there is only one reasonable way to construct the weak form. This is not the case for the Maxwell equations. Two families of formulations have proven useful: the b-conform formulations and the h-conform formulations.
(source: weakform/bconform.tex; not implemented in BELFEM, derived for completeness)
A vector potential a is postulated whose curl yields the magnetic flux density:
Its existence is justified by the fact that Gauss's magnetic law is always fulfilled: div b = div curl a = 0. The test function delta_a is applied to the Ampère-Maxwell equation (displacement currents neglected):
The Stokes corollary (stokes_theorem) is applied to the left hand side (the notes wrote the volume integral over dOmega; it is over Omega):
Inserting (transh) and (vektorpot) into the first term on the right gives the stiffness expression, and reassembling yields the weak form of the a-formulation:
If the a-formulation is coupled with another formulation such as the h-formulation, the boundary integral is used to formulate the connector element. The form above assumes that all relevant current densities j are known, which is useful, for example, when a current is imposed at a coil.
(source: weakform/bconform_voltage.tex; this file was an orphan in the notes, not compiled into the main document, and is likewise not implemented)
If currents in conducting regions are to be computed, the a-formulation must be extended with a degree of freedom for the electric voltage v. The continuity equation for charge conservation in the quasistatic case reads:
(the notes wrote div j = +rho_v,t; charge conservation carries the minus sign, which is irrelevant here since the right hand side is zero anyway). Applying the fundamental lemma and the Gauss corollary:
Combining Ohm's law with Faraday's law expresses the current in terms of both potentials:
which, inserted into (jweak) with a volume-imposed current density, gives:
and inserted into (weaka):
(source: weakform/hconform.tex; this is the implemented core, see matrices/mt_maxwell_h.cpp)
The fundamental lemma is applied to Faraday's law of induction; the magnetic field h is identified as the degree of freedom, so the test function is named delta_h:
For the first expression, the inverse of (transh) is used, with the product rule picking up a term for a time-varying permeability:
The Stokes corollary is applied to the second expression:
The boundary integral in (hstokes) is of special interest. If a coupling with the b-conform formulation is desired, the electric field can be expressed through the voltage (if it exists) and the vector potential; expressing e through the current density instead allows current coupling:
With Ohm's law e = rho curl h in the volume and v = 0, the weak form of the h-conform formulation reads:
A closer look reveals that the h-formulation is well suited for discretizing a superconducting domain, while the expression mu_,t can become very unhandy to compute if mu is not constant. In those cases the b-conform formulation avoids that term, which is why ferromagnetic domains are the classical b-conform candidates.
In BELFEM, (weakh) is discretized with edge elements in the conductors (edge_h dofs, mt_maxwell_h.cpp) and with the scalar potential h = -grad phi in the non-conducting regions (mt_maxwell_phi.cpp), following the magnetodynamic coupling of Arsenault et al. 2023 (read with the 2026 erratum). A sign-convention note: the notes' operator rule h ≈ B phi_hat is unsigned, and the phi kernel deliberately drops the Arsenault minus on the mass term because integral mu |h|^2 is even in the sign (see the comment in mt_maxwell_phi.cpp); the minus matters wherever the sign survives, for example in an interface term. The coupling at the conductor boundary is performed by hanging-edge condensation, not by an interface weak form; see the maxwell usage guide and cl_Maxwell_TMatrix.
For the element-level operators (N, E, B, C) that discretize these weak forms, see ../../interpolation/doc/nedelec_derivation.md.
1.18.0