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    *  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.  */
  66
  67 static double hypotD(const double p[], const double q[]) {
  68   int k;
  69   double pq[D3];
  70   gsl_vector v;
  71
  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 */
  77
  78   return gsl_blas_snrm2(&v);
  79 }
  80
  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
  86
  87 static const l3_epsison= 1e-6;
  88
  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;
  92
  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));
  98
  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;
 105
 106     sigma_bd2_l += b * d2 / l3;
 107   }
 108
 109
 110
 111 static gsl_multimin_fminimizer *minimiser;
 112
 113
 114
 115 {
 116   gsl_multimin_fminimizer_alloc(
 117