integration.h

Go to the documentation of this file.

#ifndef INTEGRATION_H
#define INTEGRATION_H

#include <particle.h>
#include <particledisp.h>
#include <particlesampler.h>
#include <typedefs.h>
#include <zpr.h>

// TODO rename integration to particleintbox
//   remove template F as d2toindex is used.
//   remove template D as it is part of the initial
//   setup and hence no point optimizing. 
//   The code becomes fixed, 
//   The distributions code should use inheritance 


// D is the distribution.
// F is the index function.
template< typename D, typename F >
class integration
{
  particledisp pd;
  D distribution;
  F fn; 
  cell C;
public:

  long stepcount;
  gltextmsg stepcountmsg;

  uint numParticles;
  bool uniform;
  bool printstate;
  float h;

  bool displayEk;

  bool LJ; 
  double LJpsi;
  double LJsigma;
  void LJcalc();

  // Histogram parameters
  double histlow;
  double histhigh;
  uint histn;

  particlesampler sampler;
  uint samplerstate;

  particle * pi;

  uint * pu;

  framerate<> fr;

  integration();

  void reset(int argc, char**argv);

  void draw();

  void bruteforce();

  void bruteforce
  (
    uintc * p,
    uintc sz
  );

  bool run;
  void togglestartstop();

  void step()
  {
    if (!uniform)
      bruteforce();
    else
      uniformgridtest();

    ++stepcount;

    if (samplerstate)
      ++sampler;
  }
  void steplarge();


  void uniformgridtest();

//  void steppersecupdate();

};


//---------------------------------------------------------
// Implementation

/*
template< typename D, typename F >
void integration<D,F>::steppersecupdate()
{
  static int t0;
  static uint stepcount0;

  int t1 = glutGet(GLUT_ELAPSED_TIME);
  double t = 1000.0/(t1-t0);
  t0=t1;

  uint steppersec = t*(stepcount-stepcount0);
  stepcount0 = stepcount;

  steppersecmsg.updatevalue(steppersec,"step/s");
}
*/





template< typename D, typename F >
void integration<D,F>::uniformgridtest()
{
  boxcollision const & box(distribution.box);
  for (uint i=0; i<numParticles; ++i)
  {
    pi[i].step(h);
    box.walls(pi[i],h);
  }

  C.eval(fn);

if (true&&printstate)
{
uint W(fn.W);
uint count;
int idx=0;
cout << endl;
for (uint i=0; i<W; ++i)
{
  for (uint j=0; j<W; ++j)
  {
   
    count=0;
    for (uint k=0; k<numParticles; ++k)
    {
      if (idx==fn.index(k))
      {
        cout << k << ',';
        ++count;
      }
    }
    if (count==0)
      cout << "*";
    cout << " ";
    ++idx;
  }
  cout << endl;
}
}

  uint sz;
  for (uint i=0; i<numParticles; ++i)
  {
    C.getsurroundingparticles(pu,sz,fn,i);

//cout << SHOW(i) << "   " << SHOW(sz) << "  ";
//print(pu,pu+sz);

    bruteforce(pu,sz);
  }


if (true&&printstate)
{
cout << endl;
for (uint i=0; i<numParticles; ++i)
{
  cout << fn.index(i) << " pi[" << i << "] ";
  C.getsurroundingparticles(pu,sz,fn,i);
  cout << "  ";
  for (uint k=0; k<sz; ++k)
    cout << pu[k] << " ";
  cout << "        ";
  cout << pi[i];
cout << endl;
}
}


if (true&&printstate)
{
uint W(fn.W);
uint count;
int idx=0;
cout << endl;
for (uint i=0; i<W; ++i)
{
  for (uint j=0; j<W; ++j)
  {
   
    count=0;
    for (uint k=0; k<numParticles; ++k)
    {
      if (idx==fn.index(k))
      {
        cout << k << ',';
        ++count;
      }
    }
    if (count==0)
      cout << "*";
    cout << " ";
    ++idx;
  }
  cout << endl;
}
}


}


template< typename D, typename F >
void integration<D,F>::steplarge()
{
  for (uint i=0; i<100; ++i)
    step();
}
  

template< typename D, typename F >
void integration<D,F>::togglestartstop()
{
  if (run)
    run=false;
  else
    run=true;
}

