Similar to the FDM, the FEM is a numerical method for approximating solutions to partial differential equations (PDEs). However, the FEM offers significantly greater flexibility in handling complex computational domains. To
solve a problem using the FEM, the simulation domain is discretized into a mesh composed of smaller, simpler subdomains known as finite elements. These elements are connected at nodes, which are points shared by two or more
elements, along with boundary lines or surfaces. This discretization transforms the original PDE into a system of algebraic equations by formulating the problem locally for each element and assembling the local systems into a
global system. The following sections will introduce the FEM by outlining the key steps involved in this process. To illustrate how the FEM works in practice, we consider the Poisson equation again:
where \(u\) is the unknown scalar field, \(f\) is a known source term, \(u_D\) defines the Dirichlet boundary condition on \(\Gamma _D\), and \(g\) prescribes the Neumann condition on \(\Gamma _N\). This also serves as a
blueprint for more complex equations, such as the LLG equation. For general background on the finite element method, see [137, 138, 139, 140].
4.2.1 Weak Formulation and Variational Form
The first step in the FEM procedure is to rewrite the original PDE in its weak formulation. This approach relaxes the regularity requirements on the solution and test functions by integrating the governing equation with
respect to a set of smooth test functions. We multiply both sides of the equation by a test function \(v\) and integrate over the domain \(\Omega \):
\(\seteqnumber{0}{4.}{34}\)
\begin{equation}
\int _\Omega (-\nabla ^2 u)\, v \, \mathrm {d}x = \int _\Omega f\, v \, \mathrm {d}x.
\end{equation}
Applying integration by parts, Gauss’s (divergence) theorem:
and using the product rule for differentiation, the second-order derivative term is expressed as a sum of first-order derivatives. Here, the left-hand side is a volume integral over the domain \(\Omega \), while the right-hand side
represents a surface integral over the boundary \(\partial \Omega \).
This results in the following weak formulation:
\(\seteqnumber{0}{4.}{35}\)
\begin{equation}
\int _\Omega \nabla u \cdot \nabla v \, \mathrm {d}x - \int _{\partial \Omega } \frac {\partial u}{\partial n} \, v \, \mathrm {d}s = \int _\Omega f\, v \, \mathrm {d}x,
\end{equation}
where the term \(\frac {\partial u}{\partial n} = \nabla u \cdot \mathbf {n}\) represents the normal derivative of \(u\) on the boundary \(\partial \Omega \).
We consider the boundary \(\partial \Omega \) to be partitioned into two disjoint parts:
where \(\Gamma _D\) is the portion where Dirichlet boundary conditions are prescribed, and \(\Gamma _N\) is the portion with Neumann boundary conditions. Before applying the Dirichlet condition strongly, the boundary integral
over \(\Gamma _D\) remains part of the weak form. It reflects a general formulation suitable for alternative treatments (e.g., Nitsche’s method1). The resulting weak form is:
\(\seteqnumber{0}{4.}{36}\)
\begin{equation}
\int _\Omega \nabla u \cdot \nabla v \, \mathrm {d}x - \int _{\Gamma _D} \frac {\partial u}{\partial n} \, v \, \mathrm {d}s = \int _\Omega f\, v \, \mathrm {d}x + \int _{\Gamma _N} g\, v \,
\mathrm {d}s
\end{equation}
Assuming the Dirichlet condition \(u = u_D\) is imposed strongly (i.e., incorporated directly into the trial function space), the test function \(v\) vanishes on \(\Gamma _D\). The weak form then simplifies to:
\(\seteqnumber{0}{4.}{37}\)
\begin{equation}
\int _\Omega \nabla u \cdot \nabla v \, \mathrm {d}x = \int _\Omega f\, v \, \mathrm {d}x + \int _{\Gamma _N} g\, v \, \mathrm {d}s, \label {eq::weak_form}
\end{equation}
where \(g\) denotes the prescribed Neumann boundary portion, such that \(\frac {\partial u}{\partial n} = g\) on \(\Gamma _N\).
This weak formulation is commonly expressed in an abstract variational form:
where the bilinear form \(a(u, v)\) and the linear form \(L(v)\) are defined as:
\(\seteqnumber{0}{4.}{39}\)
\begin{align}
a(u, v) & = \int _\Omega \nabla u \cdot \nabla v \, \mathrm {d}x, \\ L(v) & = \int _\Omega f\, v \, \mathrm {d}x + \int _{\Gamma _N} g\, v \, \mathrm {d}s.
\end{align}
We assume that \(f \in L^2(\Omega )\), \(g \in L^2(\Gamma _N)\), and \(u_D \in H^{1/2}(\Gamma _D)\) to ensure well-posedness of the weak formulation.
The test function and the solution are assumed to belong to Hilbert spaces. In the Galerkin method, employed for the results presented in this work, both are chosen from the same Hilbert space.
The space \(V\) of test functions is a subspace of the Sobolev space \(H^1(\Omega )\), defined as:
\(\seteqnumber{0}{4.}{41}\)
\begin{equation}
V = \left \{ v \in H^1(\Omega ) \ \middle | \ v = 0 \text { on } \Gamma _D \right \},
\end{equation}
and the trial space is given by:
\(\seteqnumber{0}{4.}{42}\)
\begin{equation}
U = \left \{ u \in H^1(\Omega ) \ \middle | \ u = u_D \text { on } \Gamma _D \right \}.
\end{equation}
More explicitly, the Sobolev space \(H^1(\Omega )\) is defined as:
This formulation is referred to as weak because it only requires equality in an integral sense, whereas the strong form demands pointwise equality and more regularity.
In practical computations, both the solution and test functions are approximated in finite-dimensional subspaces \(U_h \subset U\) and \(V_h \subset V\), which typically consist of globally continuous, piecewise affine (i.e.,
linear on each mesh element) basis functions defined over the computational mesh.
With all components defined, the final variational problem is expressed as follows:
\(\seteqnumber{0}{4.}{45}\)
\begin{equation}
\text {Find } u \in U \text { such that } a(u, v) = L(v) \quad \forall v \in V.
\end{equation}
1 Nitsche’s method enforces Dirichlet boundary conditions weakly by modifying the bilinear form. For details, see [141].
4.2.2 Finite Element Discretization
To solve the weak form numerically, the domain \(\Omega \) is subdivided into a finite number of smaller subdomains forming a triangulation \(\mathcal {T}_h\). These are typically triangles or quadrilaterals in two dimensions,
tetrahedra or hexahedra are common in three dimensions. The set of all mesh nodes is denoted by \(\mathcal {N} = \{ \mathbf {x}_i \}\), where \(\mathbf {x}_i\) are the coordinates of node \(i\).
Based on this triangulation, we define the finite element space \(V_h\), a finite-dimensional subspace of \(V\), consisting of globally continuous, piecewise affine functions:
\(\seteqnumber{0}{4.}{46}\)
\begin{equation}
V_h := \left \{ v \in C^0(\Omega ) \ \middle | \ v|_T \text { is affine } \forall T \in \mathcal {T}_h \right \}.
\end{equation}
Each function \(v_h \in V_h\) is uniquely determined by its nodal values. The set of nodes \(\mathcal {N}\) is partitioned into two disjoint subsets: \(\mathcal {N}_D\), consisting of nodes on the Dirichlet boundary \(\Gamma
_D\), and \(\mathcal {N}_f\), the set of remaining (free) nodes. A nodal basis \(\{\phi _i\}\) is chosen such that:
here, \(U_i\) are the nodal coefficients (degrees of freedom) corresponding to basis functions \(\phi _i\), and they form the components of the global solution vector \(\mathbf {U} = (U_1, U_2, \dots , U_N)^\top \).
Substituting the basis expansion \(u_h\) and \(v_h = \phi _j\) into the weak form (4.38) and requiring it to hold for all basis functions
leads to the component-wise system:
The system is then partitioned into known (Dirichlet) and unknown (free) degrees of freedom (DOF). The free DOFs correspond to the nodes in \(\mathcal {N}_f\), for which the solution values \(\mathbf {u}_f\) are not
prescribed and must be computed by solving the system. In contrast, the Dirichlet DOFs \(\mathbf {u}_D\), associated with nodes in \(\mathcal {N}_D\), are imposed directly according to the boundary conditions. Letting
\(\mathbf {U} = (\mathbf {u}_D, \mathbf {u}_f)^\top \), and partitioning \(A\) and \(\mathbf {b}\) accordingly, written as:
where \(I\) is the identity matrix used to enforce Dirichlet values.
Figure 4.3: Representation of the original solution \(u\) and the finite element approximation \(u_h\) in one (a) and two (b) spatial dimensions. The basis functions for all nodes are shown at the bottom of the 1D plot. The
highlighted basis function \(\varphi _2(x)\) and corresponding solution value \(u_2\) at node \(x_2\) are shown in orange. The 2D plot shows the smooth function \(u(x, y)\), its piecewise approximation \(u_h(x, y)\), and
the basis function \(\varphi _{5,6}(x, y)\) centered at node \((x_5, y_6)\).
The final reduced symmetric system for the unknown values \(\mathbf {u}_f\) is then:
This reduced system involves only the unknown free DOFs and forms the core of the numerical solution process. This system can be solved efficiently because of the sparsity of \(A\), which arises from the local support of the basis
functions. Moreover, the framework naturally accommodates mixed boundary conditions. It can be extended using higher-order basis functions and adaptive mesh refinement to improve accuracy in regions where the solution varies
rapidly.
Figure 4.3 illustrates the representations of the original solution and its finite element approximation in both the 1D and 2D settings. In the 1D case (Figure 4.3 (a)), the nodes are aligned along a single axis, and piecewise linear basis functions approximate the solution. In the 2D case (Figure 4.3 (b)), the domain is discretized into triangles, and the hat function associated with a single node is shown along with a smooth solution surface and projection. The label
\(\varphi _{5,6}(x_5, y_6)\) denotes the basis function and nodal value associated with the corresponding node. This construction naturally extends to higher-order approximations by increasing the polynomial degree of the
basis functions and adding more nodes per element. It also generalizes to 3D problems, where the domain is typically discretized using tetrahedral elements and basis functions are defined over a reference tetrahedron in a similar
manner [142].
