|
|
|
|
|
|
The optimization is done with help of the Metropolis algorithm. The system energy is calculated as a sum of discrepancies between an element volume Ve and assumed element volume V0 = h0321/2/12 where h0 denotes a prescribed length of the edge
| E = | ∑ | (Ve - V0)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 Ei. 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
| pi, new = pi - ks | ∑ | Δ rij |
| ij | ||
| Δ rij = | ∑ | (|pi - pj| - h0) (pi - pj) / (|pi - pj|) |
| j | ||
where ks denotes a shifting strength and |pi - pj| is the length of ij edge.
The value of ks 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/h0 for meshes in a cubic domain [0, π] ✗ [0, π] ✗ [0, π] with a) unique elements created for h0 = 0.22 and b) meshes optimized with the Metropolis algorithm with h0 = 0.2.
Within the Metropolis routine the transition probability is calculated by the formula
where kB is a Boltzmann constant (here set as 1), T temperature and ΔEi = Ei, new - Ei 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 Tmax (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/V0 for the same domain as in figure and for a) unique elements; b) elements with V0 = 9 ✗ 10-4 optimized with the Metropolis algorithm.
Note that for the regular tetrahedron a 10% uncertainty in the edge length h0 (see Fig.2b) causes uncertainty in the element volume Δ(V/V0) ± 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.