#include "Head.h"
#include <math.h>
/* Scalar times tensor */
static tri operator*(float e, tri f)
  {tri v; int i; for (i=0; i<3; i++) v.A[i] = e*f.A[i]; return v;}
/* Tensor times Vector */
static vec operator*(tri a, vec b) {vec c; int i; for (i=0; i<3; i++)
  c.a[i] = a.A[i]* b; return c;}

/* Push on a vertex. */
static void shove(int i, vec f) {u[i].mom = u[i].mom + f;}
static tri push(lt sh, tri or) {tri b=tp(or), v, w; int i, j;
  for (i=0; i<3; i++) w.A[i] = sh*b.A[i];
  for (i=0; i<3; i++) for (j=0; j<=i; j++)
    v.A[j].a[i] = v.A[i].a[j] = w.A[i]*w.A[j];
  return v;}
void physics() {/*float lam = -.2, mu = -.05;
   float mu_dt = lam*dt, lambda_dt = mu*dt;*/ int i; 
   if(damp){int i; for (i=0; i<n_vert; i++) u[i].mom = .99*u[i].mom;}
 /* Visit each zone. */ cycle++; for (i=0; i<xx; i++)
   {int n = s[i].v[0]; float compress =0.;
   vec a=cdvec(n,s[i].v[1]), b=cdvec(n,s[i].v[2]), c=cdvec(n,s[i].v[3]);
   tri strain = push(s[i].sh, ctri(a, b, c));
   {int j; for (j=0; j<3; j++) compress +=(strain.A[j].a[j] -= 1.);}
    compress *= lam[s[i].mat]*dt;
   {tri stress = mu[s[i].mat]*dt*strain; vec forth = zvec;
     int j; for (j=0; j<3; j++) stress.A[j].a[j] += compress; 
 {vec frc=stress*(a%b); forth=forth - frc; shove(s[i].v[3], frc);}
 {vec frc=stress*(b%c); forth=forth - frc; shove(s[i].v[1], frc);}
 {vec frc=stress*(c%a); forth=forth - frc; shove(s[i].v[2], frc);}
     shove(n, forth);}}
   {/* Visit each vertex. */ int i; for (i=0; i<n_vert; i++)
   {u[i].wh = u[i].wh + 1./u[i].mass*u[i].mom;}}}