chiark / gitweb /
1 /*
2  * We try to find an optimal triangle grid
3  *
4  * Vertices in strip are numbered as follows:
5  *
6  *     ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3   ___ 4   __
7  *         W-1     W-1      0       0       0       0       0
8  *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
9  *       /    \  /    \  /    \  /    \  /    \  /    \  /    \
10  *     L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___  3  ___ 4
11  *     W-2     W-2     W-2      1       1       1       1       1
12  *       \    /  \    /  \    /  \    /  \    /  \    /  \    /
13  *        \  /    \  /    \  /    \  /    \  /    \  /    \  /
14  *     ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3   __
15  *         W-3     W-3     W-3      2       2       2       2
16  *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
17  *       /    \  /    \  /    \  /    \  /    \  /    \  /    \
18  *    .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
19  *
20  *      _ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3
21  *         1       1       1       1      W-2     W-2     W-2     W-2
22  *       /  \    /  \    /  \    /  \    /  \    /  \    /  \    /
23  *      /    \  /    \  /    \  /    \  /    \  /    \  /    \  /
24  *    L-5 ___ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   __
25  *     0       0       0       0       0      W-1     W-1     W-1
26  *
27  * Node x,y for
28  *   0 <= x < L     x = distance along    L = length
29  *   0 <= y < W     y = distance across   W = width
30  *
31  * Vertices are in reading order from diagram above ie x varies fastest.
32  *
33  * Where necessary, edges are counted anticlockwise from to-the-right:
34  *
35  *                          \2   /1
36  *                           \  /
37  *                        ___ 0   __
38  *                        3    1   0
39  *                           /  \
40  *                         4/   5\
41  *
42  * When we iterate over edges, we iterate first over vertices and then
43  * over edges 0 to 2, disregarding edges 3 to 5.
44  */
46 #define L 10
47 #define W 10
48 #define N (L*W)
50 static double energy_function(const double vertices[N]) {
51   /* Energy has the following parts right now:
52    */
53   /*
54    *  Edgewise expected vertex displacement
55    *
56    *
57    *                Q `-_
58    *              / |    `-_
59    *   R' - _ _ _/_ |       `-.
60    *    .       /   M - - - - - S
61    *    .      /    |      _,-'
62    *    .     /     |  _,-'
63    *    .    /    , P '
64    *    .   /  ,-'
65    *    .  /,-'
66    *    . /'
67    *     R
68    *
69    *
70    *
71    *  Find R', the `expected' location of R, by
72    *  reflecting S in M (the midpoint of QP).
73    *
74    *  Let  d = |RR'|
75    *       m = |PQ|
76    *       l = |RS|
77    *
78    *  Giving energy contribution:
79    *
80    *                               2
81    *                            m d
82    *    E             =  F   .  ----
83    *     vd, edge PQ      vd      l
84    *
85    *  By symmetry, this calculation gives the same answer
86    *  with R and S exchanged.  Looking at the projection in
87    *  the RMS plane:
88    *
89    *
90    *                           S'
91    *                         ,'
92    *                       ,'
93    *    R'               ,'             d" = |SS'| = |RR'| = d
94    *      `-._         ,'
95    *          `-._   ,'
96    *              ` M
97    *              ,'  `-._
98    *            ,'        `-._
99    *          ,'              ` S
100    *        ,'
101    *      ,'
102    *     R
103    *
104    *  We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
105    *  because we want this symmetry and because we're happy to punish
106    *  bending more than uneveness in the metric.
107    *
108    *  In practice to avoid division by zero we'll add epsilon to l and
109    *  the huge energy ought then to be sufficient for the model to
110    *  avoid being close to R=S.
111    */
113 static gsl_multimin_fminimizer *minimiser;
117 {
118   gsl_multimin_fminimizer_alloc(