Strip extra spaces from code.

[voro++.git] / branches / 2d / src / worklist_gen.pl

#!/usr/bin/perl
# Voro++, a 2D and 3D cell-based Voronoi library
#
# Author   : Chris H. Rycroft (LBL / UC Berkeley)
# Email    : chr@alum.mit.edu
# Date     : August 30th 2011
#
# worklist_gen.pl - this perl script is used to automatically generate the
# worklists of blocks that are stored in worklist.cc, that are used by the
# compute_cell() routine to ensure maximum efficiency

# Each region is divided into a grid of subregions, and a worklist is
# constructed for each. This parameter sets is set to half the number of
# subregions that the block is divided into.
$hr=4;

# This parameter is automatically set to the the number of subregions that the
# block is divided into
$r=$hr*2;

# This parameter sets the total number of block entries in each worklist
$ls=63;

# When the worklists are being constructed, a mask array is made use of. To
# prevent the creation of array elements with negative indices, this parameter
# sets an arbitrary displacement.
$dis=8;

# Constants used mask indexing
$d=2*$dis+1;$d0=(1+$d)*$dis;

use Math::Trig;

# Construct the worklist header file, based on the parameters above
open W,">worklist_2d.hh";
print W <<EOF;
// Voro++, a 2D and 3D cell-based Voronoi library
//
// Author   : Chris H. Rycroft (LBL / UC Berkeley)
// Email    : chr\@alum.mit.edu
// Date     : August 30th 2011

/** \\file worklist_2d.hh
 * \\brief Header file for setting constants used in the block worklists that are
 * used during cell computation.
 *
 * This file is automatically generated by worklist_gen.pl and it is not
 * intended to be edited by hand. */

#ifndef VOROPP_WORKLIST_2D_HH
#define VOROPP_WORKLIST_2D_HH

namespace voro {

/** Each region is divided into a grid of subregions, and a worklist is
# constructed for each. This parameter sets is set to half the number of
# subregions that the block is divided into. */
EOF
print W "const int wl_hgrid_2d=$hr;\n";
print W <<EOF;
/** The number of subregions that a block is subdivided into, which is twice
the value of hgrid. */
EOF
print W "const int wl_fgrid_2d=$r;\n";
print W <<EOF;
/** The total number of worklists, set to the cube of hgrid. */
EOF
printf W "const int wl_hgridsq_2d=%d;\n",$hr*$hr;
print W <<EOF;
/** The number of elements in each worklist. */
EOF
printf W "const int wl_seq_length_2d=%d;\n",$ls+1;
print W <<EOF;

}
#endif
EOF
close W;

# Construct the preamble to the worklist.cc file
open W,">v_base_wl_2d.cc";
print W <<EOF;
// Voro++, a 2D and 3D cell-based Voronoi library
//
// Author   : Chris H. Rycroft (LBL / UC Berkeley)
// Email    : chr\@alum.mit.edu
// Date     : August 30th 2011

/** \\file v_base_wl_2d.cc
 * \\brief The table of block worklists that are used during the cell
 * computation, which is part of the voro_base class.
 *
 * This file is automatically generated by worklist_gen.pl and it is not
 * intended to be edited by hand. */

EOF
printf W "const unsigned int voro_base_2d::wl[wl_seq_length_2d*wl_hgridsq_2d]={\n";

# Now create a worklist for each subregion
for($jj=0;$jj<$hr;$jj++) {
        for($ii=0;$ii<$hr;$ii++) {
                worklist($ii,$jj);
        }
}

# Finish the file and close it
print W "};\n";
close W;

