branches/2d/src/v_connect/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=full_connect=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         full_connect=false;
  23         xmin*=2;xmax*=2;ymin*=2;ymax*=2;
  24         int f0=face(xmin,xmax,ymin,ymax,wx0,wy0),
  25             f1=face(xmin,xmax,ymin,ymax,wx1,wy1);
  26
  27         *pts=0;pts[1]=0;
  28         pts[2]=wx0;pts[3]=wy0;p=4;
  29         if(f0!=f1||wx0*wy1<wx1*wy0) {
  30                 do {//VTG INFO???
  31                         if(f0>1) {
  32                                 if(f0==2) {pts[p++]=xmin;pts[p++]=ymin;}
  33                                 else {pts[p++]=xmax;pts[p++]=ymin;}
  34                         } else {
  35                                 if(f0==0) {pts[p++]=xmax;pts[p++]=ymax;}
  36                                 else {pts[p++]=xmin;pts[p++]=ymax;}
  37                         }
  38                         f0++;f0&=3;
  39                 } while(f0!=f1);
  40         }
  41         pts[p++]=wx1;pts[p++]=wy1;
  42         p>>=1;
  43
  44         int i,*q=ed;
  45         *(q++)=1;*(q++)=p-1;
  46         for(i=1;i<p-1;i++) {*(q++)=i+1;*(q++)=i-1;}
  47         *(q++)=0;*q=p-2;
  48
  49
  50
  51
  52
  53
  54
  55
  56         exclude=true;
  57         if(wx0*wy1>wx1*wy0) nonconvex=false;
  58         else {
  59                 nonconvex=true;
  60                 if(wx0*wx1+wy0*wy1>0) {
  61                         *reg=wy0+wy1;reg[1]=-wx0-wx1;
  62                 } else {
  63                         *reg=wx1-wx0;reg[1]=wy1-wy0;
  64                 }
  65         }
  66         reg[2]=wx0;reg[3]=wy0;
  67         reg[4]=wx1;reg[5]=wy1;
  68 }
  69
  70 void voronoicell_nonconvex_neighbor_2d::init_nonconvex(double xmin,double xmax,double ymin,double ymax,double wx0,double wy0,double wx1,double wy1) {
  71         init_nonconvex_base(xmin,xmax,ymin,ymax,wx0,wy0,wx1,wy1);
  72         *ne=-5;
  73         for(int i=1;i<p-1;i++) ne[i]=-99;
  74         ne[p-1]=-5;
  75 }
  76
  77 /** Cuts the Voronoi cell by a particle whose center is at a separation of
  78  * (x,y) from the cell center. The value of rsq should be initially set to
  79  * \f$x^2+y^2\f$.
  80  * \param[in] (x,y) the normal vector to the plane.
  81  * \param[in] rsq the distance along this vector of the plane.
  82  * \return False if the plane cut deleted the cell entirely, true otherwise. */
  83 template<class vc_class>
  84 bool voronoicell_nonconvex_base_2d::nplane_nonconvex(vc_class &vc,double x,double y,double rsq,int p_id) {
  85         int up=0,*edd;
  86         double u,rx,ry,sx,sy;
  87
  88         if(x*(*reg)+y*reg[1]<0) {edd=ed;rx=reg[2];ry=reg[3];sx=*reg;sy=reg[1];}
  89         else {edd=ed+1;rx=reg[4];ry=reg[5];sx=-*reg;sy=-reg[1];}
  90
  91
  92         up=*edd;u=pos(x,y,rsq,up);
  93         while(u<tolerance) {
  94                 up=edd[2*up];if(up==0) return true;
  95                 if(pts[2*up]*sx+pts[2*up+1]*sy>0&&rx*pts[2*up]+ry*pts[2*up+1]>0) return true;
  96                 u=pos(x,y,rsq,up);
  97         }
  98
  99         return nplane_cut(vc,x,y,rsq,p_id,u,up);
 100 }
 101
 102
 103 inline int voronoicell_nonconvex_base_2d::face(double xmin,double xmax,double ymin,double ymax,double &wx,double &wy) {
 104         if(wy>0) {
 105                 if(xmin*wy>ymax*wx) {wy*=xmin/wx;wx=xmin;return 2;}
 106                 if(xmax*wy>ymax*wx) {wx*=ymax/wy;wy=ymax;return 1;}
 107                 wy*=xmax/wx;wx=xmax;return 0;
 108         }
 109         if(xmax*wy>ymin*wx) {wy*=xmax/wx;wx=xmax;return 0;}
 110         if(xmin*wy>ymin*wx) {wx*=ymin/wy;wy=ymin;return 3;}
 111         wy*=xmin/wx;wx=xmin;return 2;
 112 }
 113
 114 void voronoicell_nonconvex_neighbor_2d::neighbors(vector<int> &v) {
 115         v.resize(p);
 116         for(int i=0;i<p;i++) v[i]=ne[i];
 117 }
 118
 119 // Explicit instantiation
 120 template bool voronoicell_nonconvex_base_2d::nplane_nonconvex(voronoicell_nonconvex_2d&,double,double,double,int);
 121 template bool voronoicell_nonconvex_base_2d::nplane_nonconvex(voronoicell_nonconvex_neighbor_2d&,double,double,double,int);
 122
 123 }