new rim_twist_cost as yet unused

diff --git a/energy.c b/energy.c
index 72c6fd395d45c9047358d4c6fa53699fa9d1af2a..aefdd8beb6f6637c86971d3099c0e5a34bb3689d 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -497,6 +497,46 @@ double noncircular_rim_cost(const Vertices vertices, int section) {
   return cost;
 }
 
+/*---------- rim contact angle rotation ----------*/
+
+/*
+ *
+ *               vq ---- vs
+ *              /       /
+ *           e0/       /e1
+ *           vp 0--- vr
+ *
+ */
+
+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);
+
+      total_cost += cost;
+    }
+  }
+
+  return total_cost;
+}
+
 /*---------- overly sharp edge cost ----------*/
 
   /*

diff --git a/minimise.h b/minimise.h
index f4eddfb222bfe4bb165475bffcb4d7d95e7300e3..4999fd04790bab482f69acba42fa46d74814addc 100644 (file)
--- a/minimise.h
+++ b/minimise.h
@@ -27,6 +27,7 @@ double noncircular_rim_cost(const Vertices vertices, int section);
 double edge_length_variation_cost(const Vertices vertices, int section);
 double prop_edge_length_variation_cost(const Vertices vertices, int section);
 double rim_proximity_cost(const Vertices vertices, int section);
+double rim_twist_cost(const Vertices vertices, int section);
 double edge_angle_cost(const Vertices vertices, int section);
 double small_triangles_cost(const Vertices vertices, int section);
 double nonequilateral_triangles_cost(const Vertices vertices, int section);