Assembly¶

From weak form to a linear system¶

Each model in μfem states its physics as a weak (variational) form on the computational mesh. Assembly is the step that turns that continuous form into a finite system of algebraic equations for the degrees of freedom \(\mathbf{x}\):

\[ A\,\mathbf{x} = \mathbf{b} . \]

Symbol	Meaning
\(\mathbf{x}\)	Solution vector, the coefficients of the discretized field (the true degrees of freedom)
\(A\)	System (stiffness) matrix, assembled from the bilinear part of the weak form; its entries carry the material properties
\(\mathbf{b}\)	Right-hand side (load vector), assembled from the linear part: sources and boundary data

Element matrices¶

The unknown field is expanded in the basis functions of a finite-element space. On each mesh element \(e\) the weak form reduces to a small dense element matrix \(A_e\) and element vector \(\mathbf{b}_e\). A stiffness entry is an integral over the element of products of the basis functions (or their gradients), weighted by the local material property \(\kappa\):

\[ (A_e)_{ij} = \int_{\Omega_e} \kappa\,\nabla\phi_i \cdot \nabla\phi_j \; \mathrm{d}\Omega . \]

Symbol	Meaning
\(\Omega_e\)	The element domain
\(\phi_i\)	Basis (shape) function of local degree of freedom \(i\)
\(\kappa\)	Material property from the weak form (conductivity, reluctivity, ...)

The reference element¶

To integrate over an arbitrary physical element \(\Omega_e\), it is pulled back to a single fixed reference element \(\hat\Omega\) (for example the unit triangle), on which the basis functions and the quadrature points are defined once. The two elements are related by a mapping that runs in both directions: the parametrization \(\mathbf{x}(\boldsymbol{\xi})\) sends the reference element to the physical one (built from the element's geometry nodes, whose geometry order sets how curved it can be), and its inverse \(\boldsymbol{\xi}(\mathbf{x})\) maps the physical element back to the reference.

Top row: a curved physical element mapped to the reference triangle by the inverse map xi(x) with Jacobian inverse. Bottom row: the reference triangle mapped to the curved physical element by the parametrization x(xi) with Jacobian J equal to the derivative of x with respect to xi

The derivative of the parametrization is the Jacobian matrix

\[ J = \frac{\partial \mathbf{x}}{\partial \boldsymbol{\xi}} , \]

which encodes how the reference coordinates stretch and rotate into physical coordinates. Two quantities derived from it appear in every element integral:

its determinant \(|J|\), the local volume-scaling factor that converts the reference measure to the physical one, \(\mathrm{d}\Omega = |J|\, \mathrm{d}\hat\Omega\);
its inverse transpose \(J^{-\mathsf{T}}\), which maps reference-space gradients to physical ones, \(\nabla_{\mathbf{x}}\phi = J^{-\mathsf{T}}\, \nabla_{\boldsymbol{\xi}}\hat\phi\).

Because everything is pulled back to the reference element, the basis functions and quadrature rule are tabulated once and reused for every element; only \(J\) and \(|J|\), evaluated from each element's geometry, change. For a straight-sided (affine) element \(J\) is constant; for a curved element it varies from point to point.

Numerical integration¶

These integrals are almost never evaluated in closed form. Instead a quadrature rule samples the integrand at a fixed set of points \(\boldsymbol{\xi}_q\) inside the element and forms a weighted sum:

\[ (A_e)_{ij} \approx \sum_{q} w_q \; \kappa(\boldsymbol{\xi}_q)\; \nabla\phi_i(\boldsymbol{\xi}_q) \cdot \nabla\phi_j(\boldsymbol{\xi}_q)\; |J(\boldsymbol{\xi}_q)| , \]

where \(w_q\) are the quadrature weights and \(|J|\) is the Jacobian determinant introduced above, evaluated at each quadrature point. μfem uses Gauss quadrature, whose rule of order \(m\) integrates polynomials up to degree \(m\) exactly. A higher order uses more quadrature points, as shown here on a tetrahedron:

A tetrahedron shown four times with the Gauss quadrature points for integration orders one to four: one point, then four, then five, then eleven points, increasingly filling the element

The order is chosen so the integrand is captured accurately, and it grows with:

the element order \(p\): higher-order basis functions raise the polynomial degree of the integrand;
curved elements, where the Jacobian \(|J|\) varies across the element;
field-dependent properties, where \(\kappa\) is not constant over the element.

Too low an order under-integrates the system and can destabilise the solution; too high an order simply costs more. For an order-\(p\) element with a constant property, a rule exact to degree \(2p\) is a typical default.

Scattering into the global system¶

Assembly then scatters the element contributions into the global system:

\[ A \;=\; \sum_{e} P_e^{\mathsf{T}}\, A_e\, P_e , \qquad \mathbf{b} \;=\; \sum_{e} P_e^{\mathsf{T}}\, \mathbf{b}_e , \]

where \(P_e\) is the local-to-global map that sends each element degree of freedom to its global index. The entries of \(A_e\) are then added into the rows and columns of \(A\) at those global indices; a degree of freedom shared by several elements accumulates a contribution from each. The figure below shows one element, whose three local nodes carry global indices \(\{2, 5, 7\}\), scattering its \(3 \times 3\) matrix into those rows and columns of the global matrix.

A 3x3 element matrix scattered into the rows and columns of a global sparse matrix at the element's global degree-of-freedom indices

Because each element couples only the few degrees of freedom it contains, the assembled matrix \(A\) is large and sparse: most entries are zero, and only the non-zeros are stored.

Before this system is handed to a solver, the boundary and coupling constraints still have to be imposed on \(A\) and \(\mathbf{b}\).