better with more bendingness costs

[moebius2.git] / energy.c

diff --git a/energy.c b/energy.c
index 72c6fd395d45c9047358d4c6fa53699fa9d1af2a..3c1731a8e971c304c2d40bf222b741ad89011b4a 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)
+#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( 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)
@@ -57,22 +63,24 @@ 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 1e-5
+    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)
 #endif
 
 #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)
+#define STOP_EPSILON 1e-4
+    COST(  3e5,   vertex_displacement_cost)
+    COST(  3e5,   vertex_edgewise_displ_cost)
+    COST( 3e-2,  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;
 }
 
@@ -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) {
@@ -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 ----------*/
 
   /*