<p>In the design optimization of electrical machines, highly accurate numerical models are indispensable. However, due to their time-consuming nature, the computational efficiency of the optimization procedure is deteriorated. To get rid of the issue, developing both time-efficient and accurate finite element solvers is pivotal, by which a high number of function evaluations can be executed without the traditional concerns. This paper presents a comparative study of first- and second-order triangular finite element solvers applied to two-dimensional nonlinear magneto-static problems commonly encountered in the analysis of electrical machines. Accordingly, identifying the minimum mesh density required for first-order elements to produce acceptable results enables the efficient use of first-order meshes in large-scale optimization tasks, thereby balancing accuracy with computational efficiency. The finite element formulation is based on the magnetic vector potential, and a detailed implementation of the Newton–Raphson nonlinear solver is presented and thoroughly discussed. Three benchmark cases are considered, including a current-carrying conductor, an L-shaped ferromagnetic domain with nonlinear <i>B</i>(<i>H</i>) curve behavior, and a switched reluctance machine (SRM) featuring pronounced local magnetic saturations in its rotor and stator iron cores. Accuracy, convergence, and computational cost are assessed for both element orders on topologically equivalent meshes. The results demonstrate that while second-order elements are essential for accurately capturing local saturation effects, first-order elements can still deliver sufficiently accurate results in a computationally efficient manner, making them suitable for optimization tasks.</p>