The numerical solution of boundary value problems by the finite element method leads to a system of equations that depends on the triangulation mesh. A change in the geometric configuration of the mesh leads to the changes in the rate of convergence of the process of successive approximations to the numerical solution of the problem. The paper shows how it is possible to speed up the rate of convergence of the iterative solution process without changing the position of the mesh vertices. A class of kernels of boundary value problems is found whose mesh optimization leads to Delaunay triangulation. An example of using the results of the work in the numerical solution of the magnetic field is given.