chiark / gitweb /
1 /*
2  * We try to find an optimal triangle grid
3  */
5   /* Energy has the following parts right now:
6    */
7   /*
8    *  Edgewise expected vertex displacement
9    *
10    *
11    *                Q `-_
12    *              / |    `-_
13    *   R' - _ _ _/_ |       `-.
14    *    .       /   M - - - - - S
15    *    .      /    |      _,-'
16    *    .     /     |  _,-'
17    *    .    /    , P '
18    *    .   /  ,-'
19    *    .  /,-'
20    *    . /'
21    *     R
22    *
23    *
24    *
25    *  Find R', the `expected' location of R, by
26    *  reflecting S in M (the midpoint of QP).
27    *
28    *  Let 2d = |RR'|
29    *       b = |PQ|
30    *       l = |RS|
31    *
32    *  Giving energy contribution:
33    *
34    *                               2
35    *                            b d
36    *    E             =  F   .  ----
37    *     vd, edge PQ      vd      3
38    *                             l
39    *
40    *  By symmetry, this calculation gives the same answer
41    *  with R and S exchanged.  Looking at the projection in
42    *  the RMS plane:
43    *
44    *
45    *                           S'
46    *                         ,'
47    *                       ,'
48    *    R'               ,'             2d" = |SS'| = |RR'| = 2d
49    *      `-._         ,'
50    *          `-._   ,'                 By congruent triangles,
51    *              ` M                   with M' = midpoint of RS,
52    *              ,'  `-._              |MM'| = |RR'|/2 = d
53    *            ,'        `-._
54    *          ,'              ` S       So use
55    *        ,'       M' _ , - '            d = |MM'|
56    *      ,'   _ , - '
57    *     R - '
58    *
59    *  We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
60    *  because we want this symmetry and because we're happy to punish
61    *  bending more than uneveness in the metric.
62    *
63    *  In practice to avoid division by zero we'll add epsilon to l^3
64    *  and the huge energy ought then to be sufficient for the model to
65    *  avoid being close to R=S.  */
67 static double hypotD(const double p[], const double q[]) {
68   int k;
69   double pq[D3];
70   gsl_vector v;
72   K pq[_k]= p[k] - q[k];
73   v.size= D3;
74   v.stride= 1;
75   v.vector.data= pq;
76   /* owner and block ought not to be used */
78   return gsl_blas_snrm2(&v);
79 }
81 #ifdef FP_FAST_FMA
82 # define ffsqa(factor,term) fma((factor),(factor),(term))
83 #else
84 # define ffsqu(factor,term) ((factor)*(factor)+(term))
85 #endif
87 static const l3_epsison= 1e-6;
89 static double energy_function(const double vertices[N][D3]) {
90   int pi,e,qi,ri,si, k;
91   double m[D3], mprime[D3], b, d2, l, sigma_bd2_l;
93   FOR_EDGE(pi,e,qi) {
94     ri= EDGE_END2(pi,(e+1)%V6); if (r<0) continue;
95     si= EDGE_END2(pi,(e+5)%V6); if (s<0) continue;
96     assert(ri == EDGE_END2(qi,(e+2)%V6));
97     assert(si == EDGE_END2(qi,(e+4)%V6));
99     K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
100     K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
101     b= hypotD(vertices[pi], vertices[qi]);
102     d2= 0; K d2= ffsqa(m[k] - mprime[k], d2);
103     l= hypotD(vertices[ri][k] - vertices[si][k]);
104     l3 = l*l*l + l3_epsilon;
106     sigma_bd2_l += b * d2 / l3;
107   }
111 static gsl_multimin_fminimizer *minimiser;
115 {
116   gsl_multimin_fminimizer_alloc(