get rim twist right

[moebius2.git] / energy.c

diff --git a/energy.c b/energy.c
index aefdd8beb6f6637c86971d3099c0e5a34bb3689d..aa2dd47fa4364eed4fc5c616ca2336ba9573eee9 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -7,7 +7,11 @@
 #include "mgraph.h"
 #include "parallel.h"
 
-double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
+double vertex_mean_edge_lengths[N];
+
+static double vertex_areas[N];
+static double edge_lengths[N][V6];
+static double rim_vertex_angles[N];
 
 static double best_energy= DBL_MAX;
 
@@ -29,7 +33,8 @@ typedef struct {
 
 static PreComputation *const precomps[]= {
   compute_edge_lengths,
-  compute_vertex_areas
+  compute_vertex_areas,
+  compute_rim_twist_angles
 };
 
 static const CostContribution costs[]= {
@@ -45,10 +50,11 @@ static const CostContribution costs[]= {
 #endif
 
 #if XBITS==4
-#define STOP_EPSILON 1.2e-3
-    COST(  3e3,   vertex_displacement_cost)
-    COST(  3e3,   vertex_edgewise_displ_cost)
-    COST( 0.2e3,  rim_proximity_cost)
+#define STOP_EPSILON 5e-3
+    COST(  3e4,   vertex_displacement_cost) // NB this is probably wrong now
+    COST(  3e4,   vertex_edgewise_displ_cost) // we have changed the power
+    COST( 2e2,  rim_proximity_cost)
+    COST( 1e4,  rim_twist_cost)
     COST(  1e12,   noncircular_rim_cost)
     COST(  10e1,   nonequilateral_triangles_cost)
 //    COST(  1e1,   small_triangles_cost)
@@ -57,12 +63,13 @@ static const CostContribution costs[]= {
 #endif
 
 #if XBITS==5
-#define STOP_EPSILON 1.2e-4
-    COST(  3e3,   vertex_displacement_cost)
-    COST(  3e3,   vertex_edgewise_displ_cost)
+#define STOP_EPSILON 7e-4
+    COST(  3e4,   vertex_displacement_cost)
+    COST(  3e4,   vertex_edgewise_displ_cost)
     COST(  2e-1,  rim_proximity_cost)
+    COST(  3e3,  rim_twist_cost)
     COST(  1e12,   noncircular_rim_cost)
-    COST(  10e1,   nonequilateral_triangles_cost)
+    COST(   3e2,   nonequilateral_triangles_cost)
 //    COST(  1e1,   small_triangles_cost)
 //    COST(  1e6,   edge_angle_cost)
                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7)
@@ -70,9 +77,10 @@ static const CostContribution costs[]= {
 
 #if XBITS==6
 #define STOP_EPSILON 1.2e-4
-    COST(  3e3,   vertex_displacement_cost)
-    COST(  3e3,   vertex_edgewise_displ_cost)
-    COST( 0.02e0,  rim_proximity_cost)
+    COST(  3e4,   vertex_displacement_cost)
+    COST(  3e4,   vertex_edgewise_displ_cost)
+    COST( 2e-1,  rim_proximity_cost)
+    COST(  1e3,  rim_twist_cost)
     COST(  1e12,   noncircular_rim_cost)
     COST(  10e1,   nonequilateral_triangles_cost)
 //    COST(  1e1,   small_triangles_cost)
@@ -182,7 +190,7 @@ double compute_energy(const struct Vertices *vs) {
 
 static void addcost(double *energy, double tweight, double tcost, int pr) {
   double tenergy= tweight * tcost;
-  if (pr) printf(" %# e > %# e* |", tcost, tenergy);
+  if (pr) printf(/*" %# e >"*/ " %# e* |", /*tcost,*/ tenergy);
   *energy += tenergy;
 }
 
@@ -292,7 +300,7 @@ static double vertex_one_displ_cost(const double r[D3], const double s[D3],
 
 double vertex_displacement_cost(const Vertices vertices, int section) {
   const double inv_lambda= 1.0/1; //2;
-  const double delta= 4;
+  const double delta= 6;
 
   int si,e,qi,ri;
   double total_cost= 0;
@@ -317,7 +325,7 @@ double vertex_displacement_cost(const Vertices vertices, int section) {
 
 double vertex_edgewise_displ_cost(const Vertices vertices, int section) {
   const double inv_lambda= 1.0/1; //2;
-  const double delta= 4;
+  const double delta= 6;
 
   int pi,e,qi,ri,si, k;
   double m[D3];
@@ -463,7 +471,7 @@ static void find_nearest_oncircle(double oncircle[D3], const double p[D3]) {
 }
 
 double rim_proximity_cost(const Vertices vertices, int section) {
-  double oncircle[3], cost=0;
+  double oncircle[D3], cost=0;
   int v;
 
   FOR_VERTEX(v, OUTER) {
@@ -499,39 +507,33 @@ double noncircular_rim_cost(const Vertices vertices, int section) {
 
 /*---------- rim contact angle rotation ----------*/
 
-/*
- *
- *               vq ---- vs
- *              /       /
- *           e0/       /e1
- *           vp 0--- vr
- *
- */
+void compute_rim_twist_angles(const Vertices vertices, int section) {
+  double oncircle[D3], distance[D3];
+  int vpy,vpx,v,k;
+    
+  FOR_NEAR_RIM_VERTEX(vpy,vpx,v, 1,OUTER) {
+    find_nearest_oncircle(oncircle, vertices[v]);
+    /* we are interested in the angle subtended at the rim, from the
+     * rim's point of view. */
+    K distance[k]= vertices[v][k] - oncircle[k];
+    double distance_positive_z= distance[2];
+    double distance_radial_outwards= dotprod(distance, oncircle);
+    rim_vertex_angles[v]= atan2(distance_positive_z, distance_radial_outwards);
+  }
+}
 
 double rim_twist_cost(const Vertices vertices, int section) {
-  const double our_epsilon=1e-6;
-  int vpy,vpx,vpi, vqi,vri,vsi, e0,e1, k;
-  double total_cost= 0.0;
-  double pq[D3], rs[D3], pq_x_rs[D3];
-
-  FOR_RIM_VERTEX(vpy,vpx,vpi, OUTER) {
-    FOR_VPEDGE(e0) {
-      if (e0==0 || e0==3) continue;
-      vqi= EDGE_END2(vpi,e0);  if (vqi<0) continue;
-      vri= EDGE_END2(vpi,0);  assert(vri>=0);
-      e1= !vertices_span_join_p(vpi,vri) ? e0 : V6 - e0;
-      vsi= EDGE_END2(vri,e1);  assert(vsi>=0);
-
-      K pq[k]= -vertices[vpi][k] + vertices[vqi][k];
-      K rs[k]= -vertices[vri][k] + vertices[vsi][k];
-
-      xprod(pq_x_rs, pq,rs);
-      double magndiv= edge_lengths[vpi][e0] * edge_lengths[vri][e1];
-
-      double cost= magnD2(pq_x_rs) / (magndiv*magndiv + our_epsilon);
+  double total_cost= 0;
+  int vpy,vpx,v0,v1;
+  
+  FOR_NEAR_RIM_VERTEX(vpy,vpx,v0, 1,OUTER) {
+    v1= EDGE_END2(v0,0);  assert(v1!=0);
+    double delta= rim_vertex_angles[v0] - rim_vertex_angles[v1];
+    if (delta < M_PI) delta += 2*M_PI;
+    if (delta > M_PI) delta -= 2*M_PI;
 
-      total_cost += cost;
-    }
+    double cost= pow(delta, 4);
+    total_cost += cost;
   }
 
   return total_cost;