trunk/examples/no_release/ellipsoid.cc

   1 // Ellipsoidal wall example
   2 //
   3 // Author   : Simon Konstandin (UMM Mannheim)
   4
   5 #include "voro++.hh"
   6 using namespace voro;
   7
   8 #include "quartic/quartic.hpp"
   9 using namespace magnet::math;
  10
  11 // Set up constants for the container geometry
  12 const double x_min=-6,x_max=6;
  13 const double y_min=-6,y_max=6;
  14 const double z_min=-6,z_max=6;
  15
  16 // Golden ratio constants
  17 const double Phi=0.5*(1+sqrt(5.0));
  18 const double phi=0.5*(1-sqrt(5.0));
  19
  20 // Set up the number of blocks that the container is divided
  21 // into.
  22 const int n_x=5,n_y=5,n_z=5;
  23
  24 // Set the number of particles that are going to be randomly introduced
  25 const int particles=500;
  26
  27 // This function returns a random double between 0 and 1
  28 double rnd() {return double(rand())/RAND_MAX;}
  29
  30 // Ellipsoidal wall object
  31 class wall_ellipsoid : public wall {
  32         public:
  33                 const double a;
  34                 const double aa;
  35                 const double aa_inv;
  36                 const double a_cal;
  37                 const double rr;
  38                 const double rr_inv;
  39                 const double xc;
  40                 const double yc;
  41                 const double zc;
  42                 const int w_id;
  43                 wall_ellipsoid(double a_,double r_,double xc_,double yc_,double zc_)
  44                         : a(a_), aa(a_*a_), aa_inv(1/aa), a_cal(2.0+2.0/aa), rr(r_*r_),
  45                         rr_inv(1/rr), xc(xc_), yc(yc_), zc(zc_), w_id(99) {}
  46                 bool point_inside(double x,double y,double z) {
  47                         return x*x+y*y+z*z*a*a<rr;
  48                 }
  49                 virtual bool cut_cell(voronoicell &c,double x,double y,double z) {
  50                         return cut_cell_base(c,x,y,z);
  51                 }
  52                 virtual bool cut_cell(voronoicell_neighbor &c,double x,double y,double z) {
  53                         return cut_cell_base(c,x,y,z);
  54                 }
  55                 template<class v_cell>
  56                 bool cut_cell_base(v_cell &c,double x,double y,double z) {
  57                         x-=xc;y-=yc;z-=zc;
  58                         double xx=x*x,yy=y*y,zz=z*z,s[4],
  59                                 b_cal=1.0-(xx+yy)*rr_inv+(4.0-zz*rr_inv)*aa_inv+aa_inv*aa_inv,
  60                                 c_cal=2.0*(1.0-(xx+yy+zz)*rr_inv)*aa_inv+2.0*aa_inv*aa_inv,
  61                                 d_cal=(1.0-(xx+yy+zz*aa)*rr_inv)*aa_inv*aa_inv,
  62                                 xell[4],yell[4],zell[4],d2[4];
  63                         int i,j=0,q=quarticSolve(a_cal,b_cal,c_cal,d_cal,s[0],s[1],s[2],s[3]);
  64
  65                         for (i=0;i<q;i++) {
  66                                 xell[i]=x/(1+s[i]);
  67                                 yell[i]=y/(1+s[i]);
  68                                 zell[i]=z/(1+s[i]*aa);
  69                                 d2[i]=(xell[i]-x)*(xell[i]-x)+(yell[i]-y)*(yell[i]-y)+(zell[i]-z)*(zell[i]-z);
  70                         }
  71
  72                         for(i=1;i<4;i++) if ((d2[i]<=d2[j])&&(abs(xell[i]*xell[i]+yell[i]*yell[i]+aa*zell[i]*zell[i]-rr)<1e-4)) j=i;
  73
  74                         double xn=xell[j]-x,yn=yell[j]-y,zn=zell[j]-z,dn2=xn*xn+yn*yn+zn*zn;
  75
  76                         return c.nplane(xn,yn,zn,xx+yy+aa*zz<rr?2*dn2:0,w_id);
  77                 }
  78 };
  79
  80 int main() {
  81         int i;
  82         double x,y,z;
  83
  84         // Create a container with the geometry given above. This is bigger
  85         // than the particle packing itself.
  86         container con(x_min,x_max,y_min,y_max,z_min,z_max,n_x,n_y,n_z,
  87                         false,false,false,8);
  88
  89         // Create the "initial shape" wall class and add it to the container
  90         wall_ellipsoid we(1.5,5,0,0,0);
  91         con.add_wall(we);
  92
  93         // Randomly add particles into the container
  94         for(i=0;i<particles;) {
  95                 x=x_min+rnd()*(x_max-x_min);
  96                 y=y_min+rnd()*(y_max-y_min);
  97                 z=z_min+rnd()*(z_max-z_min);
  98                 if(con.point_inside(x,y,z)) con.put(i++,x,y,z);
  99         }
 100
 101         // Save the particles and Voronoi cells in POV-Ray format
 102         con.draw_particles("ellipsoid_p.gnu");
 103         con.draw_cells_gnuplot("ellipsoid_v.gnu");
 104 }