2 * We try to find an optimal triangle grid
5 /* Energy has the following parts right now:
8 * Edgewise expected vertex displacement
25 * Find R', the `expected' location of R, by
26 * reflecting S in M (the midpoint of QP).
32 * Giving energy contribution:
40 * By symmetry, this calculation gives the same answer
41 * with R and S exchanged. Looking at the projection in
48 * R' ,' 2d" = |SS'| = |RR'| = 2d
50 * `-._ ,' By congruent triangles,
51 * ` M with M' = midpoint of RS,
52 * ,' `-._ |MM'| = |RR'|/2 = d
55 * ,' M' _ , - ' d = |MM'|
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.
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[]) {
72 K pq[_k]= p[k] - q[k];
76 /* owner and block ought not to be used */
78 return gsl_blas_snrm2(&v);
82 # define ffsqa(factor,term) fma((factor),(factor),(term))
84 # define ffsqu(factor,term) ((factor)*(factor)+(term))
87 static const l3_epsison= 1e-6;
89 static double energy_function(const double vertices[N][D3]) {
91 double m[D3], mprime[D3], b, d2, l, sigma_bd2_l;
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;
111 static gsl_multimin_fminimizer *minimiser;
116 gsl_multimin_fminimizer_alloc(