init r0, r1

scontact/scontact.c

@@ -137,10 +137,6 @@ struct _PeriodicTriangle {
 };
 #endif
 
-
-
-
-
 static double dot(const double *x, const double *y) {
   #if DIMENSION == 3
   return x[0] * y[0] +  x[1] * y[1] + x[2] * y[2];
@@ -337,8 +333,7 @@ static int contact_init_disk_particle(Contact *c, ParticleProblem *p, size_t dis
   }
   c->type = PARTICLE_DISK;
   double basis[DIMENSION][DIMENSION];
-  double r0[DIMENSION], r1[DIMENSION];
-  double D = contact_update_disk_particle(c, p, x, basis, r0, r1);
+  double D = contact_update_disk_particle(c, p, x, basis, c->r0, c->r1);
   c->D0 = fmin(D, 0);
   return D < alert;
 }
@@ -412,8 +407,7 @@ static int contact_init_segment_particle(Contact *c, ParticleProblem *p, size_t
   }
   c->type = PARTICLE_SEGMENT;
   double basis[DIMENSION][DIMENSION];
-  double r0[DIMENSION],r1[DIMENSION];
-  double D = contact_update_segment_particle(c, p, x, basis, r0, r1);
+  double D = contact_update_segment_particle(c, p, x, basis, c->r0, c->r1);
   c->D0 = fmin(D, 0);
   return D < alert;
 }
@@ -804,8 +798,6 @@ void contact_tree_add_particle(ContactTree *tree, int id, const Contact *old_con
         if ((cold = findContactSorted(c, old_contacts, iold))) {
           for (int d = 0; d < DIMENSION; ++d) {
             c->reaction[d] = cold->reaction[d];
-            c->r0[d] = cold->r0[d];
-            c->r1[d] = cold->r1[d];
           }
         }
       }
@@ -855,8 +847,6 @@ static void contact_tree_gen_boundary_contact(ContactTree *tree, ContactType typ
       if ((cold = findContactSorted(c, old_contacts, iold))) {
         for (int d = 0; d < DIMENSION; ++d) {
           c->reaction[d] = cold->reaction[d];
-          c->r0[d] = cold->r0[d];
-          c->r1[d] = cold->r1[d];
         }
       }
     }
@@ -1023,10 +1013,10 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
   }
   // allow the pre-existing (from previous time steps) interpenetrations
   D -= c->D0;
-  double pold[DIMENSION] = {0};
+  double rold[DIMENSION] = {0};
   // un-apply the old contact
   for (int d = 0; d < DIMENSION; ++d) {
-    pold[d] = c->reaction[d]*dt;
+    rold[d] = c->reaction[d];
   }
   contact_apply(c, p, -dt);
   for (int d=0; d<DIMENSION; ++d) {
@@ -1044,9 +1034,9 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
   double wdet = wnn*wtt-wnt*wnt;
   double v[DIMENSION]= {0.}; // local free velocity in the contact basis
   contact_get_local_velocity(c, p, v, basis); // v[0] < 0 if particles get closer
-  double oldpn = pold[0]*basis[0][0] + pold[1]*basis[0][1];
-  double oldpt = pold[0]*basis[1][0] + pold[1]*basis[1][1];
-  v[0] = fmin(D/dt-oldpn*wnn-oldpt*wnt, 0.);
+  double oldrn = rold[0]*basis[0][0] + rold[1]*basis[0][1];
+  double oldrt = rold[0]*basis[1][0] + rold[1]*basis[1][1];
+  v[0] = fmin(D/dt-dt*(oldrn*wnn+oldrt*wnt), 0.);
   double pstick[2] = {
       -(+wtt*v[0]-wnt*v[1])/wdet, 
       -(-wnt*v[0]+wnn*v[1])/wdet
@@ -1067,22 +1057,22 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
       pt = pstick[1];
     }
   }
-  c->reaction[0] = (pn*basis[0][0] + pt*basis[1][0])/dt;
-  c->reaction[1] = (pn*basis[0][1] + pt*basis[1][1])/dt;
+  c->reaction[0] = (pn*basis[0][0]+pt*basis[1][0])/dt;
+  c->reaction[1] = (pn*basis[0][1]+pt*basis[1][1])/dt;
   
   // apply new contact and check convergence
-  double dvmax = 0;
+  double drmax = 0;
 #if DIMENSION == 2
-  double dp[DIMENSION];
+  double dr[DIMENSION];
   for (int d = 0; d <DIMENSION; ++d) {
-    dp[d] = c->reaction[d]*dt-pold[d];
+    dr[d] = c->reaction[d]-rold[d];
   }
-  dvmax = fmax(fabs(dp[0]*wnn+dp[1]*wnt),fabs(dp[0]*wnt+dp[1]*wtt));
+  drmax = fmax(fabs(dr[0]*wnn+dr[1]*wnt),fabs(dr[0]*wnt+dr[1]*wtt));
 #else
   //TOD0
 #endif
   contact_apply(c, p, dt);
-  return dvmax < tol/dt;
+  return drmax*dt*dt < tol;
 }
 
 void swap(int *a, int *b) {