From 9b87c1b77dfb1d5fc55d8c0de6d14941c0ae6c30 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 7 Oct 2015 21:45:36 +0700
Subject: [PATCH] Rewrite neighbor_lists using type(VectorInt)

---
 src/geometry.F90 | 88 +++++++++++++++---------------------------------
 1 file changed, 27 insertions(+), 61 deletions(-)

diff --git a/src/geometry.F90 b/src/geometry.F90
index e922479b1d..75165158bf 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -9,6 +9,7 @@ module geometry
   use particle_header,        only: LocalCoord, Particle
   use particle_restart_write, only: write_particle_restart
   use surface_header
+  use stl_vector,             only: VectorInt
   use string,                 only: to_str
   use tally,                  only: score_surface_current
 
@@ -871,81 +872,46 @@ contains
 
   subroutine neighbor_lists()
 
-    integer :: i          ! index in cells/surfaces array
-    integer :: j          ! index of surface in cell
-    integer :: i_surface  ! index in count arrays
-    integer, allocatable :: count_positive(:) ! # of cells on positive side
-    integer, allocatable :: count_negative(:) ! # of cells on negative side
-    logical :: positive   ! positive side specified in surface list
-    type(Cell), pointer  :: c
+    integer :: i  ! index in cells/surfaces array
+    integer :: j  ! index in region specification
+    integer :: k  ! surface half-space spec
+    type(VectorInt), allocatable :: neighbor_pos(:)
+    type(VectorInt), allocatable :: neighbor_neg(:)
 
     call write_message("Building neighboring cells lists for each surface...", &
-         &4)
+         4)
 
-    allocate(count_positive(n_surfaces))
-    allocate(count_negative(n_surfaces))
-    count_positive = 0
-    count_negative = 0
+    allocate(neighbor_pos(n_surfaces))
+    allocate(neighbor_neg(n_surfaces))
 
     do i = 1, n_cells
-      c => cells(i)
+      do j = 1, size(cells(i)%region)
+        ! Get token from region specification and skip any tokens that
+        ! correspond to operators rather than regions
+        k = cells(i)%region(j)
+        if (abs(k) >= OP_UNION) cycle
 
-      ! loop over each region specification
-      do j = 1, size(c%region)
-        i_surface = c % region(j)
-        positive = (i_surface > 0)
-
-        ! Skip any tokens that correspond to operators rather than regions
-        i_surface = abs(i_surface)
-        if (i_surface >= OP_UNION) cycle
-
-        if (positive) then
-          count_positive(i_surface) = count_positive(i_surface) + 1
+        ! Add this cell ID to neighbor list for k-th surface
+        if (k > 0) then
+          call neighbor_pos(abs(k))%push_back(i)
         else
-          count_negative(i_surface) = count_negative(i_surface) + 1
+          call neighbor_neg(abs(k))%push_back(i)
         end if
       end do
     end do
 
-    ! allocate neighbor lists for each surface
     do i = 1, n_surfaces
-      if (count_positive(i) > 0) then
-        allocate(surfaces(i)%obj%neighbor_pos(count_positive(i)))
-      end if
-      if (count_negative(i) > 0) then
-        allocate(surfaces(i)%obj%neighbor_neg(count_negative(i)))
-      end if
+      ! Copy positive neighbors to Surface instance
+      j = neighbor_pos(i)%size()
+      allocate(surfaces(i)%obj%neighbor_pos(j))
+      surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:j)
+
+      ! Copy negative neighbors to Surface instance
+      j = neighbor_neg(i)%size()
+      allocate(surfaces(i)%obj%neighbor_neg(j))
+      surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:j)
     end do
 
-    count_positive = 0
-    count_negative = 0
-
-    ! loop over all cells
-    do i = 1, n_cells
-      c => cells(i)
-
-      ! loop through the region specification
-      do j = 1, size(c%region)
-        i_surface = c % region(j)
-        positive = (i_surface > 0)
-
-        ! Skip any tokens that correspond to operators rather than regions
-        i_surface = abs(i_surface)
-        if (i_surface >= OP_UNION) cycle
-
-        if (positive) then
-          count_positive(i_surface) = count_positive(i_surface) + 1
-          surfaces(i_surface)%obj%neighbor_pos(count_positive(i_surface)) = i
-        else
-          count_negative(i_surface) = count_negative(i_surface) + 1
-          surfaces(i_surface)%obj%neighbor_neg(count_negative(i_surface)) = i
-        end if
-      end do
-    end do
-
-    deallocate(count_positive)
-    deallocate(count_negative)
-
   end subroutine neighbor_lists
 
 !===============================================================================