sub worklist {
        $ind=@_[0]+$hr*@_[1];
        $ac=0;$v+=2;
        $xp=$yp=0;
        $x=(@_[0]+0.5)/$r;
        $y=(@_[1]+0.5)/$r;
        $m[$d0]=$v;
        add(1,0,0);add(0,1,0);add(-1,0,0);add(0,-1,0);
        foreach $l (1..$ls) {
                $minwei=1e9;
                foreach (0..$ac-1) {
                        $xt=@a[2*$_];$yt=@a[2*$_+1];
                        $wei=adis($x,$y,$xt,$yt)+0.02*sqrt(($xt-$xp)**2+($yt-$yp)**2);
                        $nx=$_,$minwei=$wei if $wei<$minwei;
                }
                $xp=@a[2*$nx];$yp=@a[2*$nx+1];
                add($xp+1,$yp);add($xp,$yp+1);
                add($xp-1,$yp);add($xp,$yp-1);
#               print "=> $l $xp $yp $zp\n" if $l<4;
                push @b,(splice @a,2*$nx,2);$ac--;
        }

        # Mark all blocks that are on this worklist entry
        $m[$d0]=++$v;
        for($i=0;$i<$#b;$i+=2) {
                $xt=$b[$i];$yt=$b[$i+1];
                #print "$xt $yt\n";
                $m[$d0+$xt+$d*$yt]=$v;
        }

        # Find which neighboring outside blocks need to be marked when
        # considering this block, being as conservative as possible by
        # overwriting the marks, so that the last possible entry that can reach
        # a block is used
        for($i=$j=0;$i<$#b;$i+=2,$j++) {
                $xt=$b[$i];$yt=$b[$i+1];
                $k=$d0+$xt+$yt*$d;
                $la[$k+1]=$j, $m[$k+1]=$v+1 if $xt>=0 && $m[$k+1]!=$v;
                $la[$k+$d]=$j, $m[$k+$d]=$v+1 if $yt>=0 && $m[$k+$d]!=$v;
                $la[$k-1]=$j, $m[$k-1]=$v+1 if $xt<=0 && $m[$k-1]!=$v;
                $la[$k-$d]=$j, $m[$k-$d]=$v+1 if $yt<=0 && $m[$k-$d]!=$v;
        }

        # Scan to ensure that no neighboring blocks have been missed by the
        # outwards-looking logic in the above section
        for($i=0;$i<$#b;$i+=2) {
                wl_check($d0+$b[$i]+$b[$i+1]*$d);
        }

        # Compute the number of entries where outside blocks do not need to be
        # consider
        for($i=$j=0;$i<$#b;$i+=2,$j++) {
                $k=$d0+$b[$i]+$b[$i+1]*$d;
                last if $m[$k+1]!=$v;
                last if $m[$k+$d]!=$v;
                last if $m[$k-1]!=$v;
                last if $m[$k-$d]!=$v;
        }
        print W "\t$j";

        # Create worklist entry and save to file
        $j=0;
        while ($#b>0) {
                $xt=shift @b;$yt=shift @b;
                $k=$d0+$xt+$yt*$d;
                $o=0;
                $o|=1 if $m[$k+1]!=$v && $la[$k+1]==$j;
                $o^=3 if $m[$k-1]!=$v && $la[$k-1]==$j;
                $o|=8 if $m[$k+$d]!=$v && $la[$k+$d]==$j;
                $o^=24 if $m[$k-$d]!=$v && $la[$k-$d]==$j;
                printf W ",%#x",(($xt+64)|($yt+64)<<7|$o<<21);
                $j++;
        }
        print W "," unless $ind==$hr*$hr-1;
        print W "\n";

        # Remove list memory
        undef @a;
        undef @b;
}

sub add {
        if ($m[$d0+@_[0]+$d*@_[1]]!=$v) {
                $ac++;
                push @a,@_[0],@_[1];
                $m[$d0+@_[0]+$d*@_[1]]=$v;
        }
}

sub adis {
        $xco=$yco=0;
        $xco=@_[0]-@_[2] if @_[2]>0;
        $xco=@_[0]-@_[2]-1 if @_[2]<0;
        $yco=@_[1]-@_[3] if @_[3]>0;
        $yco=@_[1]-@_[3]-1 if @_[3]<0;
        return sqrt $xco*$xco+$yco*$yco;
}

sub wl_check {
        die "Failure in worklist construction\n" if $m[@_[0]+1]<$v||$m[@_[0]-1]<$v
                                                  ||$m[@_[0]+$d]<$v||$m[@_[0]-$d]<$v;
}