diff --git a/examples/carotid_web/D/.gittouch b/examples/carotid_web/D/.gittouch
new file mode 100644
index 0000000..e69de29
diff --git a/examples/carotid_web/Input/carotid.e b/examples/carotid_web/Input/carotid.e
new file mode 100644
index 0000000..5fe6cba
Binary files /dev/null and b/examples/carotid_web/Input/carotid.e differ
diff --git a/examples/carotid_web/Input/tube.in b/examples/carotid_web/Input/tube.in
new file mode 100644
index 0000000..a5f5e18
--- /dev/null
+++ b/examples/carotid_web/Input/tube.in
@@ -0,0 +1,30 @@
+0.44 !0.04              ! alpha_Ewld !changed for big tube
+1.E-3                   ! eps_Ewld
+8                       ! PBspln_Ewld
+
+3                       ! nCellTypes
+1 1 1                   ! viscRat -- for RBCS, approx 5, for Sickles, approximated as 7 from https://doi.org/10.1016/0167-4889(90)90025-9
+1. 1. 1.                ! refRad
+.false.                 ! Deflate
+
+0.0     !
+0.0     !
+0.0     !
+
+15000                       ! Nt
+0.0007                  ! Ts
+
+100                      ! cell_out
+-10000                   ! wall_out
+-10                      ! pGrad_out
+-10                      ! flow_out
+-1                       ! ftotout
+100                      ! restart_out
+
+'D/restart.LATEST.dat'
+
+0.02                ! epsDist
+10.                 ! forceCoef
+10.                 ! viscRatThresh
+.false.             ! rigidsep
+0. 0. 0. 0.  		! fmags
diff --git a/examples/carotid_web/Input/vertical_web.e b/examples/carotid_web/Input/vertical_web.e
new file mode 100644
index 0000000..47e5593
Binary files /dev/null and b/examples/carotid_web/Input/vertical_web.e differ
diff --git a/examples/carotid_web/Input/web.e b/examples/carotid_web/Input/web.e
new file mode 100644
index 0000000..a72fd29
Binary files /dev/null and b/examples/carotid_web/Input/web.e differ
diff --git a/examples/carotid_web/Makefile b/examples/carotid_web/Makefile
new file mode 100755
index 0000000..487e508
--- /dev/null
+++ b/examples/carotid_web/Makefile
@@ -0,0 +1,22 @@
+include ../../Makefile.in
+
+LIB = $(COMMON_LIB) $(SPHPK_LIB) $(FFTW_LIB) $(LAPACK95_LIB) $(PETSC_LIB) $(NETCDF_LIB) $(MKL_LIB) $(STATIC)
+EXECUTABLES = initcond tube
+
+all : $(EXECUTABLES)
+
+lib : $(wildcard $(WORK_DIR)/common/*.F90)
+	make -C $(WORK_DIR)/common
+
+$(EXECUTABLES) : % : %.o $(wildcard $(COMMON_DIR)/*.F90)
+	make lib
+	$(FC) $(LDFLAGS) -o $@ $< $(LIB)
+
+clean :
+	-rm -f $(EXECUTABLES) *.o *.mod *.a core
+
+# Dependency
+.depend : $(wildcard *.F90)
+	$(MAKEDEPEND_BIN) $^ > .depend
+
+include .depend
diff --git a/examples/carotid_web/initcond.F90 b/examples/carotid_web/initcond.F90
new file mode 100644
index 0000000..815fc09
--- /dev/null
+++ b/examples/carotid_web/initcond.F90
@@ -0,0 +1,407 @@
+program randomized_cell_gen
+
+  use ModDataTypes
+  use ModDataStruct
+  use ModRbc
+  use ModWall
+  use ModConf
+  use ModData
+  use ModIO
+  use ModBasicMath
+
+  integer, parameter :: ranseed = 112
+
+  !initial condition setup parameters
+  real(WP), parameter :: hematocrit = 0.2
+  real(WP) :: tuber = 4
+  real(WP), parameter :: tubelen = 30
+
+  integer :: nrbcMax ! how many cells
+
+  type(t_Rbc), pointer :: rbc
+  type(t_Wall), pointer :: wall
+  character(CHRLEN) :: fn
+  integer :: zmax
+  integer :: i
+  real(WP) :: th, actlen
+  real(WP) :: clockBgn, clockEnd
+  real(WP) :: offset(3)
+  integer :: nodeNum, numNodes
+
+  call InitMPI
+  call MPI_Comm_Rank(MPI_COMM_WORLD, nodeNum, ierr)
+  call MPI_Comm_Size(MPI_COMM_WORLD, numNodes, ierr)
+  nodeNum = nodeNum + 1
+
+  !calculate number of cells for the defined hematocrit, assuming all blood cells are healthy RBCs for volume
+  !hematocrit = 4 * nrbc / (tube_radius^2 * tube_length)
+  nrbcMax = ((3*(tubelen*tuber**2*hematocrit))/4)
+  ! nrbcMax = 500
+
+  write (*, *) "Num RBCs in simulation is ", nrbcMax
+
+  !set other initialization params
+  vBkg(1:2) = 0.; vBkg(3) = 8.
+  Nt = 0; time = 0.
+
+  !Create wall
+  nwall = 2
+  allocate (walls(nwall))
+  wall => walls(1)
+  call ReadWallMesh('Input/carotid.e', wall)
+  wall%f = 0.
+  wall => walls(2)
+  call ReadWallMesh('Input/web.e', wall)
+  wall%f = 0.
+
+  wall => walls(1)
+
+  if (rootWorld) write (*, *) walls(1)%nvert, walls(2)%nvert
+  actlen = tubelen
+
+  !recenter walls to be strictly positive
+  call recenterWalls
+
+  !set periodic boundary box based on tube shape
+  Lb(1) = maxval(wall%x(:, 1)) + 0.5
+  Lb(2) = maxval(wall%x(:, 2)) + 0.5
+  Lb(3) = tubelen
+
+  !Write Walls Out to Wall TecPlot File
+  write (fn, FMT=fn_FMT) 'D/', 'wall', 0, '.dat'
+  call WriteManyWalls(fn, nwall, walls)
+
+  !for each cell to add
+  allocate (rbcs(nrbcMax))
+  nrbc = 0
+
+  do i = 1, nrbcMax
+    if (rootWorld) write (*, *) "Adding Cell #", nrbc + 1
+    clockBgn = MPI_Wtime()
+    call place_cell(1)
+    clockEnd = MPI_Wtime()
+
+    nrbc = nrbc + 1
+    if (rootWorld) write (*, *) "Cell #", nrbc, " Placed; Time Cost = ", clockEnd - clockBgn
+
+    !Write Cells Out to Cell TecPlot File
+    if (modulo(i, 10) .eq. 0) then
+      write (fn, FMT=fn_FMT) 'D/', 'x', 0, '.dat'
+      call WriteManyRBCs(fn, nrbc, rbcs)
+      !Write restart blob
+      write (fn, FMT=fn_FMT) 'D/', 'restart', 0, '.dat'
+      call WriteRestart(fn, Nt0, time)
+      fn = 'D/restart.LATEST.dat'
+      call WriteRestart(fn, Nt0, time)
+    end if
+  end do
+
+  write (fn, FMT=fn_FMT) 'D/', 'x', 0, '.dat'
+  call WriteManyRBCs(fn, nrbc, rbcs)
+  !Write restart blob
+  write (fn, FMT=fn_FMT) 'D/', 'restart', 0, '.dat'
+  call WriteRestart(fn, Nt0, time)
+  fn = 'D/restart.LATEST.dat'
+  call WriteRestart(fn, Nt0, time)
+
+  if (nrbcMax .ne. 0) deallocate (rbcs)
+  if (nwall .ne. 0) deallocate (walls)
+  call FinalizeMPI
+  stop
+
+contains
+
+!offset everything to be positive
+  subroutine recenterWalls
+    integer :: i, ii
+    type(t_wall), pointer :: wall
+    real(WP) :: offset(3)
+
+    wall => walls(1)
+
+    offset(1) = -1*minval(wall%x(:, 1))
+    offset(2) = -1*minval(wall%x(:, 2))
+    offset(3) = -1*minval(wall%x(:, 3))
+
+    !shift wall
+    do ii = 1, 3
+      wall => walls(1)
+      wall%x(:, ii) = wall%x(:, ii) + offset(ii)
+      wall => walls(2)
+      wall%x(:, ii) = wall%x(:, ii) + offset(ii)
+    end do
+
+  end subroutine recenterWalls
+
+  recursive subroutine choose_point(pt)
+    real(WP) :: pt(3)
+
+    type(t_wall), pointer :: wall
+    real(WP) :: rad, len, maxN, minN
+    real(WP) :: threshold
+    integer :: i
+
+    wall => walls(1)
+    len = RandomNumber(ranseed)*actlen
+    threshold = 0.2
+    maxN = -1E9
+    minN = 1E9
+
+    !figure out the radius based on the location in the tube
+    do i = 1, wall%nvert
+      if (abs(wall%x(i, 3) - len) .le. threshold) then
+        maxN = maxval((/maxN, wall%x(i, 1), wall%x(i, 2)/))
+        minN = minval((/minN, wall%x(i, 1), wall%x(i, 2)/))
+      end if
+    end do
+
+    rad = (maxN - minN)/2
+
+    pt(2) = RandomNumber(ranseed)*2*PI
+    pt(3) = sqrt(RandomNumber(ranseed))*(rad)
+    pt(1) = pt(3)*cos(pt(2)) + rad + minN
+    pt(2) = pt(3)*sin(pt(2)) + rad + minN
+    pt(3) = len
+
+    !edge case for the web
+    wall => walls(2)
+    maxN = -1E9
+    minN = 1E9
+    threshold = 0.5
+    do i = 1, wall%nvert
+      if (abs(wall%x(i, 3) - len) .le. threshold) then
+        minN = min(minN, wall%x(i, 2))
+        maxN = max(maxN, wall%x(i, 2))
+      end if
+    end do
+
+    if (minN .le. pt(2) .and. pt(2) .le. maxN) then !.and. maxN .ge. pt(2)) then
+      ! if (rootWorld) write (*,*) "point is inside web"
+      pt = 0
+      call choose_point(pt)
+    end if
+
+  end subroutine choose_point
+
+!find an open spot in the simulation to place a cell
+  subroutine place_cell(celltype)
+
+    type(t_Rbc) :: newcell
+    type(t_Rbc), pointer :: cell
+    real(WP) :: tmp_xc(3)
+
+    integer :: celli, ii, ierr
+    integer :: nlat0
+    logical :: place_success_loc, place_success_glb
+
+    integer :: celltype
+
+    nlat0 = 12
+
+    !create template newcell
+    select case (celltype)
+    case (1)
+      call RBC_Create(newcell, nlat0)
+      call RBC_MakeBiconcave(newcell, 1.)
+    case (2)
+      call RBC_Create(newcell, nlat0)
+      call RBC_MakeLeukocyte(newcell, 1.)
+    case (3)
+      call ImportReadRBC('Input/SickleCell.dat', newcell)
+    case default
+      stop "bad cellcase"
+    end select
+    newcell%celltype = celltype
+
+    place_success_glb = .false.
+    tmp_xc = 0.
+
+    !choose a location & orientation for newcell to fit into the simulation
+    do while (.not. place_success_glb)
+
+      !reset cell position
+      do ii = 1, 3
+        newcell%x(:, :, ii) = newcell%x(:, :, ii) - tmp_xc(ii)
+      end do
+
+      !randomly rotate cell
+      call rotate_cell(newcell)
+
+      !randomly select a tmp_xc
+      call choose_point(tmp_xc)
+
+      !shift newcell by the xc
+      do ii = 1, 3
+        newcell%x(:, :, ii) = newcell%x(:, :, ii) + tmp_xc(ii)
+      end do
+
+      place_success_glb = .true.
+      newcell%x(:, :, 1) = newcell%x(:, :, 1) - tuber
+      newcell%x(:, :, 2) = newcell%x(:, :, 2) - tuber
+      if (any(newcell%x(:, :, 3) .ge. tubelen .or. newcell%x(:, :, 3) .lt. 0)) then
+        place_success_loc = .false.
+        place_success_glb = .false.
+      end if
+      newcell%x(:, :, 1) = newcell%x(:, :, 1) + tuber
+      newcell%x(:, :, 2) = newcell%x(:, :, 2) + tuber
+      if (.not. place_success_glb) cycle
+
+      !check if cell collides with wall
+      place_success_loc = .not. check_wall_collision(newcell)
+      place_success_glb = .true.
+      call MPI_AllReduce(place_success_loc, place_success_glb, 1, MPI_Logical, MPI_LAND, MPI_COMM_WORLD, ierr)
+      if (.not. place_success_glb) cycle
+
+      !if no wall collision, check if cell collides with other cells
+      celli = 1
+      do while (place_success_loc .and. celli .le. nrbc)
+        cell => rbcs(celli)
+        place_success_loc = .not. check_cell_collision(cell, newcell)
+        celli = celli + 1
+      end do
+
+      place_success_glb = .true.
+      call MPI_AllReduce(place_success_loc, place_success_glb, 1, MPI_Logical, MPI_LAND, MPI_COMM_WORLD, ierr)
+
+    end do !while not place_success
+
+    !the cell placement location was successful
+    rbcs(nrbc + 1) = newcell
+  end subroutine place_cell
+
+!helper for place cell, chooses a random orientation and rotates the cell
+  subroutine rotate_cell(cell)
+    type(t_Rbc) :: cell
+
+    real(WP) :: v1, v2, vsq
+    real(WP) :: zv(3), rc(3), rotmat(3, 3)
+
+    integer :: i, j
+
+    !randomly generate unit vector
+    vsq = 10
+    do while (vsq .ge. 1)
+      v1 = RandomNumber(ranseed)*2 - 1
+      v2 = RandomNumber(ranseed)*2 - 1
+      vsq = v1**2 + v2**2
+    end do
+    zv(1) = v1*2*sqrt(1 - vsq)
+    zv(2) = v2*2*sqrt(1 - vsq)
+    zv(3) = 1 - (2*vsq)
+
+    !generate rotation matrix
+    rotmat = RotateMatrix(zv) !(/ 0., 1., 0./))
+
+    !rotate cell with rotmat
+    !each point: x = Rx
+    forall (i=1:cell%nlat, j=1:cell%nlon)
+      cell%x(i, j, :) = reshape(matmul(rotmat, cell%x(i, j, :)), (/3/))
+    end forall
+
+  end subroutine rotate_cell
+
+!helper for place_cell, checks if cell1 intersects/collides with wall
+  logical function check_wall_collision(cell)
+    type(t_Rbc) :: cell
+
+    integer :: i, j, i2, j2
+    real(WP) :: c1p(3), c2p(3), sp(3)
+
+    real(WP) :: threshold
+    threshold = 0.3
+
+    check_wall_collision = .false.
+
+    !each point in cell1
+    do i = 1, cell%nlat
+    do j = 1, cell%nlon
+
+      !define 2 diagonal points c1p & c2p for collision detection
+      c1p = cell%x(i, j, :)
+      c2p = cell%x(modulo(i, cell%nlat) + 1, modulo(j, cell%nlon) + 1, :)
+
+      !each wall
+      do i2 = 1, nwall
+        !each point in the wall
+        do j2 = nodeNum, walls(i2)%nvert, numNodes
+
+          sp = walls(i2)%x(j2, :)
+
+          if (VecNorm(c1p - sp) .lt. threshold) then
+            check_wall_collision = .true.
+            return
+          end if
+        end do !j2
+      end do !i2
+
+    end do !j
+    end do !i
+  end function check_wall_collision
+
+! helper for place_cell, checks if cell1 and cell2 intersect/collide
+! creates miniature bounding boxes along cell meshes to do AABB-AABB collision detection
+! should not have any false-negatives (return false even if the cells are intersecting)
+! may have false-positives (might return true if cells aren't intersecting, but are very close)
+  logical function check_cell_collision(cell1, cell2)
+    type(t_Rbc) :: cell1, cell2
+
+    integer :: i, j, i2, j2, ii
+    integer :: p_it
+
+    real(WP) :: p1(3), p2(3)
+    real(WP) :: b1(3, 2), b2(3, 2)
+
+    check_cell_collision = .false.
+
+    !first check the centers; if distance > 4 then we are sure that they don't collide
+    if (VecNorm(cell1%xc - cell2%xc) .ge. 4) return
+
+    !each point in cell1
+    do i = 1, cell1%nlat
+    do j = 1, cell1%nlon
+
+      ! create first AABB from patch of cell1
+      p1 = cell1%x(i, j, :)
+      p2 = cell1%x(modulo(i, cell1%nlat) + 1, modulo(j, cell1%nlon) + 1, :)
+
+      do ii = 1, 3
+        b1(ii, 1) = min(p1(ii), p2(ii))
+        b1(ii, 2) = max(p1(ii), p2(ii))
+      end do
+
+      !each point in cell2
+      do p_it = nodeNum, cell2%nlat*cell2%nlon, numNodes
+        ! do i2 = 1, cell2%nlat
+        ! do j2 = 1, cell2%nlon
+        i2 = (p_it/cell2%nlon) + 1
+        j2 = modulo(p_it, cell2%nlon) + 1
+
+        !create second AABB from patch of cell2
+        p1 = cell2%x(i2, j2, :)
+        p2 = cell2%x(modulo(i2, cell2%nlat) + 1, modulo(j2, cell2%nlon) + 1, :)
+        do ii = 1, 3
+          b2(ii, 1) = min(p1(ii), p2(ii))
+          b2(ii, 2) = max(p1(ii), p2(ii))
+        end do
+
+        !check for AABB collision, we find AABB intersection iff X, Y, and Z intervals all overlap
+        if ( &
+          b1(1, 1) .le. b2(1, 2) &
+          .and. b1(1, 2) .ge. b2(1, 1) &
+          .and. b1(2, 1) .le. b2(2, 2) &
+          .and. b1(2, 2) .ge. b2(2, 1) &
+          .and. b1(3, 1) .le. b2(3, 2) &
+          .and. b1(3, 2) .ge. b2(3, 1) &
+          ) then
+          check_cell_collision = .true.
+          return
+        end if
+
+        ! end do !j2
+        ! end do !i2
+      end do !p_it
+    end do !j
+    end do !i
+  end function check_cell_collision
+
+end program randomized_cell_gen
diff --git a/examples/carotid_web/tube.F90 b/examples/carotid_web/tube.F90
new file mode 100755
index 0000000..9e2f7c3
--- /dev/null
+++ b/examples/carotid_web/tube.F90
@@ -0,0 +1,179 @@
+! cells in cylindrical tubes
+program cells_in_tube
+
+  use ModDataTypes
+  use ModDataStruct
+  use ModRBC
+  use ModWall
+  use ModConf
+  use ModData
+  use ModTimeInt
+  use ModIO
+  use ModBasicMath
+  use ModQuadRule
+
+  implicit none
+  integer :: cutoff
+  character(CHRLEN) :: fn
+
+#include "../../petsc_include.h"
+
+  call InitAll
+
+  call TimeInt_Euler
+
+  call FinalizeAll
+
+  stop
+
+contains
+
+!**********************************************************************
+  subroutine InitAll
+
+    ! System intialization
+    call InitMPI()
+    call GaussQuad_Init
+    call ReadConfig('Input/tube.in')
+    call InitSystem
+
+    ! Note: SetEwaldPrms and DomainDecomp must be called after InitSystem
+    call SetEwaldPrms
+    call DomainDecomp
+    call GlobData_Init
+
+    ! Prepare for time integration
+    call IO_Init
+    call TimeInt_Init
+
+  end subroutine InitAll
+
+!**********************************************************************
+  subroutine FinalizeAll
+
+    call IO_Finalize
+    call TimeInt_Finalize
+
+    call GlobData_Finalize
+    call FinalizeSystem
+
+    call GaussQuad_Finalize
+    call FinalizeMPI
+
+  end subroutine FinalizeAll
+
+!**********************************************************************
+  subroutine InitSystem
+
+    integer :: irbc, iwall
+    type(t_Rbc), pointer :: rbc, rbcRef
+    type(t_Wall), pointer :: wall
+    integer :: nlat0
+    real(WP) :: radEqv
+    integer :: ierr
+
+    call ReadRestart(restart_file)
+
+    ! Reference cells
+    allocate (rbcRefs(3))
+
+    if (nrbc > 0) then
+      radEqv = 1.
+      nlat0 = rbcs(1)%nlat0
+
+      rbcRef => rbcRefs(1)
+      call RBC_Create(rbcRef, nlat0)
+      call RBC_MakeBiconcave(rbcRef, radEqv)
+      call RBC_ComputeGeometry(rbcRef)
+
+      rbcRef => rbcRefs(2)
+      call RBC_Create(rbcRef, nlat0)
+      call RBC_MakeLeukocyte(rbcRef, radEqv)
+      call RBC_ComputeGeometry(rbcRef)
+
+      rbcRef => rbcRefs(3)
+      call ImportReadRBC('Input/SickleCell.dat', rbcRef)
+      call RBC_ComputeGeometry(rbcRef)
+
+    end if
+
+    ! Wall periodic boundary condition
+    do iwall = 1, nwall
+      wall => walls(iwall)
+      call Wall_Build_V2V(wall, Lb)
+    end do ! iwall
+
+    ! Mechanical properties of cells and walls
+    do irbc = 1, nrbc
+      rbc => rbcs(irbc)
+      select case (rbc%celltype)
+      case (1)
+        rbc%ES = 12.4
+        rbc%ED = 200.
+        rbc%EB = 6.69D-2
+        ! Mechanical properties of leukocytes
+        ! Calculated according to section 6 of reference paper with WBC radius of 1.4 sim units
+        ! and Es* scaled by 10^2
+        ! Reference:
+        !   Zhao, H., Isfahani, A. H., Olson, L. N., & Freund, J. B. (2010).
+        !   A spectral boundary integral method for flowing blood cells.
+        !   Journal of Computational Physics, 229(10), 3726-3744.
+      case (2)
+        rbc%ES = 887
+        rbc%ED = 200.
+        rbc%EB = 2.44D-2
+
+        ! Mechanical properties for sickle cell roughly determined to be
+        ! Es = (20 / 7.1) * Es for a healthy RBC (case(1) cell)
+        ! Ed = (49.4 / 15.4) * Ed for healthy RBC (case(1) cell)
+        ! Eb = (19.5 / 5.7) * Eb for a healthy RBC (case(1) cell)
+        ! Reference:
+        !   HeeSu Byun, Timothy R. Hillman, John M. Higgins, Monica Diez-Silva, Zhangli Peng, Ming Dao, Ramachandra R. Dasari, Subra Suresh, YongKeun Park
+        !   Optical measurement of biomechanical properties of individual erythrocytes
+        !   Acta Biomaterialia, Volume 8, Issue 11, 2012, Pages 4130-4138,
+        !   https://doi.org/10.1016/j.actbio.2012.07.011
+      case (3)
+        rbc%ES = 12.4*20/7.1
+        rbc%ED = 200*49.4/15.4
+        rbc%EB = 6.69D-2*19.5/5.7
+      case default
+        stop "bad cellcase"
+      end select
+    end do ! irbc
+
+    do iwall = 1, nwall
+      wall => walls(iwall)
+    end do ! iwall
+
+    ! Background velocity
+    vbkg(1:2) = 0.
+    vbkg(3) = 8.
+    print *, vbkg
+
+  end subroutine InitSystem
+
+!**********************************************************************
+  subroutine FinalizeSystem
+
+    integer :: irbc, iwall
+    type(t_Rbc), pointer :: rbc
+    type(t_Wall), pointer :: wall
+
+    do irbc = 1, nrbc
+      rbc => rbcs(irbc)
+      call RBC_Destroy(rbc)
+    end do ! irbc
+
+    do iwall = 1, nwall
+      wall => walls(iwall)
+      call Wall_Destroy(wall)
+    end do ! iwall
+
+    if (nrbc > 0) deallocate (rbcs)
+    if (nwall > 0) deallocate (walls)
+
+  end subroutine FinalizeSystem
+
+!**********************************************************************
+
+end program
diff --git a/examples/case/Makefile b/examples/case/Makefile
index 9bb1e09..34d05ee 100755
--- a/examples/case/Makefile
+++ b/examples/case/Makefile
@@ -1,7 +1,7 @@
 include ../../Makefile.in
 
 LIB = $(COMMON_LIB) $(SPHPK_LIB) $(FFTW_LIB) $(LAPACK95_LIB) $(PETSC_LIB) $(NETCDF_LIB) $(MKL_LIB) $(STATIC)
-EXECUTABLES = tube initcond field
+EXECUTABLES = tube initcond
 
 all : $(EXECUTABLES)
 
diff --git a/examples/randomized_case/Makefile b/examples/randomized_case/Makefile
index 19c68e1..487e508 100755
--- a/examples/randomized_case/Makefile
+++ b/examples/randomized_case/Makefile
@@ -1,7 +1,7 @@
 include ../../Makefile.in
 
 LIB = $(COMMON_LIB) $(SPHPK_LIB) $(FFTW_LIB) $(LAPACK95_LIB) $(PETSC_LIB) $(NETCDF_LIB) $(MKL_LIB) $(STATIC)
-EXECUTABLES = initcond_random tube
+EXECUTABLES = initcond tube
 
 all : $(EXECUTABLES)
 
diff --git a/examples/randomized_case/initcond_random.F90 b/examples/randomized_case/initcond.F90
similarity index 78%
rename from examples/randomized_case/initcond_random.F90
rename to examples/randomized_case/initcond.F90
index 1454a6a..bd6715b 100644
--- a/examples/randomized_case/initcond_random.F90
+++ b/examples/randomized_case/initcond.F90
@@ -26,12 +26,19 @@ program randomized_cell_gen
   real(WP) :: th, actlen
   real(WP) :: clockBgn, clockEnd
 
+  integer :: numNodes, nodeNum
+
   call InitMPI
+  call MPI_Comm_Rank(MPI_COMM_WORLD, nodeNum, ierr)
+  call MPI_Comm_Size(MPI_COMM_WORLD, numNodes, ierr)
+  nodeNum = nodeNum + 1
 
   !calculate number of cells for the defined hematocrit, assuming all blood cells are healthy RBCs for volume
   !hematocrit = 4 * nrbc / (tube_radius^2 * tube_length)
   nrbcMax = ((3*(tubelen*tuber**2*hematocrit))/4)
 
+  if (rootWorld) write (*, *) "Num RBCs in simulation is ", nrbcMax
+
   !set periodic boundary box based on tube shape
   Lb(1) = tuber*2 + 0.5
   Lb(2) = tuber*2 + 0.5
@@ -138,9 +145,9 @@ subroutine place_cell(celltype)
     type(t_Rbc), pointer :: cell
     real(WP) :: tmp_xc(3)
 
-    integer :: celli, ii
+    integer :: celli, ii, ierr
     integer :: nlat0
-    logical :: place_success
+    logical :: place_success_loc, place_success_glb
 
     integer :: celltype
 
@@ -161,10 +168,16 @@ subroutine place_cell(celltype)
     end select
     newcell%celltype = celltype
 
-    place_success = .false.
+    place_success_glb = .false.
+    tmp_xc = 0.
 
     !choose a location & orientation for newcell to fit into the simulation
-    do while (.not. place_success)
+    do while (.not. place_success_glb)
+
+      !reset cell position
+      do ii = 1, 3
+        newcell%x(:, :, ii) = newcell%x(:, :, ii) - tmp_xc(ii)
+      end do
 
       !randomly rotate cell
       call rotate_cell(newcell)
@@ -181,32 +194,29 @@ subroutine place_cell(celltype)
         newcell%x(:, :, ii) = newcell%x(:, :, ii) + tmp_xc(ii)
       end do
 
-      !check if newcell collides with wall, assume cylinder wall for now
+      !check if newcell collides with wall, assume cylinder wall here
       if (any(abs(newcell%x(:, :, 3)) .ge. tubelen/2 .or. &
               newcell%x(:, :, 1)**2 + newcell%x(:, :, 2)**2 .ge. tuber**2)) then
 
-        !placement collides with wall, so try again
-        do ii = 1, 3
-          newcell%x(:, :, ii) = newcell%x(:, :, ii) - tmp_xc(ii)
-        end do
+        !placement collides with wall, try again
         cycle
       end if
 
-      !set boolean success to true (reset to false if we find a collision)
-      place_success = .true.
-      do celli = 1, nrbc
-        cell => rbcs(celli)
+      !if your wall isn't a cylinder, you can do this instead
+      ! place_success_loc = .not. check_wall_collision(newcell)
+      ! place_success_glb = .true.
+      ! call MPI_AllReduce(place_success_loc, place_success_glb, 1, MPI_Logical, MPI_LAND, MPI_COMM_WORLD, ierr)
+      ! if (.not. place_success_glb) cycle
 
-        !we find a collision
-        if (check_cell_collision(cell, newcell)) then
-          !placement collides with existing cell, so try again
-          do ii = 1, 3
-            newcell%x(:, :, ii) = newcell%x(:, :, ii) - tmp_xc(ii)
-          end do
-          place_success = .false.
-          exit
-        end if
-      end do !celli
+      place_success_loc = .true.
+      celli = 1
+      do while (place_success_loc .and. celli .le. nrbc)
+        cell => rbcs(celli)
+        place_success_loc = .not. check_cell_collision(cell, newcell)
+        celli = celli + 1
+      end do
+      place_success_glb = .true.
+      call MPI_AllReduce(place_success_loc, place_success_glb, 1, MPI_Logical, MPI_LAND, MPI_COMM_WORLD, ierr)
 
     end do !while not place_success
 
@@ -248,9 +258,10 @@ end subroutine rotate_cell
 !helper for place_cell, checks if cell1 intersects/collides with wall
   logical function check_wall_collision(cell)
     type(t_Rbc) :: cell
+    real(WP), parameter :: threshold = 0.2 !'magic' adjustable number cell-to-wall distance tolerance
 
     integer :: i, j, i2, j2
-    real(WP) :: c1p(3), c2p(3), sp(3)
+    real(WP) :: cp(3), sp(3)
 
     check_wall_collision = .false.
 
@@ -258,28 +269,21 @@ logical function check_wall_collision(cell)
     do i = 1, cell%nlat
     do j = 1, cell%nlon
 
-      !define 2 diagonal points c1p & c2p for collision detection
-      c1p = cell%x(i, j, :)
-      c2p = cell%x(modulo(i, cell%nlat) + 1, modulo(j, cell%nlon) + 1, :)
+      !each point in cell
+      cp = cell%x(i, j, :)
 
       !each wall
       do i2 = 1, nwall
         !each point in the wall
-        do j2 = 1, walls(i2)%nvert
+        do j2 = nodeNum, walls(i2)%nvert, numNodes
 
           sp = walls(i2)%x(j2, :)
 
-          !check if the cell2 point lies inside the cube delineated by c1p & c2p
-          if ( &
-            ((c1p(1) .le. sp(1) .and. sp(1) .le. c2p(1)) .or. (c2p(1) .le. sp(1) .and. sp(1) .le. c1p(1))) .and. &  !x direction overlap
-            ((c1p(2) .le. sp(2) .and. sp(2) .le. c2p(2)) .or. (c2p(2) .le. sp(2) .and. sp(2) .le. c1p(2))) .and. &  !y direction overlap
-            ((c1p(3) .le. sp(3) .and. sp(3) .le. c2p(3)) .or. (c2p(3) .le. sp(3) .and. sp(3) .le. c1p(3))) &        !z direction overlap
-            ) then
-            ! since we found a collision, return true
+          !just do a distance check
+          if (VecNorm(sp - cp) .le. threshold) then
             check_wall_collision = .true.
             return
           end if
-
         end do !j2
       end do !i2
 
@@ -295,12 +299,16 @@ logical function check_cell_collision(cell1, cell2)
     type(t_Rbc) :: cell1, cell2
 
     integer :: i, j, i2, j2, ii
+    integer :: p_it
 
     real(WP) :: p1(3), p2(3)
     real(WP) :: b1(3, 2), b2(3, 2)
 
     check_cell_collision = .false.
 
+    !first check the centers; if distance > 4 then we are sure that they don't collide
+    if (VecNorm(cell1%xc - cell2%xc) .ge. 4) return
+
     !each point in cell1
     do i = 1, cell1%nlat
     do j = 1, cell1%nlon
@@ -315,8 +323,11 @@ logical function check_cell_collision(cell1, cell2)
       end do
 
       !each point in cell2
-      do i2 = 1, cell2%nlat
-      do j2 = 1, cell2%nlon
+      do p_it = nodeNum, cell2%nlat*cell2%nlon, numNodes
+        ! do i2 = 1, cell2%nlat
+        ! do j2 = 1, cell2%nlon
+        i2 = (p_it/cell2%nlon) + 1
+        j2 = modulo(p_it, cell2%nlon) + 1
 
         !create second AABB from patch of cell2
         p1 = cell2%x(i2, j2, :)
@@ -326,7 +337,7 @@ logical function check_cell_collision(cell1, cell2)
           b2(ii, 2) = max(p1(ii), p2(ii))
         end do
 
-        !check for AABB collision, we find AABB intersection iff X, Y, and Z intervals all overlap
+        !AABB collision check, we find AABB intersection iff X, Y, and Z intervals all overlap
         if ( &
           b1(1, 1) .le. b2(1, 2) &
           .and. b1(1, 2) .ge. b2(1, 1) &
@@ -339,10 +350,12 @@ logical function check_cell_collision(cell1, cell2)
           return
         end if
 
-      end do !j2
-      end do !i2
+        ! end do !j2
+        ! end do !i2
+      end do !p_it
     end do !j
     end do !i
+
   end function check_cell_collision
 
 end program randomized_cell_gen
diff --git a/examples/velocity_field/D/.gittouch b/examples/velocity_field/D/.gittouch
new file mode 100644
index 0000000..e69de29
diff --git a/examples/velocity_field/Input/tube.in b/examples/velocity_field/Input/tube.in
new file mode 100644
index 0000000..a5f5e18
--- /dev/null
+++ b/examples/velocity_field/Input/tube.in
@@ -0,0 +1,30 @@
+0.44 !0.04              ! alpha_Ewld !changed for big tube
+1.E-3                   ! eps_Ewld
+8                       ! PBspln_Ewld
+
+3                       ! nCellTypes
+1 1 1                   ! viscRat -- for RBCS, approx 5, for Sickles, approximated as 7 from https://doi.org/10.1016/0167-4889(90)90025-9
+1. 1. 1.                ! refRad
+.false.                 ! Deflate
+
+0.0     !
+0.0     !
+0.0     !
+
+15000                       ! Nt
+0.0007                  ! Ts
+
+100                      ! cell_out
+-10000                   ! wall_out
+-10                      ! pGrad_out
+-10                      ! flow_out
+-1                       ! ftotout
+100                      ! restart_out
+
+'D/restart.LATEST.dat'
+
+0.02                ! epsDist
+10.                 ! forceCoef
+10.                 ! viscRatThresh
+.false.             ! rigidsep
+0. 0. 0. 0.  		! fmags
diff --git a/examples/velocity_field/Makefile b/examples/velocity_field/Makefile
new file mode 100755
index 0000000..8d1211e
--- /dev/null
+++ b/examples/velocity_field/Makefile
@@ -0,0 +1,22 @@
+include ../../Makefile.in
+
+LIB = $(COMMON_LIB) $(SPHPK_LIB) $(FFTW_LIB) $(LAPACK95_LIB) $(PETSC_LIB) $(NETCDF_LIB) $(MKL_LIB) $(STATIC)
+EXECUTABLES = field
+
+all : $(EXECUTABLES)
+
+lib : $(wildcard $(WORK_DIR)/common/*.F90)
+	make -C $(WORK_DIR)/common
+
+$(EXECUTABLES) : % : %.o $(wildcard $(COMMON_DIR)/*.F90)
+	make lib
+	$(FC) $(LDFLAGS) -o $@ $< $(LIB)
+
+clean :
+	-rm -f $(EXECUTABLES) *.o *.mod *.a core
+
+# Dependency
+.depend : $(wildcard *.F90)
+	$(MAKEDEPEND_BIN) $^ > .depend
+
+include .depend
diff --git a/examples/velocity_field/README.md b/examples/velocity_field/README.md
new file mode 100644
index 0000000..8fc9a62
--- /dev/null
+++ b/examples/velocity_field/README.md
@@ -0,0 +1,30 @@
+# Velocity Fields for RBC3D
+
+While other provided testcases simulate blood cells in a wall geometry, this example will calculate and output the velocities in a grid within an already-completed simulation.
+
+## What Does This Do?
+
+After building with `make`, the `field` program will generate a list of coordinates as a 3D grid inside the simulation's bounding box (without any regard to the location of the wall or the cells).
+Then, it calculates velocities (using the same method to calculate blood cell motion) at each point in the list.
+
+### How To Run?
+
+After running your blood flow simulation, you should have a series of restart files named in the format: `restart{timestep}.dat` corresponding to the state of the simulation at the timestep, and `restart.LATEST.dat` is the most recent restart file.
+Pick the restart file corresponding to the timestep for which you want the velocity field, and insert into the `./D/` directory.
+Then, modify the `./Input/tube.in` configuration file to match the restart file's path; the current path is `D/restart.LATEST.dat`.
+
+
+### Output File
+
+After runnning the simulation, you will receive an output file `./D/field.csv`. This file contains the following comma-separated columns:
+
+X | Y | Z | Vx | Vy | Vz
+
+where, for each row, (X, Y, Z) denotes the location, and (Vx, Vy, Vz) is the velocity at that location.
+
+#### Visualization in ParaView
+
+You can easily visualize this velocity field in ParaView with the following steps after Importing and Parsing `field.csv`:
+1. Apply the "Table to Points" Filter, selecting X=X, Y=Y, and Z=Z as your coordinates.
+2. Apply the "Merge Vector Components" Filter, selecting X=Vx, Y=Vy, and Z=Vz for your vectors. Rename this to "velocity."
+3. Apply the "Glyph" filter, selecting the orientation to be "velocity." Scale and color as necessary. 
\ No newline at end of file
diff --git a/examples/case/field.F90 b/examples/velocity_field/field.F90
similarity index 71%
rename from examples/case/field.F90
rename to examples/velocity_field/field.F90
index a230e7c..f050df6 100644
--- a/examples/case/field.F90
+++ b/examples/velocity_field/field.F90
@@ -38,23 +38,15 @@ SUBROUTINE GetVelocity
     type(t_RBC), pointer :: rbc
     type(t_TargetList) :: tlist
     real(WP), allocatable :: xs(:, :), vs(:, :)
-    real(WP), allocatable :: x(:, :, :), v(:, :, :)
-    real(WP) :: xmin(3), xmax(3), szcell(3)
-    integer ::  N1, N2, N3, npoint, i1, i2, i3, idx
+    integer :: npoint
     character(CHRLEN) :: fn
     integer, parameter :: funit = 90
     integer :: ierr
 
-    real(WP) :: Xo
-    integer  :: Ny, Nz
+    integer  :: Nx, Ny, Nz
     real(WP), parameter :: dx0 = 0.3
     real(WP), dimension(3) :: vl
 
-    !======================================================================
-    ! Set up background velocity and RBC surface density
-    !    vBkg(1) = 0.
-    !    vBkg(2) = 0.
-    !    vBkg(3) = 2.5
     call MPI_Bcast(vBkg, 3, MPI_WP, 0, MPI_COMM_WORLD, ierr); 
     do irbc = 1, nrbc
       rbc => rbcs(irbc)
@@ -87,13 +79,11 @@ SUBROUTINE GetVelocity
       end if
     end if
 
-    Xo = Lb(1)/2.
-    Ny = NINT(Lb(2)/dx0); print *, "Ny = ", Ny
-    Nz = NINT(Lb(3)/dx0); print *, "Nz = ", Nz
-    allocate (x(Ny, Nz, 3))
+    Nx = NINT(Lb(1)/dx0)
+    Ny = NINT(Lb(2)/dx0)
+    Nz = NINT(Lb(3)/dx0)
 
-    call MakeMesh(Xo, Ny, Nz, x)
-    npoint = Ny*Nz
+    npoint = Nx*Ny*Nz
     allocate (xs(npoint, 3), vs(npoint, 3))
 
     ! Evolve cells
@@ -102,9 +92,7 @@ SUBROUTINE GetVelocity
     ! Enforce no-slip condition on the wall
     call NoSlipWall
 
-    call leukVel(rbcs(1), vl)
-
-    xs = RESHAPE(x, SHAPE(xs))
+    call MakeMesh(Nx, Ny, Nz, xs)
 
     print *, "A"
     call TargetList_CreateFromRaw(tlist, xs)
@@ -112,74 +100,53 @@ SUBROUTINE GetVelocity
     print *, "B"
     call CalcVelocityField(tlist, vs)
 
-    allocate (v(Ny, Nz, 3))
-
-    v = RESHAPE(vs, SHAPE(v))
-    v(:, :, 3) = v(:, :, 3) - vBkg(3)
-    !    v(:,:,:,3) = v(:,:,:,3) - vl(3)  ! substract leuk velocity
-
-    call writeVel(Ny, Nz, v, x)
+    call WriteVelCsv(npoint, vs, xs)
 
     call TargetList_Destroy(tlist)
     deallocate (xs, vs)
 
   END SUBROUTINE GetVelocity
 
-  SUBROUTINE writeVel(Ny, Nz, v, x)
-    integer                          :: Ny, Nz
-    real(WP), dimension(Ny, Nz, 3)     :: v, x
-
-    integer   :: i, j, k, m
-
-    open (1, file='v.f')
-    write (1, *) Ny, Nz, 3
-    write (1, *) v
-    close (1)
+  subroutine WriteVelCsv(npoint, vels, pts)
+    integer :: npoint
+    real(WP), dimension(npoint, 3) :: vels, pts
 
-    open (1, file='v.g')
-    write (1, *) Ny, Nz
-    write (1, *) x(:, :, 2:3)
-    close (1)
+    integer :: i
 
-  END SUBROUTINE writeVel
-
-  SUBROUTINE MakeMesh(Xo, Ny, Nz, x)
-    real(WP)                         :: Xo
-    integer                          :: Ny, Nz
-    real(WP), dimension(:, :, :)      :: x
-
-    real(WP)  :: y, z, dy, dz
-    integer   :: i, j, k
-
-    dy = Lb(2)/REAL(Ny); print *, "dy = ", dy
-    dz = Lb(3)/REAL(Nz); print *, "dz = ", dz
-
-    do k = 1, Nz
-      z = 0.5*dz + REAL(k - 1)*dz
-      do j = 1, Ny
-        y = 0.5*dy + REAL(j - 1)*dy
-        x(j, k, 1) = Lb(1)/2
-        x(j, k, 2) = y
-        x(j, k, 3) = z
-        !  do k = 1,Nz
-        !     z = 0.5*dz+REAL(k-1)*dz
-        !     do j = 1,Ny
-        !        y = 0.5*dy+REAL(j-1)*dy
-        !       !  x(j,k,1) = Lb(1)/2.
-        !       !  x(j,k,2) = y
-        !        x(j,k,1) = x
-        !        x(j,k,3) = z
-
-        !  print *, "COORDINATE", x(j, k, 1), x(j, k, 2), x(j, k, 3)
-      end do
-    end do
-
-    open (1, file='v.g')
-    write (1, *) Ny, Nz
-    write (1, *) x(:, :, 2:3)
+    fn = "./D/field.csv"
+    open (1, file=fn)
+    write (1, *) "X,Y,Z,Vx,Vy,Vz"
+    do i = 1, npoint
+      write (1, *) pts(i, 1), ",", pts(i, 2), ",", pts(i, 3), ",", vels(i, 1), ",", vels(i, 2), ",", vels(i, 3)
+    end do !i
     close (1)
-
-  END SUBROUTINE MakeMesh
+  end subroutine writeVelCsv
+
+  subroutine MakeMesh(Nx, Ny, Nz, pts)
+    integer :: Nx, Ny, Nz
+    real(WP), dimension(:, :) :: pts ! has dimensions (Nx * Ny * Nz , 3)
+
+    real(WP) :: x, y, z, dx, dy, dz
+    integer :: x_it, y_it, z_it, p_it
+
+    dx = Lb(1)/REAL(Nx); 
+    dy = Lb(2)/REAL(Ny); 
+    dz = Lb(3)/REAL(Nz); 
+    p_it = 1
+    do x_it = 1, Nx
+      x = 0.5*dx + REAL(x_it - 1)*dx
+      do y_it = 1, Ny
+        y = 0.5*dy + REAL(y_it - 1)*dy
+        do z_it = 1, Nz
+          z = 0.5*dz + REAL(z_it - 1)*dz
+
+          pts(p_it, :) = (/x, y, z/)
+          p_it = p_it + 1
+
+        end do !z_it
+      end do !y_it
+    end do !x_it
+  end subroutine MakeMesh
 
   SUBROUTINE leukVel(rbc, vl)
     type(t_RBC)            :: rbc
diff --git a/sample_files/sample_walls/README.md b/sample_files/sample_walls/README.md
index 180034c..08c9556 100644
--- a/sample_files/sample_walls/README.md
+++ b/sample_files/sample_walls/README.md
@@ -18,6 +18,17 @@ This wall is "U"-shaped, curving downward and back upward. The wall has a z-leng
 The mesh a very simplified blood vessel 'valve', as a cylinder with 2 'flaps' protruding into the center. The cylinder is 13.33 units long, and has a diameter of 6 units.
 
 
+## Carotid Artery and Web Meshes
+We include a series of meshes for simulating the geometry of the carotid artery. These meshes can be combined by including multiple walls in simulations. A Carotid Web example case is provided  in `examples\carotid_web\`.
+
+### `carotid.e`
+This is a simplified model of the Carotid bulb area, with narrow entry and exit points, and a larger center bulge.
+
+### `web.e` and `vertical_web.e`
+These two meshes are models which approximate the "web" in the carotid artery, a condition in which fibrous tissue builds up inside the bulb.
+They can be loaded as a second wall along with the `carotid.e` geometry to simulate a carotid artery with a web.
+
+
 ## Creating a Wall Mesh
 
 To create your own wall-mesh for RBC3D, simply use the Coreform Cubit software to create a tri-mesh, and export the mesh. This will create an Exodus-II formatted file, which can be imported into your simulation using:
diff --git a/sample_files/sample_walls/wall_journals/carotid_and_web.e b/sample_files/sample_walls/wall_journals/carotid_and_web.e
new file mode 100644
index 0000000..19bae64
Binary files /dev/null and b/sample_files/sample_walls/wall_journals/carotid_and_web.e differ
diff --git a/sample_files/sample_walls/wall_journals/vertical_web.e b/sample_files/sample_walls/wall_journals/vertical_web.e
new file mode 100644
index 0000000..47e5593
Binary files /dev/null and b/sample_files/sample_walls/wall_journals/vertical_web.e differ
diff --git a/sample_files/sample_walls/wall_meshes/carotid.e b/sample_files/sample_walls/wall_meshes/carotid.e
new file mode 100644
index 0000000..5fe6cba
Binary files /dev/null and b/sample_files/sample_walls/wall_meshes/carotid.e differ
diff --git a/sample_files/sample_walls/wall_meshes/vertical_web.e b/sample_files/sample_walls/wall_meshes/vertical_web.e
new file mode 100644
index 0000000..47e5593
Binary files /dev/null and b/sample_files/sample_walls/wall_meshes/vertical_web.e differ
diff --git a/sample_files/sample_walls/wall_meshes/web.e b/sample_files/sample_walls/wall_meshes/web.e
new file mode 100644
index 0000000..a72fd29
Binary files /dev/null and b/sample_files/sample_walls/wall_meshes/web.e differ