template< typename D, typename F >
void integration<D,F>::LJcalc()
{
  double rx;
  double ry;

  double r;
  double z;
  double f;
  double fx;
  double fy;

  for (uint i=0; i<numParticles; ++i)
    pi[i].acczero();

//cout << "Before" << endl;
//  for (uint i=0; i<numParticles; ++i)
//    cout << SHOW(pi[i]) << endl;
//cout << endl;

  for (uint i=0; i+1<numParticles; ++i)
  {
    for (uint k=i+1; k<numParticles; ++k)
    {
      rx = pi[i].pos[0] - pi[k].pos[0];
      ry = pi[i].pos[1] - pi[k].pos[1];
    
      r = sqrt(rx*rx+ry*ry);
//cout << SHOW(r) << " ";
      z = LJsigma/r;
//cout << SHOW(z) << " ";
      z *= z;
      z *= z*z;

      // Force divided by r.
      f = LJpsi*24.0/(r*r)*z*(2.0*z-1.0);

if (abs(f)>10)
cout << SHOW(f) << endl;

//cout << SHOW(f) << " ";

      fx = f*rx/r;
      fy = f*ry/r;

//cout << SHOW(fx) << " ";
//cout << SHOW(fy) << " ";
//cout << endl;

      pi[i].acc[0] -= fx;
      pi[i].acc[1] -= fy;

      pi[k].acc[0] += fx;
      pi[k].acc[1] += fy;
    }

  }

//  for (uint i=0; i<numParticles; ++i)
//    cout << SHOW(pi[i]) << endl;

}

template< typename D, typename F >
void integration<D,F>::bruteforce
(
  uintc * x,
  uintc sz
)
{
  boxcollision const & box(distribution.box);

//cout << SHOW(sz) "  ";
//print(p,p+sz);


  uint a;
  uint b;
  for (uint i=0; i+1<sz; ++i)
  {
    for (uint k=i+1; k<sz; ++k)
    {
      a = x[i];
      b = x[k];

      if (pi[a].isIntersecting(pi[b]))
      {
        pi[a].step(-h);
        pi[b].step(-h);

        pi[a].collisionDynamics(pi[b]);

        pi[a].step(h);
        pi[b].step(h);

        box.walls(pi[a],h);
        box.walls(pi[b],h);
      }
    }
  }
}


template< typename D, typename F >
void integration<D,F>::bruteforce()
{
  boxcollision const & box(distribution.box);
  for (uint i=0; i<numParticles; ++i)
  {
    pi[i].step(h);
    box.walls(pi[i],h);
  }

  if (LJ)
    LJcalc();
  
  for (uint i=0; i+1<numParticles; ++i)
  {
    for (uint k=i+1; k<numParticles; ++k)
    {

      if (pi[i].isIntersecting(pi[k]))
      {
        pi[i].step(-h);
        pi[k].step(-h);

        pi[i].collisionDynamics(pi[k]);

        pi[i].step(h);
        pi[k].step(h);

if (!LJ)
{
        box.walls(pi[i],h);
        box.walls(pi[k],h);
}
      }
    }
  }

    

#ifdef DEBUG_PART
if (false)
{
cout << endl;
for (uint i=0; i<numParticles; ++i)
  cout << pi[i] << endl;
}
#endif

}


template< typename D, typename F >
void integration<D,F>::draw()
{
  // Display the particles

  glPushAttrib(GL_CURRENT_BIT);
  glPushAttrib(GL_LIGHTING_BIT);
  glColor3f(1.0,0.0,0.0);
  pd.draw(pi[0]);
  glPopAttrib();
  glPopAttrib();

  for (uint i=1; i+1<numParticles; ++i)
  {
    if ((i%10)!=0)
      pd.draw(pi[i]);
    else
    {
      glPushAttrib(GL_CURRENT_BIT);
      glPushAttrib(GL_LIGHTING_BIT);
      glColor3f(0.1,0.2,1.0);
      pd.draw(pi[i]);
      glPopAttrib();
      glPopAttrib();
    }
  }

  glPushAttrib(GL_CURRENT_BIT);
  glPushAttrib(GL_LIGHTING_BIT);
  glColor3f(0.0,1.0,0.0);
  pd.draw(pi[numParticles-1]);
  glPopAttrib();
  glPopAttrib();


  // Verify energy in system.
  if (displayEk)
  {

    static int ekcounter;
    ++ekcounter;
    if ((ekcounter%10)==0)
    {

      double sum(0.0);
      double energy;
      for (uint i=0; i<numParticles; ++i)
      {
        particle const & p(pi[i]);
        energy = p.vel[0] * p.vel[0] + p.vel[1] * p.vel[1];
     
        sum += energy;
      }
      cout << "Ek*2/m=" << sum << endl;
    }
  }

  // Draw Box
  glPushAttrib(GL_CURRENT_BIT);
  glPushAttrib(GL_LIGHTING_BIT);
  glDisable(GL_LIGHTING);
  glColor3f(.6,.03,.03);
  glBegin(GL_LINES);

  doublec x0 = distribution.box.x0;
  doublec x1 = distribution.box.x1;
  doublec y0 = distribution.box.y0;
  doublec y1 = distribution.box.y1;

  glVertex3f(x0,y0,0.0);
  glVertex3f(x1,y0,0.0);
  glVertex3f(x1,y0,0.0);
  glVertex3f(x1,y1,0.0);
  glVertex3f(x1,y1,0.0);
  glVertex3f(x0,y1,0.0);
  glVertex3f(x0,y1,0.0);
  glVertex3f(x0,y0,0.0);

  glEnd();
  glPopAttrib();
  glPopAttrib();

  // Update Counters

  fr.display();
  stepcountmsg.display();
}




