cell.h

Go to the documentation of this file.

#ifndef CELL_H
#define CELL_H

#include <cassert>
#include <vector>
#include <cmath>
#include <iostream>
using namespace std;

#include <typedefs.h>


//#define DEBUG_CELL

//
//  F Interface
//
//  uint F::W;
//
//  int const F::index( P const & ) const;
//
//  void getsurroundingcells
//  (
//    uint * & _nb,
//    uint & sz,
//    uintc k
//  );
//


class cell
{
public:

  uint * count;  // 
  uint * index;  
  uint * vi;
  uint * fi;    

  uint W;
  uint WW;
  uint N;
public:

  cell();
  cell(uintc N_);
  ~cell();

  void reset(uintc N_);

  template <typename F>
  void eval(F & f);

  template <typename F>
  void getsurroundingparticles
  (
    uint * const nb,
    uint & sz,
    F & f,
    uintc pti
  );

};


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

template< typename T >
void print(T const * v, T const * vend)
{
  for (;v!=vend; ++v)
    cout << *v << " "; 

  cout << endl;
}

template <typename F>
void cell::eval(F & f)
{
  f.W = W;
assert(WW==W*W);

  for (uint i=0; i<WW; ++i)
    assert(count[i]==0);

  for (uint i=0; i<N; ++i)
  {
    fi[i] = f.index(i);
    ++count[fi[i]];
  }

  index[0] = 0;
  for (uint i=1; i<WW; ++i)
    index[i] = index[i-1]+count[i-1];


#ifdef DEBUG_CELL 
cout << "fi     ";
print(fi,fi+N);

cout << "count  ";
print(count,count+WW);

cout << "index  ";
print(index,index+WW);
#endif

  uint k;
  for (uint i=0; i<N; ++i)
  {
    k = fi[i];
    assert(k<WW);
    vi[index[k] + count[k] - 1] = i;
    --count[k];
  }

  // count[] should be all zero.
  for (uint i=0; i<WW; ++i)
    assert(count[i]==0);

#ifdef DEBUG_CELL
cout << "vi      ";
print(vi,vi+N);

cout << "count  ";
print(count,count+WW);
#endif

}

template <typename F>
void cell::getsurroundingparticles
(
  uint * const nb,
  uint & sz,
  F & f,
  uintc pti
)
{
  uint * ci;
  uint n(0);
  uint c;
  uint c0;
  uint z(0);
  f.getsurroundingcells(ci,n,f.index(pti));
  for (uint i=0; i<n; ++i)
  {
    // Determine how many particles in cell k
    c0 = index[ci[i]];
    c = index[ci[i]+1]-c0;

    for (uint k=0; k<c; ++k)
      nb[z++] = vi[c0+k];
  }
  sz = z;
}





#endif