branches/dynamic/src/dynamic.cc

   1 // Voro++, a 3D cell-based Voronoi library
   2 //
   3 // Author   : Chris H. Rycroft (LBL / UC Berkeley)
   4 // Email    : chr@alum.mit.edu
   5 // Date     : July 1st 2008
   6
   7 /** \file dynamic.hh
   8  * \brief Header file for the dynamic extension classes, which add
   9  * functionality for a variety of dynamic particle motions. */
  10
  11 #include "dynamic.hh"
  12
  13 template<class r_option>
  14 container_dynamic_base<r_option>::container_dynamic_base(fpoint xa,fpoint xb,fpoint ya,
  15                 fpoint yb,fpoint za,fpoint zb,int xn,int yn,int zn,const bool xper,const bool yper,
  16                 const bool zper,int memi) : container_base<r_option> (xa,xb,ya,yb,za,zb,xn,yn,zn,xper,yper,zper,memi), v_inter(ve) {
  17         gh=new int[nxyz];
  18         ve=new fpoint*[nxyz];
  19         for(int i=0;i<nxyz;i++) ve[i]=new fpoint[3*memi];
  20 }
  21
  22 template<class r_option>
  23 container_dynamic_base<r_option>::~container_dynamic_base() {
  24         for(int i=0;i<nxyz;i++) delete [] ve[i];
  25         delete [] ve;
  26         delete [] gh;
  27 }
  28
  29 template<class r_option>
  30 void container_dynamic_base<r_option>::wall_diagnostic() {
  31         fpoint x,y,z,dx,dy,dz;
  32         int i,j,k,ijk=0,l,w;
  33         for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
  34                 for(l=0;l<co[ijk];l++) {
  35                         x=p[ijk][sz*l];y=p[ijk][sz*l+1];z=p[ijk][sz*l+2];
  36                         cout << id[ijk][l] << " " << x << " " << y << " " << z;
  37                         radius.print(cout,ijk,l);
  38                         for(w=0;w<wall_number;w++) {
  39                                 walls[w]->min_distance(x,y,z,dx,dy,dz);
  40                                 cout << " " << dx << " " << dy << " " << dz << "\n";
  41                         }
  42                 }
  43         }
  44 }
  45
  46 template<class r_option>
  47 int container_dynamic_base<r_option>::count(fpoint x,fpoint y,fpoint z,fpoint r) {
  48         fpoint xx,yy,zz,px,py,pz;
  49         int cou=0,s,q;
  50         voropp_loop l1(this);
  51         s=l1.init(x,y,z,r,px,py,pz);r*=r;
  52         do {
  53                 for(q=0;q<co[s];q++) {
  54                         xx=p[s][sz*q]+px-x;yy=p[s][sz*q+1]+py-y;zz=p[s][sz*q+2]+pz-z;
  55                         if (xx*xx+yy*yy+zz*zz<r) cou++;
  56                 }
  57         } while((s=l1.inc(px,py,pz))!=-1);
  58         return cou;
  59 }
  60
  61 template<class r_option>
  62 void container_dynamic_base<r_option>::spot(fpoint cx,fpoint cy,fpoint cz,fpoint dx,fpoint dy,fpoint dz,fpoint rad) {
  63         velocity_constant vcl(dx,dy,dz);
  64         local_move(vcl,cx,cy,cz,rad);
  65 }
  66
  67 template<class r_option>
  68 void container_dynamic_base<r_option>::gauss_spot(fpoint cx,fpoint cy,fpoint cz,fpoint dx,fpoint dy,fpoint dz,fpoint dec,fpoint rad) {
  69         velocity_gaussian vcl(cx,cy,cz,dx,dy,dz,dec);
  70         local_move(vcl,cx,cy,cz,rad);
  71 }
  72
  73 template<class r_option>
  74 template<class v_class>
  75 void container_dynamic_base<r_option>::local_move(v_class &vcl,fpoint cx,fpoint cy,fpoint cz,fpoint rad) {
  76         int l,ll,s,ni,nj,nk;
  77         fpoint x,y,z,ex,ey,ez,px,py,pz;
  78         voropp_loop l1(this);
  79         s=l1.init(cx,cy,cz,rad,px,py,pz);rad*=rad;
  80         do {gh[s]=co[s];} while ((s=l1.inc(px,py,pz))!=-1);
  81
  82         s=l1.reset(px,py,pz);
  83         do {
  84                 l=0;
  85                 while(l<gh[s]) {
  86                         x=p[s][sz*l]+px;ex=x-cx;
  87                         y=p[s][sz*l+1]+py;ey=y-cy;
  88                         z=p[s][sz*l+2]+pz;ez=z-cz;
  89                         if(ex*ex+ey*ey+ez*ez>=rad) {l++;continue;}
  90                         vcl.vel(s,l,x,y,z);
  91
  92                         ni=step_int((x-ax)*xsp);
  93                         nj=step_int((y-ay)*ysp);
  94                         nk=step_int((z-az)*zsp);
  95
  96                         if(ni==l1.ip&&nj==l1.jp&&nk==l1.kp) {
  97                                 p[s][sz*l]=x;
  98                                 p[s][sz*l+1]=y;
  99                                 p[s][sz*l+2]=z;
 100                                 l++;
 101                         } else {
 102                                 if((xperiodic||(ni>=0&&ni<nx))
 103                                  &&(yperiodic||(nj>=0&&nj<ny))
 104                                  &&(zperiodic||(nk>=0&&nk<nz))) {
 105                                         if(xperiodic) {x-=step_div(ni,nx)*(bx-ax);ni=step_mod(ni,nx);}
 106                                         if(yperiodic) {y-=step_div(nj,ny)*(by-ay);nj=step_mod(nj,ny);}
 107                                         if(zperiodic) {z-=step_div(nk,nz)*(bz-az);nk=step_mod(nk,nz);}
 108                                         ni+=nx*(nj+ny*nk);
 109                                         if(co[ni]==mem[ni]) add_particle_memory(ni);
 110                                         id[ni][co[ni]]=id[s][l];
 111                                         p[ni][co[ni]*sz]=x;
 112                                         p[ni][co[ni]*sz+1]=y;
 113                                         p[ni][co[ni]*sz+2]=z;
 114                                         if(sz==4) p[ni][co[ni]*sz+3]=p[s][l*sz+3];
 115                                         co[ni]++;
 116                                 }
 117                                 co[s]--;
 118                                 id[s][l]=id[s][co[s]];
 119                                 for(ll=0;ll<sz;ll++) p[s][sz*l+ll]=p[s][sz*co[s]+ll];
 120                                 if(co[s]+1==gh[s]) {
 121                                         if (vcl.track_ve) {
 122                                                 ve[s][3*l]=ve[s][3*co[s]];
 123                                                 ve[s][3*l+1]=ve[s][3*co[s]+1];
 124                                                 ve[s][3*l+2]=ve[s][3*co[s]+2];
 125                                         }
 126                                         gh[s]--;
 127                                 } else l++;
 128                         }
 129                 }
 130         } while ((s=l1.inc(px,py,pz))!=-1);
 131 }
 132
 133 template<class r_option>
 134 void container_dynamic_base<r_option>::relax(fpoint cx,fpoint cy,fpoint cz,fpoint rad,fpoint alpha) {
 135         int l,q,s1,s2;
 136         fpoint x,y,z,ox,oy,oz,px,py,pz,ex,ey,ez;
 137         fpoint dx,dy,dz,dd,rr,radsq=rad*rad;
 138         voropp_loop l1(this),l2(this);
 139         s1=l1.init(cx,cy,cz,rad,ox,oy,oz);
 140         do {
 141                 for(l=0;l<3*co[s1];l++) ve[s1][l]=0;
 142         } while ((s1=l1.inc(ox,oy,oz))!=-1);
 143
 144         s1=l1.reset(ox,oy,oz);
 145         do {
 146                 for(l=0;l<co[s1];l++) {
 147                         ex=p[s1][sz*l];if(xperiodic) ex+=ox;dx=ex-cx;
 148                         ey=p[s1][sz*l+1];if(yperiodic) ey+=oy;dy=ey-cy;
 149                         ez=p[s1][sz*l+2];if(zperiodic) ez+=oz;dz=ez-cz;
 150                         dd=dx*dx+dy*dy+dz*dz;
 151                         if(dd>=radsq) {l++;continue;}
 152
 153                         s2=l2.init(ex,ey,ez,1,px,py,pz);
 154                         while(!l2.reached(l1.i,l1.j,l1.k)) {
 155                                 for(q=0;q<co[s2];q++) {
 156                                         x=p[s2][sz*q]-ex;if(xperiodic) x+=px;
 157                                         y=p[s2][sz*q+1]-ey;if(yperiodic) y+=py;
 158                                         z=p[s2][sz*q+2]-ez;if(zperiodic) z+=pz;
 159                                         rr=x*x+y*y+z*z;
 160                                         if (rr<1) {
 161                                                 if(dd+rr+2*(dx*x+dy*y+dz*z)<radsq) {
 162                                                         rr=alpha*(0.5-0.5/sqrt(rr));
 163                                                         ve[s2][3*q]-=x*rr;
 164                                                         ve[s2][3*q+1]-=y*rr;
 165                                                         ve[s2][3*q+2]-=z*rr;
 166                                                 } else rr=alpha*(1-1/sqrt(rr));
 167                                                 ve[s1][3*l]+=x*rr;
 168                                                 ve[s1][3*l+1]+=y*rr;
 169                                                 ve[s1][3*l+2]+=z*rr;
 170                                         }
 171                                 }
 172                                 s2=l2.inc(px,py,pz);
 173                         }
 174                         for(q=0;q<l;q++) {
 175                                 x=p[s2][sz*q]-ex;if(xperiodic) x+=px;
 176                                 y=p[s2][sz*q+1]-ey;if(yperiodic) y+=py;
 177                                 z=p[s2][sz*q+2]-ez;if(zperiodic) z+=pz;
 178                                 rr=x*x+y*y+z*z;
 179                                 if (rr<1) {
 180                                         if(dd+rr+2*(dx*x+dy*y+dz*z)<radsq) {
 181                                                 rr=alpha*(0.5-0.5/sqrt(rr));
 182                                                 ve[s2][3*q]-=x*rr;
 183                                                 ve[s2][3*q+1]-=y*rr;
 184                                                 ve[s2][3*q+2]-=z*rr;
 185                                         } else rr=alpha*(1-1/sqrt(rr));
 186                                         ve[s1][3*l]+=x*rr;
 187                                         ve[s1][3*l+1]+=y*rr;
 188                                         ve[s1][3*l+2]+=z*rr;
 189                                 }
 190                         }
 191                         q++;
 192                         while(q<co[s2]) {
 193                                 x=p[s2][sz*q]-ex;if(xperiodic) x+=px;
 194                                 y=p[s2][sz*q+1]-ey;if(yperiodic) y+=py;
 195                                 z=p[s2][sz*q+2]-ez;if(zperiodic) z+=pz;
 196                                 rr=x*x+y*y+z*z;
 197                                 if (rr<1) {
 198                                         if(dd+rr+2*(dx*x+dy*y+dz*z)>=radsq) {
 199                                                 rr=alpha*(1-1/sqrt(rr));
 200                                                 ve[s1][3*l]+=x*rr;
 201                                                 ve[s1][3*l+1]+=y*rr;
 202                                                 ve[s1][3*l+2]+=z*rr;
 203                                         }
 204                                 }
 205                                 q++;
 206                         }
 207                         while((s2=l2.inc(px,py,pz))!=-1) {
 208                                 for(q=0;q<co[s2];q++) {
 209                                         x=p[s2][sz*q]-ex;if(xperiodic) x+=px;
 210                                         y=p[s2][sz*q+1]-ey;if(yperiodic) y+=py;
 211                                         z=p[s2][sz*q+2]-ez;if(zperiodic) z+=pz;
 212                                         rr=x*x+y*y+z*z;
 213                                         if (rr<1) {
 214                                                 if(dd+rr+2*(dx*x+dy*y+dz*z)>=radsq) {
 215                                                         rr=alpha*(1-1/sqrt(rr));
 216                                                         ve[s1][3*l]+=x*rr;
 217                                                         ve[s1][3*l+1]+=y*rr;
 218                                                         ve[s1][3*l+2]+=z*rr;
 219                                                 }
 220                                         }
 221                                 }
 222                         }
 223                         wall_contribution(s1,l,ex,ey,ez,alpha);
 224                 }
 225         } while ((s1=l1.inc(ox,oy,oz))!=-1);
 226
 227         local_move(v_inter,cx,cy,cz,rad);
 228 }
 229
 230 template<class r_option>
 231 inline void container_dynamic_base<r_option>::wall_contribution(int s,int l,fpoint cx,fpoint cy,fpoint cz,fpoint alpha) {
 232         fpoint x,y,z,rr;
 233         for(int q=0;q<wall_number;q++) {
 234                 walls[q]->min_distance(cx,cy,cz,x,y,z);
 235                 rr=x*x+y*y+z*z;
 236                 if (rr<0.25) {
 237                         rr=0.5*alpha*(1-0.5/sqrt(rr));
 238                         ve[s][3*l]+=x*rr;
 239                         ve[s][3*l+1]+=y*rr;
 240                         ve[s][3*l+2]+=z*rr;
 241                 }
 242         }
 243 #ifdef VOROPP_AUTO_X_WALL
 244         if(!xperiodic) {
 245                 if(cx-ax<0.5) ve[s][3*l]+=alpha*(0.5+ax-cx);
 246                 if(bx-cx<0.5) ve[s][3*l]+=alpha*(bx-cx-0.5);
 247         }
 248 #endif
 249 #ifdef VOROPP_AUTO_Y_WALL
 250         if(!yperiodic) {
 251                 if(cy-ay<0.5) ve[s][3*l+1]+=alpha*(0.5+ay-cy);
 252                 if(by-cy<0.5) ve[s][3*l+1]+=alpha*(by-cy-0.5);
 253         }
 254 #endif
 255 #ifdef VOROPP_AUTO_Z_WALL
 256         if(!zperiodic) {
 257                 if(cz-az<0.5) ve[s][3*l+2]+=alpha*(0.5+az-cz);
 258                 if(bz-cz<0.5) ve[s][3*l+2]+=alpha*(bz-cz-0.5);
 259         }
 260 #endif
 261 }
 262
 263 template<class r_option>
 264 inline int container_dynamic_base<r_option>::full_count() {
 265         int ijk=0,count=0;
 266         while(ijk<nxyz) count+=co[ijk++];
 267         return count;
 268 }
 269
 270 template<class r_option>
 271 template<class cond_class>
 272 void container_dynamic_base<r_option>::neighbor_distribution(int *bb,fpoint dr,int max) {
 273         int i,j,k,ijk,l,s,q,ll;
 274         cond_class c1;
 275         voropp_loop l1(this);
 276         fpoint idr=1/dr,maxr=dr*max,maxrsq=maxr*maxr;
 277         fpoint x,y,z,px,py,pz,cx,cy,cz,rr;
 278         for(ijk=k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
 279                 for(l=0;l<co[ijk];l++) {
 280                         cx=p[ijk][sz*l];
 281                         cy=p[ijk][sz*l+1];
 282                         cz=p[ijk][sz*l+2];
 283                         if(c1.test(cx,cy,cz)) {
 284                                 s=l1.init(cx,cy,cz,maxr,px,py,pz);
 285                                 do {
 286                                         for(q=0;q<co[s];q++) {
 287                                                 x=p[s][sz*q]+px-cx;
 288                                                 y=p[s][sz*q+1]+py-cy;
 289                                                 z=p[s][sz*q+2]+pz-cz;
 290                                                 rr=x*x+y*y+z*z;
 291                                                 if(rr<maxrsq) {ll=int(idr*sqrt(rr));if(ll<max) bb[ll]++;}
 292                                         }
 293                                 } while ((s=l1.inc(px,py,pz))!=-1);
 294                         }
 295                 }
 296         }
 297 }
 298
 299
 300 template<class r_option>
 301 template<class cond_class>
 302 fpoint container_dynamic_base<r_option>::packing_badness() {
 303         int i,j,k,ijk,l,s,q;
 304         cond_class c1;
 305         voropp_loop l1(this);
 306         int pcount=0;
 307         fpoint x,y,z,px,py,pz,cx,cy,cz,rr,badcount=0;
 308         for(ijk=k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
 309                 for(l=0;l<co[ijk];l++) {
 310                         cx=p[ijk][sz*l];
 311                         cy=p[ijk][sz*l+1];
 312                         cz=p[ijk][sz*l+2];
 313                         if(c1.test(cx,cy,cz)) {
 314                                 pcount++;
 315                                 s=l1.init(cx,cy,cz,1,px,py,pz);
 316                                 do {
 317                                         for(q=0;q<co[s];q++) {
 318                                                 x=p[s][sz*q]+px-cx;
 319                                                 y=p[s][sz*q+1]+py-cy;
 320                                                 z=p[s][sz*q+2]+pz-cz;
 321                                                 rr=x*x+y*y+z*z;
 322                                                 if(rr<1&&rr>tolerance) badcount+=1-2*sqrt(rr)+rr;
 323                                         }
 324                                 } while ((s=l1.inc(px,py,pz))!=-1);
 325                                 wall_badness(cx,cy,cz,badcount);
 326                         }
 327                 }
 328         }
 329         return pcount==0?0:sqrt(badcount/pcount);
 330 }
 331
 332 template<class r_option>
 333 inline void container_dynamic_base<r_option>::wall_badness(fpoint cx,fpoint cy,fpoint cz,fpoint &badcount) {
 334         fpoint x,y,z,rr;
 335         for(int q=0;q<wall_number;q++) {
 336                 walls[q]->min_distance(cx,cy,cz,x,y,z);
 337                 rr=x*x+y*y+z*z;
 338                 if (rr<0.25) badcount+=0.25-sqrt(rr)+rr;
 339         }
 340 #ifdef VOROPP_AUTO_X_WALL
 341         if(!xperiodic) {
 342                 if(cx-ax<0.5) badcount+=(0.5+ax-cx)*(0.5+ax-cx);
 343                 if(bx-cx<0.5) badcount+=(bx-cx-0.5)*(bx-cx-0.5);
 344         }
 345 #endif
 346 #ifdef VOROPP_AUTO_Y_WALL
 347         if(!yperiodic) {
 348                 if(cy-ay<0.5) badcount+=(0.5+ay-cy)*(0.5+ay-cy);
 349                 if(by-cy<0.5) badcount+=(by-cy-0.5)*(by-cy-0.5);
 350         }
 351 #endif
 352 #ifdef VOROPP_AUTO_Z_WALL
 353         if(!zperiodic) {
 354                 if(cz-az<0.5) badcount+=(0.5+az-cz)*(0.5+az-cz);
 355                 if(bz-cz<0.5) badcount+=(bz-cz-0.5)*(bz-cz-0.5);
 356         }
 357 #endif
 358 }
 359
 360 #ifdef YEAST_ROUTINES
 361 /** An overloaded version of the draw_yeast_pov() routine, that outputs
 362  * the particle positions to a file.
 363  * \param[in] filename the file to write to. */
 364 template<class r_option>
 365 void container_dynamic_base<r_option>::draw_yeast_pov(const char *filename) {
 366         ofstream os;
 367         os.open(filename,ofstream::out|ofstream::trunc);
 368         draw_yeast_pov(os);
 369         os.close();
 370 }
 371
 372 /** Dumps all the particle positions in the POV-Ray format. */
 373 template<class r_option>
 374 void container_dynamic_base<r_option>::draw_yeast_pov(ostream &os) {
 375         int l,c;
 376         for(l=0;l<nxyz;l++) {
 377                 for(c=0;c<co[l];c++) {
 378                         os << "// id " << id[l][c] << "\n";
 379                         os << "sphere{<" << p[l][sz*c] << "," << p[l][sz*c+1] << ","
 380                            << p[l][sz*c+2] << ">,";
 381                         radius.rad(os,l,c);
 382                         if(id[l][c]<stickycut) {
 383                                 os << " pigment{rgb <0.9,0.3,0.6>} finish{reflection 0.15 ambient 0.4 specular 0.3}}\n";
 384                         } else {
 385                                 os << " pigment{rgb <0.9,0.85,0.35>} finish{reflection 0.15 specular 0.3 ambient 0.42}}\n";
 386                         }
 387                 }
 388         }
 389 }
 390
 391 template<class r_option>
 392 void container_dynamic_base<r_option>::stick(fpoint alpha) {
 393         int i,j,k,ijk,l,q,s;
 394         fpoint bigrad,smallrad,currad;
 395         fpoint bigfac,smallfac,curfac;
 396         fpoint x,y,z,px,py,pz,cx,cy,cz,rr;
 397         voropp_loop l1(this);
 398
 399         for(ijk=0;ijk<nxyz;ijk++) for(l=0;l<3*co[ijk];l++) ve[ijk][l]=0;
 400
 401         for(ijk=k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
 402                 for(l=0;l<co[ijk];l++) {
 403                         cx=p[ijk][sz*l];
 404                         cy=p[ijk][sz*l+1];
 405                         cz=p[ijk][sz*l+2];
 406                         if(id[ijk][l]<stickycut) {
 407                                 bigrad=1.6;smallrad=1.2;
 408                                 bigfac=1/0.8;smallfac=1/0.4;
 409                         } else {
 410                                 bigrad=1.2;smallrad=0.799999999;
 411                                 bigfac=1/0.4;smallfac=0;
 412                         }
 413                         s=l1.init(cx,cy,cz,bigrad,px,py,pz);
 414                         while(!l1.reached(i,j,k)) {
 415                                 for(q=0;q<co[s];q++) {
 416                                         x=p[s][sz*q]+px-cx;
 417                                         y=p[s][sz*q+1]+py-cy;
 418                                         z=p[s][sz*q+2]+pz-cz;
 419                                         rr=x*x+y*y+z*z;
 420                                         if(id[s][q]<stickycut) {
 421                                                 currad=bigrad;curfac=bigfac;
 422                                         } else {
 423                                                 currad=smallrad;curfac=smallfac;
 424                                         }
 425                                         if (rr<currad*currad) {
 426                                                 if(rr<0.8*0.8) rr=0.5*alpha*(1-0.8/sqrt(rr));
 427                                                 else {rr=sqrt(rr);rr=0.5*alpha*(rr-0.8)*(1-(rr-0.8)*curfac*(2-(rr-0.8)*curfac))/rr;}
 428                                                 ve[ijk][3*l]+=x*rr;
 429                                                 ve[ijk][3*l+1]+=y*rr;
 430                                                 ve[ijk][3*l+2]+=z*rr;
 431                                                 ve[s][3*q]-=x*rr;
 432                                                 ve[s][3*q+1]-=y*rr;
 433                                                 ve[s][3*q+2]-=z*rr;
 434                                         }
 435                                 }
 436                                 s=l1.inc(px,py,pz);
 437                         }
 438                         if (s!=ijk) throw fatal_error("s and ijk should be aligned and they're not");
 439                         for(q=0;q<l;q++) {
 440                                 x=p[s][sz*q]+px-cx;
 441                                 y=p[s][sz*q+1]+py-cy;
 442                                 z=p[s][sz*q+2]+pz-cz;
 443                                 rr=x*x+y*y+z*z;
 444                                 if(id[s][q]<stickycut) {
 445                                         currad=bigrad;curfac=bigfac;
 446                                 } else {
 447                                         currad=smallrad;curfac=smallfac;
 448                                 }
 449                                 if (rr<currad*currad) {
 450                                         if(rr<0.8*0.8) rr=0.5*alpha*(1-0.8/sqrt(rr));
 451                                         else {rr=sqrt(rr);rr=0.5*alpha*(rr-0.8)*(1-(rr-0.8)*curfac*(2-(rr-0.8)*curfac))/rr;}
 452                                         ve[ijk][3*l]+=x*rr;
 453                                         ve[ijk][3*l+1]+=y*rr;
 454                                         ve[ijk][3*l+2]+=z*rr;
 455                                         ve[s][3*q]-=x*rr;
 456                                         ve[s][3*q+1]-=y*rr;
 457                                         ve[s][3*q+2]-=z*rr;
 458                                 }
 459                         }
 460                         wall_contribution(ijk,l,cx,cy,cz,alpha);
 461                 }
 462         }
 463
 464         move(v_inter);
 465 }
 466 #endif
 467
 468 template<class r_option>
 469 inline void container_dynamic_base<r_option>::clear_velocities() {
 470         for(int l,ijk=0;ijk<nxyz;ijk++) for(l=0;l<3*co[ijk];l++) ve[ijk][l]=0;
 471 }
 472
 473 template<class r_option>
 474 inline void container_dynamic_base<r_option>::damp_velocities(fpoint damp) {
 475         for(int l,ijk=0;ijk<nxyz;ijk++) for(l=0;l<3*co[ijk];l++) ve[ijk][l]*=damp;
 476 }
 477
 478 template<class r_option>
 479 void container_dynamic_base<r_option>::full_relax(fpoint alpha) {
 480         int i,j,k,ijk,l,q,s;
 481         fpoint x,y,z,px,py,pz,cx,cy,cz,rr;
 482         voropp_loop l1(this);
 483
 484         clear_velocities();
 485
 486         for(ijk=k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
 487                 for(l=0;l<co[ijk];l++) {
 488                         cx=p[ijk][sz*l];
 489                         cy=p[ijk][sz*l+1];
 490                         cz=p[ijk][sz*l+2];
 491                         s=l1.init(cx,cy,cz,1,px,py,pz);
 492                         while(!l1.reached(i,j,k)) {
 493                                 for(q=0;q<co[s];q++) {
 494                                         x=p[s][sz*q]+px-cx;
 495                                         y=p[s][sz*q+1]+py-cy;
 496                                         z=p[s][sz*q+2]+pz-cz;
 497                                         rr=x*x+y*y+z*z;
 498                                         if (rr<1) {
 499                                                 rr=alpha*(0.5-0.5/sqrt(rr));
 500                                                 ve[ijk][3*l]+=x*rr;
 501                                                 ve[ijk][3*l+1]+=y*rr;
 502                                                 ve[ijk][3*l+2]+=z*rr;
 503                                                 ve[s][3*q]-=x*rr;
 504                                                 ve[s][3*q+1]-=y*rr;
 505                                                 ve[s][3*q+2]-=z*rr;
 506                                         }
 507                                 }
 508                                 s=l1.inc(px,py,pz);
 509                         }
 510                         if (s!=ijk) throw fatal_error("s and ijk should be aligned and they're not");
 511                         for(q=0;q<l;q++) {
 512                                 x=p[s][sz*q]+px-cx;
 513                                 y=p[s][sz*q+1]+py-cy;
 514                                 z=p[s][sz*q+2]+pz-cz;
 515                                 rr=x*x+y*y+z*z;
 516                                 if (rr<1) {
 517                                         rr=alpha*(0.5-0.5/sqrt(rr));
 518                                         ve[ijk][3*l]+=x*rr;
 519                                         ve[ijk][3*l+1]+=y*rr;
 520                                         ve[ijk][3*l+2]+=z*rr;
 521                                         ve[s][3*q]-=x*rr;
 522                                         ve[s][3*q+1]-=y*rr;
 523                                         ve[s][3*q+2]-=z*rr;
 524                                 }
 525                         }
 526                         wall_contribution(ijk,l,cx,cy,cz,alpha);
 527                 }
 528         }
 529
 530         move(v_inter);
 531 }
 532
 533 /** Increase memory for a particular region. */
 534 template<class r_option>
 535 void container_dynamic_base<r_option>::add_particle_memory(int i) {
 536         int *idp;fpoint *pp,*vep;
 537         int l,nmem=2*mem[i];
 538 #if VOROPP_VERBOSE >=3
 539         cerr << "Particle memory in region " << i << " scaled up to " << nmem << endl;
 540 #endif
 541         if(nmem>max_particle_memory) throw fatal_error("Absolute maximum memory allocation exceeded");
 542         idp=new int[nmem];
 543         for(l=0;l<co[i];l++) idp[l]=id[i][l];
 544         pp=new fpoint[sz*nmem];
 545         for(l=0;l<sz*co[i];l++) pp[l]=p[i][l];
 546         vep=new fpoint[3*nmem];
 547         for(l=0;l<3*co[i];l++) vep[l]=ve[i][l];
 548         mem[i]=nmem;
 549         delete [] id[i];id[i]=idp;
 550         delete [] p[i];p[i]=pp;
 551         delete [] ve[i];ve[i]=vep;
 552 }
 553
 554 /** Custom int function, that gives consistent stepping for negative numbers.
 555  * With normal int, we have (-1.5,-0.5,0.5,1.5) -> (-1,0,0,1).
 556  * With this routine, we have (-1.5,-0.5,0.5,1.5) -> (-2,-1,0,1).*/
 557 template<class r_option>
 558 inline int container_dynamic_base<r_option>::step_int(fpoint a) {
 559         return a<0?int(a)-1:int(a);
 560 }
 561
 562 /** Custom mod function, that gives consistent stepping for negative numbers. */
 563 template<class r_option>
 564 inline int container_dynamic_base<r_option>::step_mod(int a,int b) {
 565         return a>=0?a%b:b-1-(b-1-a)%b;
 566 }
 567
 568 /** Custom div function, that gives consistent stepping for negative numbers. */
 569 template<class r_option>
 570 inline int container_dynamic_base<r_option>::step_div(int a,int b) {
 571         return a>=0?a/b:-1+(a+1)/b;
 572 }
 573
 574 template<class r_option>
 575 template<class v_class>
 576 void container_dynamic_base<r_option>::move(v_class &vcl) {
 577         int i,j,k,ijk,l,ll,ni,nj,nk;
 578         fpoint x,y,z;
 579
 580         // For each block, introduce a second counter that counts the number of
 581         // particles initially in this block, without including any that have
 582         // moved into the block
 583         for(ijk=0;ijk<nxyz;ijk++) gh[ijk]=co[ijk];
 584
 585         // Loop over the blocks within the counter
 586         for(ijk=k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
 587                 l=0;
 588
 589                 // Scan over the particles that started within this block
 590                 while(l<gh[ijk]) {
 591
 592                         // Get the particle's initial position
 593                         x=p[ijk][sz*l];y=p[ijk][sz*l+1];z=p[ijk][sz*l+2];
 594
 595                         // Call an external routine to compute where this
 596                         // particle should move to
 597                         vcl.vel(ijk,l,x,y,z);
 598
 599                         // Calculate which block the particle's new position is
 600                         // within
 601                         ni=step_int((x-ax)*xsp);
 602                         nj=step_int((y-ay)*ysp);
 603                         nk=step_int((z-az)*zsp);
 604
 605                         // See if the particle is within the same block
 606                         if(ni==i&&nj==j&&nk==k) {
 607
 608                                 // If so, just update its position
 609                                 p[ijk][sz*l]=x;
 610                                 p[ijk][sz*l+1]=y;
 611                                 p[ijk][sz*l+2]=z;l++;
 612                         } else {
 613
 614                                 // Check whether the particle's new position is
 615                                 // within the container
 616                                 if((xperiodic||(ni>=0&&ni<nx))&&(yperiodic||(nj>=0&&nj<ny))&&(zperiodic||(nk>=0&&nk<nz))) {
 617
 618                                         // If periodic boundary conditions are
 619                                         // being used, then remap the particle
 620                                         // back into the container
 621                                         if(xperiodic) {x-=step_div(ni,nx)*(bx-ax);ni=step_mod(ni,nx);}
 622                                         if(yperiodic) {y-=step_div(nj,ny)*(by-ay);nj=step_mod(nj,ny);}
 623                                         if(zperiodic) {z-=step_div(nk,nz)*(bz-az);nk=step_mod(nk,nz);}
 624
 625                                         // Calculate the index the new block
 626                                         // where the particle is to be stored,
 627                                         // and allocate memory if necessary
 628                                         ni+=nx*(nj+ny*nk);
 629                                         if(co[ni]==mem[ni]) add_particle_memory(ni);
 630
 631                                         // Add the particle to the new block
 632                                         id[ni][co[ni]]=id[ijk][l];
 633                                         p[ni][co[ni]*sz]=x;
 634                                         p[ni][co[ni]*sz+1]=y;
 635                                         p[ni][co[ni]*sz+2]=z;
 636
 637                                         // This class can support a velocity
 638                                         // array. Move it too, if it is being
 639                                         // used.
 640                                         if(vcl.track_ve) {
 641                                                 ve[ni][3*co[ni]]=ve[ijk][3*l];
 642                                                 ve[ni][3*co[ni]+1]=ve[ijk][3*l+1];
 643                                                 ve[ni][3*co[ni]+2]=ve[ijk][3*l+2];
 644                                         }
 645
 646                                         // Copy particle radius information if
 647                                         // it being used, and add one to the
 648                                         // number of particles in the block
 649                                         if(sz==4) p[ni][co[ni]*sz+3]=p[ijk][l*sz+3];
 650                                         co[ni]++;
 651                                 }
 652
 653                                 // Delete the particle from the original block,
 654                                 // by copying the final particle in the block's
 655                                 // memory on top of it
 656                                 co[ijk]--;
 657                                 id[ijk][l]=id[ijk][co[ijk]];
 658                                 for(ll=0;ll<sz;ll++) p[ijk][sz*l+ll]=p[ijk][sz*co[ijk]+ll];
 659                                 if (vcl.track_ve) {
 660                                         ve[ijk][3*l]=ve[ijk][3*co[ijk]];
 661                                         ve[ijk][3*l+1]=ve[ijk][3*co[ijk]+1];
 662                                         ve[ijk][3*l+2]=ve[ijk][3*co[ijk]+2];
 663                                 }
 664
 665                                 // If the copied particle was one that had
 666                                 // initially started off in this block, then
 667                                 // prepare to check that one next. If it was a
 668                                 // particle that had moved into this block,
 669                                 // then skip it.
 670                                 if(co[ijk]+1==gh[ijk]) gh[ijk]--;
 671                                 else l++;
 672                         }
 673                 }
 674         }
 675 }