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 ----------*/
/*
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);