## 2 Numerical Model
### Vectorized Right-Hand Side Assembly Routine  
Transient solvers based on time-stepping algorithms require a more careful implementation with respect to stationary solvers, since many operations must be performed multiple times; typically, once every time step. The postprocessing of the computed variables may also be critical for both the execution time and required memory, especially when dealing with meshes with many nodes or time steps. Aside from the solution of the sparse linear system arising from the finite element discretization of the physical formulation, one of the most time-consuming operations of a FEM solver is the assembly of the global matrix and the right-hand side from the element matrices. Modern interpreted languages such as Python or MATLAB(r) usually provide wrappers for optimized C, C++, or FORTRAN routines for the solution of linear systems; naive assembly techniques are instead based on for loops, which notoriously exhibit low computational efficiencies in interpreted languages [43]. For this reason, assembly operations often constitute the main performance bottleneck in MATLAB(r) finite element codes. In order to make the assembly procedure efficient, the idea introduced in [44] is exploited in FLARE. The terms of all the local element matrices (i.e., the local contributions of each element to the global solution) are computed with a vectorized procedure instead of a for-based procedure. Since the local element matrices are \(3\times 3\) matrices, the contributions must be stored into a \(9\times n_{el}\) array, where \(n_{el}\) is the total number of elements of the given mesh. Two auxiliary arrays with the same size are also used to provide a mapping between the local and global views. The arrays store the row and column indices of the global matrix corresponding to each element of the local values. In this way, the entries of the local contribution array are assembled into the global matrix with a single call to the sparse MATLAB(r) Matlab function. A comprehensive performance comparison of vectorized assembly procedures in different languages can be found in the work of [PERSON] and colleagues [45]. By preassembling the two matrices \([K]\) and \([S]\) in (9), the square sparse matrix \([M_{1}]\) in (8) can be obtained at each time step by summing \([K]\) and \([S]\) after calculating the element-wise products by \(1/2\) and \(1/\Delta t\), respectively. The same technique can be adopted for \([M_{2}]\) multiplying the solution array at the previous time instant. Unfortunately, the array \([T]^{(k+1/2)}\) must be reassembled at each time step with a time-varying \(J_{0,z}\) source term. In order to avoid this operation, we define an auxiliary matrix \([K_{t}]\) with size \(n_{p}\times n_{el}\), where \(n_{p}\) is the number of points (vertices) of the mesh. The aim of \([K_{t}]\) is to allow the array \([T]^{(k+1/2)}\) to be computed as a matrix-vector product after computing the source term at the considered time instant. The array \([Q]\) of size \(n_{el}\times 1\) contains the quantity that each element contributes to each of the three nodes representing its vertices. \([Q]\) is defined by the (weighted) element-wise product between the column-arrays \(\left[I_{0,z}^{(k+1/2)}\right]\) and \([Aa]\) :  
\[[Q]=\frac{1}{3}\mu_{0}\Big{[}I_{0,z}^{(k+1/2)}\Big{]}\odot[Aa], \tag{14}\]  
where \(\left[I_{0,z}^{(k+1/2)}\right]\) and \([Aa]\) store the element current densities evaluated at the time instant \((k+1/2)\) and the element areas, respectively. Let us consider the simple two-element and four-node mesh in Figure 1. Given the represented element and node numbering, the connectivity matrix can be written as  
\[[C]= \begin{bmatrix}1&2&4\\ 2&3&4\end{bmatrix} . \tag{15}\]  
The entries of the auxiliary matrix \([K_{t}]\) are either zeros or ones. The nonzero elements of the \(i\)th column of \([K_{t}]\) correspond to the \(i\)th row of the connectivity matrix \([C]\). In this way, for the considered example, one has\[\left[T^{(k+1/2)}\right]=[K_{t}][Q]=