template< typename D, typename F >
integration<D,F>::integration()
  : pd(gluNewQuadric(),10,3),
  distribution 
  (
    0.02,  //radius
    5.0,  //vmax
    boxcollision(0.0,10.0,0.0,10.0)  //box dimensions
  ),
  fn
  (
    distribution.box
  ),
  stepcount(0),
  stepcountmsg("0",0.5,-0.8,0.0,0.0,1.0),
  histlow(0.0),
  histhigh(7.0),
  histn(81),
  sampler
  (
    (uintc &)samplerstate,
    (uintc)20,     // Sample every x steps all particle speeds
    (uintc)5,      // Inner loop: sample every x2 steps
    histogram<double>
    (
      histlow,   // Lowest speed
      histhigh, // Highest speed
      histn      // Intervals
    ),
    string("data.txt"),
    (particle const * &)pi,
    (uintc & )numParticles
  )
/*
  (
    samplerstate,
    (uintc)20,     // Sample every x steps all particle speeds
    (uintc)5,      // Inner loop: sample every x2 steps
    histogram<double>
    (
      histlow,   // Lowest speed
      histhigh, // Highest speed
      histn      // Intervals
    ),
    string("data.txt"),
    (particle const *)pi,
    numParticles
  )
*/
{
  uniform=false;
  printstate=false;
  pi=0;
  numParticles=10;
  h=0.01;
  //zpr=true;
  run=false;
  displayEk=false;
  LJ=false;
  LJsigma=1.0;
  LJpsi=1.0;
}


template< typename D, typename F >
void integration<D,F>::reset(int argc, char**argv)
{
  commandline c(argc,argv);
  c.mapvar(numParticles,"numParticles");
  if(pi!=0)
    delete[] pi;
  fn.pi = pi = new particle[numParticles];

  C.reset(numParticles);


// Test Code
pu = new uint [numParticles];
for (uint i=0; i<numParticles; ++i)
  pu[i]=i;

  c.mapvar(uniform,"uniform");
  c.mapvar(distribution.vmax,"vmax");
  c.mapvar(distribution.radius,"radius");
  c.mapvar(h,"h");
  c.mapvar(printstate,"printstate");

// <TODO> put this option back in.
//  c.mapvar(zpr,"zpr");
//  if (zpr)
//    zprInit();

  c.mapvar(distribution.box.x0,"x0");
  c.mapvar(distribution.box.x1,"x1");
  c.mapvar(distribution.box.y0,"y0");
  c.mapvar(distribution.box.y1,"y1");

  c.mapvar(displayEk,"displayEk");

  c.mapvar(samplerstate,"samplerstate");

  c.mapvar(LJpsi,"LJpsi");
  c.mapvar(LJsigma,"LJsigma");
  c.mapvar(LJ,"LJ");
  if (LJ)
    uniform=false;

  c.mapvar(histlow,"histlow");
  c.mapvar(histhigh,"histhigh");
  c.mapvar(histn,"histn");
  
  sampler.set(histogram<double>(histlow,histhigh,histn));

  c.enablehelp(cout);

  fn.reset();

  // Initial conditions: no intersecting particles.
  for (uint i=0; i<numParticles; ++i)
  {
    distribution.eval(pi[i]);
    bool flagnoint(true);
    for (uint k=0; (k<i)&&flagnoint; ++k)
    {
      if( pi[i].isIntersecting( pi[k] ) )
      {
        flagnoint=false;
        --i;
      }
    }
  } 

}





#endif