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    *  (The dimensions of this are those of F_vd.)
41    *
42    *  By symmetry, this calculation gives the same answer with R and S
43    *  exchanged.  Looking at the projection in the RMS plane:
44    *
45    *
46    *                           S'
47    *                         ,'
48    *                       ,'
49    *    R'               ,'             2d" = |SS'| = |RR'| = 2d
50    *      `-._         ,'
51    *          `-._   ,'                 By congruent triangles,
52    *              ` M                   with M' = midpoint of RS,
53    *              ,'  `-._              |MM'| = |RR'|/2 = d
54    *            ,'        `-._
55    *          ,'              ` S       So use
56    *        ,'       M' _ , - '            d = |MM'|
57    *      ,'   _ , - '
58    *     R - '
59    *
60    *  We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
61    *  because we want this symmetry and because we're happy to punish
62    *  bending more than uneveness in the metric.
63    *
64    *  In practice to avoid division by zero we'll add epsilon to l^3
65    *  and the huge energy ought then to be sufficient for the model to
66    *  avoid being close to R=S.
67    */
69 double hypotD(const double p[D3], const double q[D3]) {
70   int k;
71   double pq[D3];
72   gsl_vector v;
74   K pq[_k]= p[k] - q[k];
75   v.size= D3;
76   v.stride= 1;
77   v.vector.data= pq;
78   /* owner and block ought not to be used */
80   return gsl_blas_snrm2(&v);
81 }
83 double hypotD2(const double p[D3], const double q[D3]) {
84   double d2= 0;
85   K d2= ffsqa(p[k] - q[k], d2);
86   return d2;
87 }
89 #ifdef FP_FAST_FMA
90 # define ffsqa(factor,term) fma((factor),(factor),(term))
91 #else
92 # define ffsqa(factor,term) ((factor)*(factor)+(term))
93 #endif
95 static const l3_epsison= 1e-6;
97 static double energy_function(const double vertices[N][D3]) {
98   int pi,e,qi,ri,si, k;
99   double m[D3], mprime[D3], b, d2, l, sigma_bd2_l;
101   FOR_EDGE(pi,e,qi) {
102     ri= EDGE_END2(pi,(e+1)%V6); if (r<0) continue;
103     si= EDGE_END2(pi,(e+5)%V6); if (s<0) continue;
104     assert(ri == EDGE_END2(qi,(e+2)%V6));
105     assert(si == EDGE_END2(qi,(e+4)%V6));
107     K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
108     K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
109     b= hypotD(vertices[pi], vertices[qi]);
110     d2= hypotD2(m, mprime);
111     l= hypotD(vertices[ri][k] - vertices[si][k]);
112     l3 = l*l*l + l3_epsilon;
114     sigma_bd2_l += b * d2 / l3;
115   }
119 static gsl_multimin_fminimizer *minimiser;
123 {
124   gsl_multimin_fminimizer_alloc(