/*
* We try to find an optimal triangle grid
- *
- * Vertices in strip are numbered as follows:
- *
- * ___ X-2 ___ X-1 ___ 0 ___ 1 ___ 2 ___ 3 ___ 4 __
- * Y-1 Y-1 0 0 0 0 0
- * / \ / \ / \ / \ / \ / \ / \
- * / \ / \ / \ / \ / \ / \ / \
- * X-3 ___ X-2 ___ X-1 ___ 0 ___ 1 ___ 2 ___ 3 ___ 4
- * Y-2 Y-2 Y-2 1 1 1 1 1
- * \ / \ / \ / \ / \ / \ / \ /
- * \ / \ / \ / \ / \ / \ / \ /
- * ___ X-3 ___ X-2 ___ X-1 ___ 0 ___ 1 ___ 2 ___ 3 __
- * Y-3 Y-3 Y-3 2 2 2 2
- *
- * . . . . . . . . . . . . . . .
- *
- * X-4 ___ X-3 ___ X-2 ___ X-1 ___ 0 ___ 1 ___ 2 ___ 3
- * 1 1 1 1 Y-2 Y-2 Y-2 Y-2
- * \ / \ / \ / \ / \ / \ / \ /
- * \ / \ / \ / \ / \ / \ / \ /
- * ___ X-4 ___ X-3 ___ X-2 ___ X-1 ___ 0 ___ 1 ___ 2 __
- * 0 0 0 0 Y-1 Y-1 Y-1
- *
- * Node x,y for
- * 0 <= x < X x = distance along
- * 0 <= y < Y y = distance across
- *
- * Vertices are in reading order from diagram above ie x varies fastest.
- *
- * Y must be even. The actual location opposite (0,0) is (X-(Y-1)/2,0),
- * and likewise opposite (0,Y-1) is ((Y-1)/2,0).
- *
- * To label edges, we counte anticlockwise[*] from to-the-right:
- *
- * \2 /1
- * \ /
- * ___ 0 __
- * 3 1 0
- * / \
- * 4/ 5\
- *
- * [*] That is, in the direction represented as anticlockwise for
- * the vertices (0,*)..(4,*) in the diagram above; and of course
- * that is clockwise for the vertices (X-5,*)..(X-1,*). The
- * numbering has an actual discontinuity between (X-1,*) and (0,*).
- *
- * When we iterate over edges, we iterate first over vertices and then
- * over edges 0 to 2, disregarding edges 3 to 5.
*/
-#define XBITS 4
-#define X (1<<XBITS)
-#define YBITS 4
-#define Y (1<<YBITS)
-
-/* vertex number: 0000 | y | x
- * YBITS XBITS
- */
-
-#define N (X*Y)
-#define XMASK (X-1)
-#define YSHIFT XBITS
-#define Y1 (1 << YSHIFT)
-#define YMASK (Y-1 << YSHIFT)
-
/* Energy has the following parts right now:
*/
/*
* 2
* b d
* E = F . ----
- * vd, edge PQ vd l
+ * vd, edge PQ vd 3
+ * l
*
- * By symmetry, this calculation gives the same answer
- * with R and S exchanged. Looking at the projection in
- * the RMS plane:
+ * (The dimensions of this are those of F_vd.)
+ *
+ * By symmetry, this calculation gives the same answer with R and S
+ * exchanged. Looking at the projection in the RMS plane:
*
*
* S'
* because we want this symmetry and because we're happy to punish
* bending more than uneveness in the metric.
*
- * In practice to avoid division by zero we'll add epsilon to l and
- * the huge energy ought then to be sufficient for the model to
+ * In practice to avoid division by zero we'll add epsilon to l^3
+ * and the huge energy ought then to be sufficient for the model to
* avoid being close to R=S.
*/
-#define V6 6
-
-static int edge_end2(unsigned v1, int e) {
- /* The topology is equivalent to that of a square lattice with only
- * half of the diagonals. Ie, the result of shearing the triangular
- * lattice to make the lines of constant x vertical. This gives
- * these six directions:
- *
- * 2 1
- * | /
- * |/
- * 3--*--0
- * /|
- * / |
- * 4 5
- *
- * This also handily makes vertical the numbering discontinuity,
- * where the join happens.
- */
- static const unsigned dx[V6]= { +1, +1, 0, -1, -1, 0 },
- dy[V6]= { 0, +Y1, +Y1, 0, -Y1, -Y1 };
- unsigned x, y;
-
- y= (v1 & YMASK) + dy[e];
- if (y & ~YMASK) return -1;
-
- x= (v1 & XMASK) + dx[e];
- if (x & ~XMASK) {
- y= (Y-1)*Y1 - y;
- x &= XMASK;;
- }
- return x | y;
-}
-
-#define D3 3
-
-#define FOR_VERTEX(v) \
- for ((v)=0; (v)<N; (v)++)
-
-#define FOR_VPEDGE(v,e) \
- for ((e)=0; (e)<V6; (e)++)
-
-#define EDGE_END2 edge_end2
-
-#define FOR_VEDGE(v1,e,v2) \
- FOR_VPEDGE((v1),(e))
- if (((v2)= EDGE_END2((v1),(e))) < 0) ; else
-
-#define FOR_EDGE(v1,e,v2) \
- FOR_VERTEX((v1)) \
- FOR_VEDGE((v1),(e),(v2))
-
-#define FOR_COORD(k) \
- for ((k)=0; (k)<D3; (k)++)
-
-#define K FOR_COORD(k)
-
-static double hypotD(const double p[], const double q[]) {
+double hypotD(const double p[D3], const double q[D3]) {
int k;
double pq[D3];
gsl_vector v;
return gsl_blas_snrm2(&v);
}
+double hypotD2(const double p[D3], const double q[D3]) {
+ double d2= 0;
+ K d2= ffsqa(p[k] - q[k], d2);
+ return d2;
+}
+
#ifdef FP_FAST_FMA
# define ffsqa(factor,term) fma((factor),(factor),(term))
#else
-# define ffsqu(factor,term) ((factor)*(factor)+(term))
+# define ffsqa(factor,term) ((factor)*(factor)+(term))
#endif
+static const l3_epsison= 1e-6;
+
static double energy_function(const double vertices[N][D3]) {
int pi,e,qi,ri,si, k;
double m[D3], mprime[D3], b, d2, l, sigma_bd2_l;
K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
b= hypotD(vertices[pi], vertices[qi]);
- d2= 0; K d2= ffsqa(m[k] - mprime[k], d2);
+ d2= hypotD2(m, mprime);
l= hypotD(vertices[ri][k] - vertices[si][k]);
+ l3 = l*l*l + l3_epsilon;
- sigma_bd2_l += b * d2 / l;
+ sigma_bd2_l += b * d2 / l3;
}
~ian / moebius2.git / blobdiff