#include "nbin_ssa.h"
#include "atom.h"
#include "update.h"
#include "group.h"
#include "domain.h"
using namespace LAMMPS_NS;
NBinSSA::NBinSSA(LAMMPS *lmp) : NBinStandard(lmp)
{
for (int i = 0; i < 8; i++) {
gairhead_ssa[i] = -1;
}
}
NBinSSA::~NBinSSA()
{
}
void NBinSSA::bin_atoms()
{
int i,ibin;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
if (includegroup) nlocal = atom->nfirst;
double **x = atom->x;
int *mask = atom->mask;
int xbin,ybin,zbin;
last_bin = update->ntimestep;
bboxlo_[0] = bboxlo[0]; bboxlo_[1] = bboxlo[1]; bboxlo_[2] = bboxlo[2];
bboxhi_[0] = bboxhi[0]; bboxhi_[1] = bboxhi[1]; bboxhi_[2] = bboxhi[2];
for (i = 0; i < 8; i++) {
gairhead_ssa[i] = -1;
}
for (i = 0; i < mbins; i++) {
binhead[i] = -1;
}
if (includegroup) {
int bitmask = group->bitmask[includegroup];
int nowned = atom->nlocal; for (i = nall-1; i >= nowned; i--) {
ibin = coord2ssaAIR(x[i]);
if (ibin < 1) continue; if (mask[i] & bitmask) {
bins[i] = gairhead_ssa[ibin];
gairhead_ssa[ibin] = i;
}
}
} else {
for (i = nall-1; i >= nlocal; i--) {
ibin = coord2ssaAIR(x[i]);
if (ibin < 1) continue; bins[i] = gairhead_ssa[ibin];
gairhead_ssa[ibin] = i;
}
}
for (i = nlocal-1; i >= 0; i--) {
ibin = coord2bin(x[i][0], x[i][1], x[i][2], xbin, ybin, zbin);
if (xbin < lbinxlo) lbinxlo = xbin;
if (xbin >= lbinxhi) lbinxhi = xbin + 1;
if (ybin < lbinylo) lbinylo = ybin;
if (ybin >= lbinyhi) lbinyhi = ybin + 1;
if (zbin < lbinzlo) lbinzlo = zbin;
if (zbin >= lbinzhi) lbinzhi = zbin + 1;
bins[i] = binhead[ibin];
binhead[ibin] = i;
}
}
void NBinSSA::bin_atoms_setup(int nall)
{
NBinStandard::bin_atoms_setup(nall);
lbinxlo = mbinx - 1; lbinylo = mbiny - 1; lbinzlo = mbinz - 1; lbinxhi = 0; lbinyhi = 0; lbinzhi = 0; }
bigint NBinSSA::memory_usage()
{
bigint bytes = NBinStandard::memory_usage();
return bytes;
}
int NBinSSA::coord2ssaAIR(const double *x)
{
int ix, iy, iz;
ix = iy = iz = 0;
if (x[2] < domain->sublo[2]) iz = -1;
if (x[2] >= domain->subhi[2]) iz = 1;
if (x[1] < domain->sublo[1]) iy = -1;
if (x[1] >= domain->subhi[1]) iy = 1;
if (x[0] < domain->sublo[0]) ix = -1;
if (x[0] >= domain->subhi[0]) ix = 1;
if(iz < 0){
return -1;
} else if(iz == 0){
if (iy<0) return -1; if ((iy==0) && (ix<0) ) return -1; if ((iy==0) && (ix==0)) return 0; if ((iy==0) && (ix>0) ) return 2; if ((iy>0) && (ix==0)) return 1; if ((iy>0) && (ix!=0)) return 3; } else { if((ix==0) && (iy==0)) return 4; if((ix==0) && (iy!=0)) return 5; if((ix!=0) && (iy==0)) return 6; if((ix!=0) && (iy!=0)) return 7; }
return -2;
}