1#!/usr/bin/perl 2# Voro++, a 3D cell-based Voronoi library 3# 4# Author : Chris H. Rycroft (LBL / UC Berkeley) 5# Email : chr@alum.mit.edu 6# Date : August 30th 2011 7# 8# worklist_gen.pl - this perl script is used to automatically generate the 9# worklists of blocks that are stored in worklist.cc, that are used by the 10# compute_cell() routine to ensure maximum efficiency 11 12# Each region is divided into a grid of subregions, and a worklist is 13# constructed for each. This parameter sets is set to half the number of 14# subregions that the block is divided into. 15$hr=4; 16 17# This parameter is automatically set to the the number of subregions that the 18# block is divided into 19$r=$hr*2; 20 21# This parameter sets the total number of block entries in each worklist 22$ls=63; 23 24# When the worklists are being constructed, a mask array is made use of. To 25# prevent the creation of array elements with negative indices, this parameter 26# sets an arbitrary displacement. 27$dis=8; 28 29# Constants used mask indexing 30$d=2*$dis+1;$dd=$d*$d;$d0=(1+$d+$dd)*$dis; 31 32use Math::Trig; 33 34# Construct the worklist header file, based on the parameters above 35open W,">worklist.hh"; 36print W <<EOF; 37// Voro++, a 3D cell-based Voronoi library 38// 39// Author : Chris H. Rycroft (LBL / UC Berkeley) 40// Email : chr\@alum.mit.edu 41// Date : August 30th 2011 42 43/** \\file worklist.hh 44 * \\brief Header file for setting constants used in the block worklists that are 45 * used during cell computation. 46 * 47 * This file is automatically generated by worklist_gen.pl and it is not 48 * intended to be edited by hand. */ 49 50#ifndef VOROPP_WORKLIST_HH 51#define VOROPP_WORKLIST_HH 52 53namespace voro { 54 55/** Each region is divided into a grid of subregions, and a worklist is 56# constructed for each. This parameter sets is set to half the number of 57# subregions that the block is divided into. */ 58EOF 59print W "const int wl_hgrid=$hr;\n"; 60print W <<EOF; 61/** The number of subregions that a block is subdivided into, which is twice 62the value of hgrid. */ 63EOF 64print W "const int wl_fgrid=$r;\n"; 65print W <<EOF; 66/** The total number of worklists, set to the cube of hgrid. */ 67EOF 68printf W "const int wl_hgridcu=%d;\n",$hr*$hr*$hr; 69print W <<EOF; 70/** The number of elements in each worklist. */ 71EOF 72printf W "const int wl_seq_length=%d;\n",$ls+1; 73print W <<EOF; 74 75} 76#endif 77EOF 78close W; 79 80# Construct the preamble to the worklist.cc file 81open W,">v_base_wl.cc"; 82print W <<EOF; 83// Voro++, a 3D cell-based Voronoi library 84// 85// Author : Chris H. Rycroft (LBL / UC Berkeley) 86// Email : chr\@alum.mit.edu 87// Date : August 30th 2011 88 89/** \\file v_base_wl.cc 90 * \\brief The table of block worklists that are used during the cell 91 * computation, which is part of the voro_base class. 92 * 93 * This file is automatically generated by worklist_gen.pl and it is not 94 * intended to be edited by hand. */ 95 96EOF 97printf W "const unsigned int voro_base::wl[wl_seq_length*wl_hgridcu]={\n"; 98 99# Now create a worklist for each subregion 100for($kk=0;$kk<$hr;$kk++) { 101 for($jj=0;$jj<$hr;$jj++) { 102 for($ii=0;$ii<$hr;$ii++) { 103 worklist($ii,$jj,$kk); 104 } 105 } 106} 107 108# Finish the file and close it 109print W "};\n"; 110close W; 111 112# A subroutine 113sub worklist { 114# print "@_[0] @_[1] @_[2]\n"; 115 $ind=@_[0]+$hr*(@_[1]+$hr*@_[2]); 116 $ac=0;$v+=2; 117 $xp=$yp=$zp=0; 118 $x=(@_[0]+0.5)/$r; 119 $y=(@_[1]+0.5)/$r; 120 $z=(@_[2]+0.5)/$r; 121 $m[$d0]=$v; 122 add(1,0,0);add(0,1,0);add(0,0,1); 123 add(-1,0,0);add(0,-1,0);add(0,0,-1); 124 foreach $l (1..$ls) { 125 $minwei=1e9; 126 foreach (0..$ac-1) { 127 $xt=@a[3*$_];$yt=@a[3*$_+1];$zt=@a[3*$_+2]; 128# $wei=dis($x,$y,$z,$xt,$yt,$zt)+1*acos(($xt*$xp+$yt*$yp+$zt*$zp)/($xt*$xt+$yt*$yt+$zt*$zt)*($xp*$xp+$yp*$yp+$zp*$zp)); 129 $wei=adis($x,$y,$z,$xt,$yt,$zt)+0.02*sqrt(($xt-$xp)**2+($yt-$yp)**2+($zt-$zp)**2); 130 $nx=$_,$minwei=$wei if $wei<$minwei; 131 } 132 $xp=@a[3*$nx];$yp=@a[3*$nx+1];$zp=@a[3*$nx+2]; 133 add($xp+1,$yp,$zp);add($xp,$yp+1,$zp);add($xp,$yp,$zp+1); 134 add($xp-1,$yp,$zp);add($xp,$yp-1,$zp);add($xp,$yp,$zp-1); 135# print "=> $l $xp $yp $zp\n" if $l<4; 136 push @b,(splice @a,3*$nx,3);$ac--; 137 } 138 139 # Mark all blocks that are on this worklist entry 140 $m[$d0]=++$v; 141 for($i=0;$i<$#b;$i+=3) { 142 $xt=$b[$i];$yt=$b[$i+1];$zt=$b[$i+2]; 143 $m[$d0+$xt+$d*$yt+$dd*$zt]=$v; 144 } 145 146 # Find which neighboring outside blocks need to be marked when 147 # considering this block, being as conservative as possible by 148 # overwriting the marks, so that the last possible entry that can reach 149 # a block is used 150 for($i=$j=0;$i<$#b;$i+=3,$j++) { 151 $xt=$b[$i];$yt=$b[$i+1];$zt=$b[$i+2]; 152 $k=$d0+$xt+$yt*$d+$zt*$dd; 153 $la[$k+1]=$j, $m[$k+1]=$v+1 if $xt>=0 && $m[$k+1]!=$v; 154 $la[$k+$d]=$j, $m[$k+$d]=$v+1 if $yt>=0 && $m[$k+$d]!=$v; 155 $la[$k+$dd]=$j, $m[$k+$dd]=$v+1 if $zt>=0 && $m[$k+$dd]!=$v; 156 $la[$k-1]=$j, $m[$k-1]=$v+1 if $xt<=0 && $m[$k-1]!=$v; 157 $la[$k-$d]=$j, $m[$k-$d]=$v+1 if $yt<=0 && $m[$k-$d]!=$v; 158 $la[$k-$dd]=$j, $m[$k-$dd]=$v+1 if $zt<=0 && $m[$k-$dd]!=$v; 159 } 160 161 # Scan to ensure that no neighboring blocks have been missed by the 162 # outwards-looking logic in the above section 163 for($i=0;$i<$#b;$i+=3) { 164 wl_check($d0+$b[$i]+$b[$i+1]*$d+$b[$i+2]*$dd); 165 } 166 167 # Compute the number of entries where outside blocks do not need to be 168 # consider 169 for($i=$j=0;$i<$#b;$i+=3,$j++) { 170 $k=$d0+$b[$i]+$b[$i+1]*$d+$b[$i+2]*$dd; 171 last if $m[$k+1]!=$v; 172 last if $m[$k+$d]!=$v; 173 last if $m[$k+$dd]!=$v; 174 last if $m[$k-1]!=$v; 175 last if $m[$k-$d]!=$v; 176 last if $m[$k-$dd]!=$v; 177 } 178 print W "\t$j"; 179 180 # Create worklist entry and save to file 181 $j=0; 182 while ($#b>0) { 183 $xt=shift @b;$yt=shift @b;$zt=shift @b; 184 $k=$d0+$xt+$yt*$d+$zt*$dd; 185 $o=0; 186 $o|=1 if $m[$k+1]!=$v && $la[$k+1]==$j; 187 $o^=3 if $m[$k-1]!=$v && $la[$k-1]==$j; 188 $o|=8 if $m[$k+$d]!=$v && $la[$k+$d]==$j; 189 $o^=24 if $m[$k-$d]!=$v && $la[$k-$d]==$j; 190 $o|=64 if $m[$k+$dd]!=$v && $la[$k+$dd]==$j; 191 $o^=192 if $m[$k-$dd]!=$v && $la[$k-$dd]==$j; 192 printf W ",%#x",(($xt+64)|($yt+64)<<7|($zt+64)<<14|$o<<21); 193 $j++; 194 } 195 print W "," unless $ind==$hr*$hr*$hr-1; 196 print W "\n"; 197 198 # Remove list memory 199 undef @a; 200 undef @b; 201} 202 203sub add { 204 if ($m[$d0+@_[0]+$d*@_[1]+$dd*@_[2]]!=$v) { 205 $ac++; 206 push @a,@_[0],@_[1],@_[2]; 207 $m[$d0+@_[0]+$d*@_[1]+$dd*@_[2]]=$v; 208 } 209} 210 211sub dis { 212 $xl=@_[3]+0.3-@_[0];$xh=@_[3]+0.7-@_[0]; 213 $yl=@_[4]+0.3-@_[1];$yh=@_[4]+0.7-@_[1]; 214 $zl=@_[5]+0.3-@_[2];$zh=@_[5]+0.7-@_[2]; 215 $dis=(abs($xl)<abs($xh)?$xl:$xh)**2 216 +(abs($yl)<abs($yh)?$yl:$yh)**2 217 +(abs($zl)<abs($zh)?$zl:$zh)**2; 218 return sqrt $dis; 219} 220 221sub adis { 222 $xco=$yco=$zco=0; 223 $xco=@_[0]-@_[3] if @_[3]>0; 224 $xco=@_[0]-@_[3]-1 if @_[3]<0; 225 $yco=@_[1]-@_[4] if @_[4]>0; 226 $yco=@_[1]-@_[4]-1 if @_[4]<0; 227 $zco=@_[2]-@_[5] if @_[5]>0; 228 $zco=@_[2]-@_[5]-1 if @_[5]<0; 229 return sqrt $xco*$xco+$yco*$yco+$zco*$zco; 230} 231 232sub wl_check { 233 die "Failure in worklist construction\n" if $m[@_[0]+1]<$v||$m[@_[0]-1]<$v 234 ||$m[@_[0]+$d]<$v||$m[@_[0]-$d]<$v 235 ||$m[@_[0]+$dd]<$v||$m[@_[0]-$dd]<$v; 236} 237