wip before increase power of bendings

[moebius2.git] / energy.c

diff --git a/energy.c b/energy.c
index 72c6fd395d45c9047358d4c6fa53699fa9d1af2a..8bc125af99f53df80d45f60e507cea337cf7efa7 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
+#define STOP_EPSILON 5e-3
     COST(  3e3,   vertex_displacement_cost)
     COST(  3e3,   vertex_edgewise_displ_cost)
     COST( 0.2e3,  rim_proximity_cost)
+    COST( 1e4,  rim_twist_cost)
     COST(  1e12,   noncircular_rim_cost)
     COST(  10e1,   nonequilateral_triangles_cost)
 //    COST(  1e1,   small_triangles_cost)
@@ -61,8 +67,9 @@ static const CostContribution costs[]= {
     COST(  3e3,   vertex_displacement_cost)
     COST(  3e3,   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)
@@ -73,6 +80,7 @@ static const CostContribution costs[]= {
     COST(  3e3,   vertex_displacement_cost)
     COST(  3e3,   vertex_edgewise_displ_cost)
     COST( 0.02e0,  rim_proximity_cost)
+    COST(  1e4,  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;
 }
 
@@ -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) {
@@ -497,6 +505,40 @@ double noncircular_rim_cost(const Vertices vertices, int section) {
   return cost;
 }
 
+/*---------- rim contact angle rotation ----------*/
+
+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[3];
+    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) {
+  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;
+
+    double cost= pow(delta, 4);
+    total_cost += cost;
+  }
+
+  return total_cost;
+}
+
 /*---------- overly sharp edge cost ----------*/
 
   /*