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 * (The dimensions of this are those of F_vd.)
42 * By symmetry, this calculation gives the same answer with R and S
43 * exchanged. Looking at the projection in the RMS plane:
49 * R' ,' 2d" = |SS'| = |RR'| = 2d
51 * `-._ ,' By congruent triangles,
52 * ` M with M' = midpoint of RS,
53 * ,' `-._ |MM'| = |RR'|/2 = d
56 * ,' M' _ , - ' d = |MM'|
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.
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.
69 double hypotD(const double p[D3], const double q[D3]) {
74 K pq[_k]= p[k] - q[k];
78 /* owner and block ought not to be used */
80 return gsl_blas_snrm2(&v);
83 double hypotD2(const double p[D3], const double q[D3]) {
85 K d2= ffsqa(p[k] - q[k], d2);
89 double hypotD2plus(const double p[D3], const double q[D3], double d2) {
90 K d2= ffsqa(p[k] - q[k], d2);
97 # define fma_fast(f1,f2,t) ((f1)*(f2)+(t))
99 #define ffsqa(factor,term) fma_fast((factor),(factor),(term))
101 static const l3_epsison= 1e-6;
103 static double energy_function(const double vertices[N][D3]) {
104 int pi,e,qi,ri,si, k;
105 double m[D3], mprime[D3], b, d2, l, sigma_bd2_l;
108 ri= EDGE_END2(pi,(e+1)%V6); if (r<0) continue;
109 si= EDGE_END2(pi,(e+5)%V6); if (s<0) continue;
110 assert(ri == EDGE_END2(qi,(e+2)%V6));
111 assert(si == EDGE_END2(qi,(e+4)%V6));
113 K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
114 K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
115 b= hypotD(vertices[pi], vertices[qi]);
116 d2= hypotD2(m, mprime);
117 l= hypotD(vertices[ri][k] - vertices[si][k]);
118 l3 = l*l*l + l3_epsilon;
120 sigma_bd2_l += b * d2 / l3;
125 static gsl_multimin_fminimizer *minimiser;
130 gsl_multimin_fminimizer_alloc(