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. Historical development and present – day concepts of finite element analysis are widely described in references . 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 . 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 . Then the corresponding nodal loads 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 , the advancing front method  and tree methods . 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 . 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 , 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).