- //double a2= vertex_areas[v2];
- double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
- double s= vertex_distances[v2];
- double s2= s*s + d2_epsilon;
- double sd2= s2 / d2;
+
+ double d= hypotD(v[v1],v[v2]);
+
+ int dist= distances[v1][v2];
+ assert(dist>=0);
+
+ double s= dist * meanedgelength * 0.03;
+
+ double enoughdistance= d - s;
+ if (enoughdistance > 1e-6) continue;
+
+ /* energy = 1/2 stiffness deviation^2
+ * where stiffness = 1/d
+ */
+
+ double cost= pow(enoughdistance,4);
+
+ //double s2= s*s + d2_epsilon;
+ //double sd2= s2 / d2;