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.

The variety of problems in physics or engineering is formulated by appropriate differential equations with some boundary conditions imposed on the desired unknown function or the set of functions. There exists a large literature which demonstrates numerical accuracy of the finite element method to deal with such issues^{[1]}.
Historical development and present – day concepts of finite element analysis are widely described in references ^{[1]}. In Sec. 2 of the paper and in its Appendixes A – D, the mathematical concept of the Finite Element Method is presented.
In presented article the well – known Laplace and Poisson equations will be examined by means of the finite element method applied to *an appropriate mesh*. The class of physical situations in which we meet these equations is really broad. Let us recall such problems like heat conduction, seepage through porous media, irrotational flow of ideal fluids, distribution of electrical or magnetic potential, torsion of prismatic shafts, lubrication of pad bearings and others ^{[2]}. Therefore, in physics and engineering arises a need of some computational methods that allow to solve accurately such a large variety of physical situations.
The considered method completes the above-mentioned task. Particularly, it refers to a standard discrete pattern allowing to find an approximate solution to continuum problem. At the beginning, the continuum domain is discretized by dividing it into a finite number of elements which properties must be determined from an analysis of the physical problem (e. g. as a result of experiments). These studies on particular problem allow to construct so – called *the stiffness matrix* for each element that, for instance, in elasticity comprising material properties like stress-strain relationships ^{[3]}. Then the corresponding *nodal loads*^{[4]} associated with elements must be found.
The construction of accurate elements constitutes the subject of a mesh generation recipe proposed by the author within the presented article. In many realistic situations, mesh generation is a time – consuming and error – prone process because of various levels of geometrical complexity. Over the years, there were developed both semi – automatic and fully automatic mesh generators obtained, respectively, by using the mapping methods or, on the contrary, algorithms based on the Delaunay triangulation method ^{[5]}, the advancing front method ^{[6]} and tree methods ^{[7]}. It is worth mentioning that the first attempt to create fully automatic mesh generator capable to produce valid finite element meshes over arbitrary domains has been made by Zienkiewicz and Phillips ^{[8]}.
The advancing front method (AFM) starts from an *initial* node distribution formed on a basis of the domain boundary, and proceeds through a sequential creation of elements within the domain until its whole region is completely covered by them. The presented mesh algorithm takes advantage from the AFM method as it is demonstrated in Sec. 3. After a node generation along the domain boundary (Sec. 3.1), in next steps interior of the domain is discretized by adding *internal* nodes that are generated at the same time together with corresponding elements which is similar to Peraire *et al.* methodology ^{[9]}, however, positions of these new nodes are chosen differently according to the manner described in Sec. 3.2. Further steps improve the quality of the mesh by applying the Delaunay criterion to triangular elements (Appendix E) and by a node shifting based on the Metropolis rule (Sec. 4).

^{^}O. C. Zienkiewicz, R. L. Taylor and J. Z. Zhu,*The Finite Element Method: Its Basis and Fundamentals, Sixth edition.*, Elsevier 2005^{^}O. C. Zienkiewicz and Y. K. Cheung,*Finite elements in the solution of field problems*, The Engineer, pp. 507-510 1965; O. C. Zienkiewicz, P. Mayer and Y. K. Cheung,*Solution of anisotropic seepage problems by finite elements*, J. Eng. Mech., ASCE,**92**, pp. 111-120, 1966; O. C. Zienkiewicz, P. L. Arlett, and A. K. Bahrani,*Solution of three – dimensional field problems by the finite element method*, The Engineer, 1967; L. R. Herrmann,*Elastic torsion analysis of irregular shapes*, J. Eng. Mech., ASCE,**91**, pp. 11-19, 1965; A. M. Winslow,*Numerical solution of the quasi-linear Poisson equation in a non-uniform triangle 'mesh'*, J. Comp. Phys.,**1**, pp. 149-172, 1966; M. M. Reddi,*Finite element solution of the incompressible lubrication problem*, Trans. Am. Soc. Mech. Eng.,**91**:524 1969^{^}R. W. Clough,*The finite element method in structural mechanics*, In O. C. Zienkiewicz and G. S. Holister, editors, Stress Analysis, Chapter 7. John Wiley & Sons, Chichester, 1965^{^}Nodes are mainly situated on the boundaries of elements, however, can also be present in their interior.^{^}A. Bowyer,*Computing Dirichlet tessellations*, Comp. J.,**24**(2), pp. 162-166, 1981; D. F. Watson,*Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes*, Comput. J.,**24**(2), pp. 167-172 1981; J. C. Cavendish, D. A. Field and W. H. Frey,*An approach to automatic three-dimensional finite element mesh generation*, Int. J. Numer. Meth. Eng.,**21**, pp. 329-347 1985;N. P. Weatherill,*A method for generating irregular computation grids in multiply connected planar domains*, Int. J. Numer. Meth. Eng.,**8**, pp. 181-197 1988; W. J. Schroeder, M. S. Shephard,*Geometry-based fully automatic mesh generation and the Delaunay triangulation*, Int. J. Numer. Meth. Eng.,**26**, pp. 2503-2515 2005; T. J. Baker,*Automatic mesh generation for complex three-dimensional regions using a constrained Delaunay triangulation*, Eng. Comp.,**5**, pp. 161-175 1989; P. L. Georgea, F. Hechta and E. Saltela,*Automatic mesh generator with specific boundary*, Comp. Meth. Appl. Mech. Eng.,**92**, pp. 269-288 1991^{^}S. H. Lo,*A new mesh generation scheme for arbitrary planar domains*, Int. J. Numer. Meth. Eng.,**21**, pp. 1403-1426 1985; J. Peraire, J. Peiro, L. Formaggia, K. Morgan, O. C. Zienkiewicz,*Finite element Euler computations in three dimensions*,**26**, pp. 2135-2159 2005; R. Löhner, P. Parikh,*Three – dimensional grid generation by the advancing front method*, Int. J. Num. Meth. Fluids**8**, pp. 1135-1149 1988^{^}M. A. Yerry, M. S. Shephard,*Automatic three-dimensional mesh generation by the modified-octree technique*, Int. J. Numer. Meth. Eng.,**20**, pp. 1965-1990 1984; P. L. Baehmann, S. L. Wittchen, M. S. Shephard, K. R. Grice, and M. A. Yerry,*Robust, geometrically-based, automatic two-dimensional mesh generation*, Int. J. Numer. Meth. Eng.,**24**, pp. 1043-1078 1987; M. S. Shephard and M. K. Georges,*Automatic three-dimensional mesh generation by the finite octree technique*, Int. J. Numer. Meth. Eng.,**32**, pp. 709-749 1991^{^}O. C. Zienkiewicz and D. V. Phillips,*An automatic mesh generation scheme for plane and curved surfaces by isoparametric coordinates*, Int. J. Numer. Meth. Eng.,**3**, pp. 519-528 1971^{^}J. Peraire, M. Vahdati, K. Morgan, and O. C. Zienkiewicz,*Adaptative remeshing for compressible flow computations*, J. Comp. Phys.**72**, pp. 449-466, 1987