fix connecting periodic boundaries for tree3d; fix place_boxes algorithm

src/tree.jl

@@ -439,23 +439,39 @@ function _connect_periodic!(t::Tree3d, ax::Cart3d.Axis)
   dirmin, dirmax = directions(ax)
   ortho_ax1, ortho_ax2 = orthogonal(ax)
 
-  # filter boundary cells
+  # filter boundary cells, returns a list of indices for the t.cells array
   i_bdrycells = findall(t.cells) do cell
     cell.neighbors[dirmin] == 0 || cell.neighbors[dirmax] == 0
   end
+  length(i_bdrycells) == 0 && error("Could not find any unconnected boundaries")
+  # determine those is for which a dirmin neighbor is missing
+  is_of_i_bdrycells_dirmin = findall(i_bdrycells) do i
+    cell = t.cells[i]
+    cell.neighbors[dirmin] == 0
+  end
+  i_bdrycells_dirmin = i_bdrycells[is_of_i_bdrycells_dirmin]
+  # demand that the number of boundary cells connected between the periodic sides is the same
+  2*length(i_bdrycells_dirmin) != length(i_bdrycells) &&
+    error("Encountered an unbalanced number of boundaries to connect")
+  i_bdrycells_dirmax = setdiff(i_bdrycells, i_bdrycells_dirmin)
+
 
   # connect pairs of cells with same coordinates in orthogonal directions
-  for i in i_bdrycells
-    cell = t.cells[i]
-    coords = t.coords[i]
-    j = findfirst(i_bdrycells) do jj
-      jj != i && coords[ortho_ax1] == t.coords[jj][ortho_ax1] && 
-                 coords[ortho_ax2] == t.coords[jj][ortho_ax2] && return true
+  # only iterating the ones which miss a dirmin neighbor
+  for i_dirmin in i_bdrycells_dirmin
+    cell = t.cells[i_dirmin]
+    coords = t.coords[i_dirmin]
+    j = findfirst(i_bdrycells_dirmax) do jj
+      coords[ortho_ax1] == t.coords[jj][ortho_ax1] &&
+      coords[ortho_ax2] == t.coords[jj][ortho_ax2] && return true
+      return false
+      coords[ortho_ax] == t.coords[jj][ortho_ax] && return true
       return false
     end
-    isnothing(j) && error("Could not find opposing cell to '$cell'")
-    cell.neighbors[dirmin]       = j
-    t.cells[j].neighbors[dirmax] = i
+    isnothing(j) && error("Could not find opposing cell for '$cell'")
+    i_dirmax = i_bdrycells_dirmax[j]
+    cell.neighbors[dirmin]              = i_dirmax
+    t.cells[i_dirmax].neighbors[dirmax] = i_dirmin
   end
 
 end
@@ -596,9 +612,9 @@ function place_boxes(t::Tree2d, x_extends, y_extends)
     prev_cell, cell = t.cells[i-1], t.cells[i]
     !is_neighbor(prev_cell, cell) && error_failed_placement(cell, prev_cell)
 
-    prev_coords, coords = t.coords[i-1], t.coords[1]
+    prev_coords, coords = t.coords[i-1], t.coords[i]
     touch_at = touch_direction(prev_coords, coords)
-    prev_box = boxes[i-1]
+    prev_box = last(boxes)
     (pxmin, pxmax), (pymin, pymax) = prev_box.extends
     xmin, xmax, ymin, ymax = if touch_at === Cart2d.Xmax
       pxmax, pxmax+dx, pymin, pymax
@@ -634,9 +650,9 @@ function place_boxes(t::Tree3d, x_extends, y_extends, z_extends)
     prev_cell, cell = t.cells[i-1], t.cells[i]
     !is_neighbor(prev_cell, cell) && error_failed_placement(cell, prev_cell)
 
-    prev_coords, coords = t.coords[i-1], t.coords[1]
+    prev_coords, coords = t.coords[i-1], t.coords[i]
     touch_at = touch_direction(prev_coords, coords)
-    prev_box = boxes[i-1]
+    prev_box = last(boxes)
     (pxmin, pxmax), (pymin, pymax), (pzmin, pzmax) = prev_box.extends
     xmin, xmax, ymin, ymax, zmin, zmax = if touch_at === Cart3d.Xmax
       pxmax, pxmax+dx, pymin, pymax, pzmin, pzmax

test/dg1dTests/src/test_tree.jl

@@ -187,13 +187,16 @@ end
 
 @testset "Tree" begin
 
-  t1d = Tree1d(5)
-  t2d = Tree2d(5,6)
-  t3d = Tree3d(5,6,7)
+  Nx, Ny, Nz = 5, 6, 7
 
-  @test n_cells(t1d) == 5
-  @test n_cells(t2d) == 5*6
-  @test n_cells(t3d) == 5*6*7
+  # test with periodicity here, because that requires an extra step for setup
+  t1d = Tree1d(Nx, periodic=(true,))
+  t2d = Tree2d(Nx,Ny, periodic=(true,true))
+  t3d = Tree3d(Nx,Ny,Nz, periodic=(true,true,true))
+
+  @test n_cells(t1d) == Nx
+  @test n_cells(t2d) == Nx*Ny
+  @test n_cells(t3d) == Nx*Ny*Nz
 
   # t1d = Tree1d(5, periodic=(true,))
   # t2d = Tree2d(5,6, periodic=(true,false))