Bugfix in Cartesian model with periodic BCs

This bug shows up in Julia 1.6 but not in 1.5

Author: fverdugo

src/Geometry/CartesianDiscreteModels.jl

@@ -408,7 +408,7 @@ function _cartesian_grid_topology_with_periodic_bcs(grid::UnstructuredGrid,
   isperiodic::NTuple,
   partition)
 
-  cell_to_vertices, vertex_to_node =
+  cell_to_vertices, vertex_to_node, =
     _generate_cell_to_vertices_from_grid(grid, isperiodic, partition)
   _generate_grid_topology_from_grid(grid,cell_to_vertices,vertex_to_node)
 end
@@ -418,8 +418,6 @@ function _generate_cell_to_vertices_from_grid(grid::UnstructuredGrid,
 
   if is_first_order(grid)
     nodes = get_cell_node_ids(grid)
-    cell_to_vertices = copy(nodes)
-
     nnodes = num_nodes(grid)
     num_nodes_x_dir = [partition[i]+1 for i in 1:length(partition)]
     point_to_isperiodic, slave_point_to_point, slave_point_to_master_point =
@@ -443,37 +441,31 @@ end
 function _generate_slave_to_master_point(num_nodes_x_dir::Vector{Int},
   isperiodic::NTuple, num_nodes::Int)
 
-  point_to_isperiodic = fill(false,num_nodes)
-
-  slave_ijk_bounds = Array{Any,1}(undef,length(isperiodic))
-  master_ijk_bounds = Array{Any,1}(undef,length(isperiodic))
-
-  linear_indices = LinearIndices(Tuple(num_nodes_x_dir))
   periodic_dirs = findall(x->x==true, isperiodic)
-  for periodic_dir in periodic_dirs
-    for i in 1:length(isperiodic)
-      if i == periodic_dir
-        slave_ijk_bounds[i] = num_nodes_x_dir[i]
-      else
-        slave_ijk_bounds[i] = 1:num_nodes_x_dir[i]
+  linear_indices = LinearIndices(Tuple(num_nodes_x_dir))
+  cartesian_indices = CartesianIndices(Tuple(num_nodes_x_dir))
+
+  point_to_isperiodic = fill(false,num_nodes)
+  for point in 1:length(point_to_isperiodic)
+    ci = Tuple(cartesian_indices[point])
+    for dir in periodic_dirs
+      if ci[dir] == num_nodes_x_dir[dir]
+        point_to_isperiodic[point] = true
       end
     end
-    slave_ijk_set = CartesianIndices(Tuple(slave_ijk_bounds))
-    point_to_isperiodic[linear_indices[slave_ijk_set]] .= true
   end
 
-  slave_point_to_point = findall( point_to_isperiodic)
+  slave_point_to_point = findall(point_to_isperiodic)
   slave_point_to_master_point = Array{Int32,1}(undef,length(slave_point_to_point))
 
-  cartesian_indices = CartesianIndices(Tuple(num_nodes_x_dir))
+  ijk = zeros(Int,length(isperiodic))
   for (i,point) in enumerate(slave_point_to_point)
-    ijk = collect(cartesian_indices[point].I)
+    ijk .= Tuple(cartesian_indices[point])
     for i in periodic_dirs
       if ijk[i] == num_nodes_x_dir[i]
         ijk[i] = 1
       end
     end
-
     master_point_ijk = CartesianIndex(Tuple(ijk))
     slave_point_to_master_point[i] = linear_indices[master_point_ijk]
   end