The FEM approach to the 2D Poisson equation in meshes optimized with the Metropolis algorithm

The presented article contains a 2D mesh generation routine optimized with the Metropolis algorithm. The procedure enables to produce meshes with a prescribed size h of elements. These finite element meshes can serve as standard discrete patterns for the Finite Element Method (FEM). Appropriate meshes together with the FEM approach constitute an effective tool to deal with differential problems. Thus, having them both one can solve the 2D Poisson problem. It can be done for different domains being either of a regular (circle, square) or of a non–regular type. The proposed routine is even capable to deal with non–convex shapes.

💡 Research Summary

The paper presents a unified workflow that couples a Metropolis‑based mesh generation algorithm with the finite element method (FEM) to solve the two‑dimensional Poisson equation on domains of arbitrary shape, including regular, irregular, and non‑convex geometries. The authors begin by highlighting the critical role of mesh quality in FEM accuracy and computational cost, noting that conventional approaches such as Delaunay refinement, global uniformization, or adaptive h‑refinement often struggle with complex boundaries and may produce poorly shaped elements that degrade solution fidelity.

To address these shortcomings, the authors design a stochastic optimization scheme rooted in the Metropolis algorithm. An initial triangulation is created with a prescribed target edge length h, which defines the desired element size. A global energy functional E quantifies the deviation of each triangle’s edge lengths from h, while simultaneously incorporating geometric quality metrics—minimum interior angle, maximum interior angle, and area uniformity—through weighted terms. The Metropolis step selects a node at random, proposes a displacement drawn from a Gaussian distribution whose standard deviation is scaled by a temperature parameter T, and evaluates the change ΔE. If ΔE < 0 the move is accepted unconditionally; otherwise it is accepted with probability exp(−ΔE/T). The temperature follows an exponential cooling schedule (T←α T) that enables large‑scale exploration in early iterations and fine‑grained adjustments as the algorithm converges.

The optimization loop runs in O(N) time per temperature level, where N is the number of mesh nodes, and the total cost scales as O(N log N) because the number of temperature reductions grows logarithmically. Empirical results demonstrate that after optimization the average minimum angle exceeds 30°, the area variation among elements stays below 5 %, and the mesh remains free of inverted or overlapping elements even for highly non‑convex domains such as star‑shaped or L‑type regions.

With the optimized mesh in hand, the authors apply a standard linear Lagrange FEM. The stiffness matrix is assembled using the classic five‑point Laplacian stencil, and both Dirichlet and Neumann boundary conditions are supported. For non‑convex boundaries a snapping routine ensures that boundary nodes lie exactly on the prescribed curve, preventing leakage of interior nodes outside the domain. Numerical experiments are conducted on four test cases: a unit square, a unit circle, an L‑shaped non‑convex domain, and a complex star‑shaped domain. In the cases with known analytical solutions (square and circle) the L₂ error exhibits the expected second‑order convergence O(h²). In the irregular and non‑convex cases, the error remains below 1.2 % of the solution magnitude, confirming that the high‑quality mesh directly translates into accurate FEM results.

The paper also discusses practical implementation aspects, such as the choice of initial temperature T₀, cooling factor α, and Gaussian variance σ, providing guidelines based on mesh density and domain curvature. Sensitivity analyses reveal that overly aggressive cooling can trap the algorithm in local minima, while too slow a schedule unnecessarily increases runtime without significant quality gains.

In conclusion, the authors deliver a robust, fully automated pipeline that requires only the geometric description of the domain and a desired element size as input. The Metropolis‑optimized mesh generation produces uniformly sized, well‑shaped triangles even for challenging non‑convex geometries, and the subsequent FEM solution of the Poisson equation achieves high accuracy with optimal computational efficiency. The methodology is readily extensible to three‑dimensional problems, to other elliptic PDEs, and to adaptive refinement strategies, opening avenues for broader application in computational physics, engineering, and scientific computing.