branches/2d/src/cell_nc_2d.cc

   1 // Voro++, a cell-based Voronoi library
   2 //
   3 // Authors  : Chris H. Rycroft (LBL / UC Berkeley)
   4 //            Cody Robert Dance (UC Berkeley)
   5 // Email    : chr@alum.mit.edu
   6 // Date     : August 30th 2011
   7
   8 /** \file cell_nc_2d.cc
   9  * \brief Function implementations for the non-convex 2D Voronoi classes. */
  10
  11 #include "cell_nc_2d.hh"
  12
  13 namespace voro {
  14
  15 void voronoicell_nonconvex_neighbor_2d::init(double xmin,double xmax,double ymin,double ymax) {
  16         nonconvex=exclude=false;
  17         init_base(xmin,xmax,ymin,ymax);
  18         *ne=-3;ne[1]=-2;ne[2]=-4;ne[3]=-1;
  19 }
  20
  21 void voronoicell_nonconvex_base_2d::init_nonconvex_base(double xmin,double xmax,double ymin,double ymax,double wx0,double wy0,double wx1,double wy1) {
  22         xmin*=2;xmax*=2;ymin*=2;ymax*=2;
  23         int f0=face(xmin,xmax,ymin,ymax,wx0,wy0),
  24             f1=face(xmin,xmax,ymin,ymax,wx1,wy1);
  25
  26         *pts=0;pts[1]=0;
  27         pts[2]=wx0;pts[3]=wy0;p=4;
  28         if(f0!=f1||wx0*wy1<wx1*wy0) {
  29                 do {
  30                         if(f0>1) {
  31                                 if(f0==2) {pts[p++]=xmin;pts[p++]=ymin;}
  32                                 else {pts[p++]=xmax;pts[p++]=ymin;}
  33                         } else {
  34                                 if(f0==0) {pts[p++]=xmax;pts[p++]=ymax;}
  35                                 else {pts[p++]=xmin;pts[p++]=ymax;}
  36                         }
  37                         f0++;f0&=3;
  38                 } while(f0!=f1);
  39         }
  40         pts[p++]=wx1;pts[p++]=wy1;
  41         p>>=1;
  42
  43         int i,*q=ed;
  44         *(q++)=1;*(q++)=p-1;
  45         for(i=1;i<p-1;i++) {*(q++)=i+1;*(q++)=i-1;}
  46         *(q++)=0;*q=p-2;
  47
  48         exclude=true;
  49         if(wx0*wy1>wx1*wy0) nonconvex=false;
  50         else {
  51                 nonconvex=true;
  52                 if(wx0*wx1+wy0*wy1>0) {
  53                         *reg=wy0+wy1;reg[1]=-wx0-wx1;
  54                 } else {
  55                         *reg=wx1-wx0;reg[1]=wy1-wy0;
  56                 }
  57         }
  58         reg[2]=wx0;reg[3]=wy0;
  59         reg[4]=wx1;reg[5]=wy1;
  60 }
  61
  62 void voronoicell_nonconvex_neighbor_2d::init_nonconvex(double xmin,double xmax,double ymin,double ymax,double wx0,double wy0,double wx1,double wy1) {
  63         init_nonconvex_base(xmin,xmax,ymin,ymax,wx0,wy0,wx1,wy1);
  64         *ne=-5;
  65         for(int i=1;i<p-1;i++) ne[i]=-99;
  66         ne[p-1]=-5;
  67 }
  68
  69 /** Cuts the Voronoi cell by a particle whose center is at a separation of
  70  * (x,y) from the cell center. The value of rsq should be initially set to
  71  * \f$x^2+y^2\f$.
  72  * \param[in] (x,y) the normal vector to the plane.
  73  * \param[in] rsq the distance along this vector of the plane.
  74  * \return False if the plane cut deleted the cell entirely, true otherwise. */
  75 template<class vc_class>
  76 bool voronoicell_nonconvex_base_2d::nplane_nonconvex(vc_class &vc,double x,double y,double rsq,int p_id) {
  77         int up=0,*edd;
  78         double u,rx,ry,sx,sy;
  79
  80         if(x*(*reg)+y*reg[1]<0) {edd=ed;rx=reg[2];ry=reg[3];sx=*reg;sy=reg[1];}
  81         else {edd=ed+1;rx=reg[4];ry=reg[5];sx=-*reg;sy=-reg[1];}
  82
  83
  84         up=*edd;u=pos(x,y,rsq,up);
  85         while(u<tolerance) {
  86                 up=edd[2*up];if(up==0) return true;
  87                 if(pts[2*up]*sx+pts[2*up+1]*sy>0&&rx*pts[2*up]+ry*pts[2*up+1]>0) return true;
  88                 u=pos(x,y,rsq,up);
  89         }
  90
  91         return nplane_cut(vc,x,y,rsq,p_id,u,up);
  92 }
  93
  94
  95 inline int voronoicell_nonconvex_base_2d::face(double xmin,double xmax,double ymin,double ymax,double &wx,double &wy) {
  96         if(wy>0) {
  97                 if(xmin*wy>ymax*wx) {wy*=xmin/wx;wx=xmin;return 2;}
  98                 if(xmax*wy>ymax*wx) {wx*=ymax/wy;wy=ymax;return 1;}
  99                 wy*=xmax/wx;wx=xmax;return 0;
 100         }
 101         if(xmax*wy>ymin*wx) {wy*=xmax/wx;wx=xmax;return 0;}
 102         if(xmin*wy>ymin*wx) {wx*=ymin/wy;wy=ymin;return 3;}
 103         wy*=xmin/wx;wx=xmin;return 2;
 104 }
 105
 106 void voronoicell_nonconvex_neighbor_2d::neighbors(vector<int> &v) {
 107         v.resize(p);
 108         for(int i=0;i<p;i++) v[i]=ne[i];
 109 }
 110
 111 // Explicit instantiation
 112 template bool voronoicell_nonconvex_base_2d::nplane_nonconvex(voronoicell_nonconvex_2d&,double,double,double,int);
 113 template bool voronoicell_nonconvex_base_2d::nplane_nonconvex(voronoicell_nonconvex_neighbor_2d&,double,double,double,int);
 114
 115 }