anneal.c

   1 /*
   2  * We try to find an optimal triangle grid
   3  */
   4
   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    */
  68
  69 double hypotD(const double p[D3], const double q[D3]) {
  70   int k;
  71   double pq[D3];
  72   gsl_vector v;
  73
  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 */
  79
  80   return gsl_blas_snrm2(&v);
  81 }
  82
  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 }
  88
  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
  94
  95 static const l3_epsison= 1e-6;
  96
  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;
 100
 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));
 106
 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;
 113
 114     sigma_bd2_l += b * d2 / l3;
 115   }
 116
 117
 118
 119 static gsl_multimin_fminimizer *minimiser;
 120
 121
 122
 123 {
 124   gsl_multimin_fminimizer_alloc(
 125