The optimization is done with help of the Metropolis algorithm. The system energy is calculated as a sum of discrepancies between an element volume V^{e} and assumed element volume V_{0} = h_{0}^{3}2^{1/2}/12 where h_{0} denotes a prescribed length of the edge
E = | ∑ | (V^{e} - V_{0})^{2}. |
e |
Thus the smaller a degeneracy from a designed volume distribution the more optimal state. The Metropolis routine starts from a nodal configuration given by described above procedure. The main point is to reach the optimal global configuration by ascertaining local optimal states. They arise from such a configuration of i-th node and its neighboring nodes which gives smaller energy E^{i}. This partial energy is calculated from the sum equation taken over elements containing the node of interest. To compute new positions for each node (giving new configuration) the following expression is put forward
p_{i, new} = p_{i} - k_{s} | ∑ | Δ r_{ij} |
ij |
Δ r_{ij} = | ∑ | (|p_{i} - p_{j}| - h_{0}) (p_{i} - p_{j}) / (|p_{i} - p_{j}|) |
j |
where k_{s} denotes a shifting strength and |p_{i} - p_{j}| is the length of ij edge.
The value of k_{s} determines the strength of a nodal shift and varies from 0 to 1. It is also worthy considering to choose its value as a random number from uniform distribution U(0,1).
Fig. 2. The figure presents a distribution of a normalized mesh edge length L/h_{0} for meshes in a cubic domain [0, π] ✗ [0, π] ✗ [0, π] with a) unique elements created for h_{0} = 0.22 and b) meshes optimized with the Metropolis algorithm with h_{0} = 0.2.
Within the Metropolis routine the transition probability is calculated by the formula
where k_{B} is a Boltzmann constant (here set as 1), T temperature and ΔE_{i} = E_{i, new} - E_{i} is a difference between energies of these two states. If a value of P is greater than a random number from U(0,1) new state is accepted. Otherwise the old one is preserved.
All above-described local Metropolis steps can lead to different global configurations. Therefore, for each division number this distribution of nodes which gives a lower energy of the total system should be kept. To find this, again the Metropolis rule is employed, but this time changes in the total energy of the whole system are examined. To estimate maximal temperature T_{max} (in expression) a range of changes in potential energy corresponding to current number of elements must be found. Moreover, in each global Metropolis step temperature might decrease according to T = ηT where parameter η < 1.
Fig. 3. The picture presents a distribution of a normalized mesh element volume V/V_{0} for the same domain as in figure and for a) unique elements; b) elements with V_{0} = 9 ✗ 10^{-4} optimized with the Metropolis algorithm.
Note that for the regular tetrahedron a 10% uncertainty in the edge length h_{0} (see Fig.2b) causes uncertainty in the element volume Δ(V/V_{0}) ± 33%. It justifies shapes of histograms in Figs 3, 4d), 5d). They show quality of finally obtained meshes characterized by distributions of elements volume. A weakness in the presented method is fact that it can produce a number of very small elements.