Finite Element Method  2D Mesh Generator  Metfem2D
Generator Optimization via the Metropolis method
Internal mesh elements
Presented method is based on the following algorithm^{[1}^{, 2}^{]}:

Define the element size h and consequently the element area A.

Initialize the configuration of triangles and then select the internal nodes P_{int} = { p_{i}: p_{i} ∈ P ∧ p_{i} ∉ Γ } i. e. these nodes does not belong to the domain boundary Γ.

For each node p_{i} in P_{int} find its subdomain Ω^{i} defined as a set of triangles T_{i} to which the node p_{i} belongs.

Perform the Metropolis approach to every internal node p_{i} within its subdomain Ω^{i}. The Metropolis algorithm is adopted in order to adjust an area of each triangle in the node's subdomain to prescribed value A by shifting the position of the node p_{i} (Fig. 5 demonstrates robustness of the Metropolis approach; compare the node distribution in a) and in b)). That adjustment is governed by the following rules:

Find an area of each triangle A_{k} (where k = 1, 2, …, K) in Ω^{i} together with the vectors r_{ji} = p_{i}  p_{j} for each p_{j} ∈ Ω^{i} connected to node p_{i}

Calculate the length of each triangle edge r_{ji} and its deviation δ r_{ji} from the designed element size h i. e. δ r_{ji} = r_{ji}  h

Calculate the new position of the node p^{new}_{i} as
p^{new}_{i} = p_{i} 

∑
_{j}

F_{j} δ r_{ji} versor(r_{ji})

where versor(r_{ji}) (i. e. an unit vector) has the standard meaning as r_{ji}/r_{ji} and F_{j} are weights corresponding to magnitude of jth force applying to node p_{i}. Finally, they were set to the constant value F.

Find an area of each triangle A^{new}_{k} ∈ Ω^{i} after shifting p_{i} → p^{new}_{i}

Apply an energetic measure E to a sub – mesh Ω^{i}. That quantity could be understand in terms of a square deviation of a mesh element area from the prescribed element area A. Therefore, in the presented paper the δ E is defined as a sum of a discrepancy between each triangle area A_{k} and A after moving node p_{i} and prior it, respectively
δE =

∑
_{k}

(
(A^{new}_{k}  A)^{2}  (A_{k}  A)^{2}
)

k = 1, 2, …, K.

If the obtained value of an energetic change is lower than zero the change is accepted. Otherwise, the Metropolis rule is applied i. e. the following condition is checked
e^{δE⁄T} > r
where r is an uniformly distributed random number on the unit interval (0, 1) and T denotes temperature.

The above – presented algorithm is repeated unless an assumed tolerance will be achieved.
In order to reach a better convergence of the presented method several other improvements could be adopted. For instance, the change in the length of the triangle edge could be an additional measure of mesh approximation goodness. That condition will ensure a lack of elongated mesh elements i. e. elements with very high ratio of its edge lengths (to see such skinny elements look at Fig. 5a)).
Figure shows an application of the Metropolis algorithm. Picture a) presents initial positions of new nodes just after generating them whereas picture b) shows their positions after node shifts according to the procedure described in Sec. 4.1 with the following two values of the force strength F_{j}: 0.006 and 0.1 applied to each internal node j = 1, …, 7 and temperature set as T = 0.01. The table presents the total number of Metropolis steps that was required to obtain the final result shown in b).
References