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