add bulkidxs, bdryidxs to Mesh structs

src/mesh.jl

@@ -13,6 +13,8 @@ struct Mesh{N_Dim,N_Sides} <: AbstractMesh
   element::SpectralElement
   cache::Cache
   offsets::Vector{Int64} # data index offsets for each cell
+  bulkidxs::Vector{Int64}
+  bdryidxs::NTuple{N_Sides,Vector{Int64}}
   # TODO Move to evolution part!
   CFL::Float64
 end
@@ -64,7 +66,10 @@ function Mesh1d(; N=5, K=10, range=[-1.0,1.0], CFL=0.5, basis="lgl_lumped", peri
     nx[2*(i-1)+2] = 1.0
   end
 
-  return Mesh1d(tree, boxes, extends, element, cache, offsets, CFL)
+  bulkidxs = collect(2:K*Npts-1)
+  bxmin, bxmax = 1, K*Npts
+
+  return Mesh1d(tree, boxes, extends, element, cache, offsets, bulkidxs, (bxmin, bxmax), CFL)
 end
 
 
@@ -108,6 +113,8 @@ function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
   # outward facing normals on rectangular grid
   @unpack nx, ny = get_bdry_variables(cache)
 
+  bulkidxs = Int64[]
+  bxmin, bxmax, bymin, bymax = Int64[], Int64[], Int64[], Int64[]
   for (i,b) in enumerate(boxes)
     (xmin, xmax), (ymin,ymax) = extrema(b)
     xi = @. 0.5 * ((xmax - xmin) * element_x.z + (xmax + xmin))
@@ -136,9 +143,32 @@ function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
     @views ny[rng_down] .= -1.0
     @views nx[rng_up]   .= 0.0
     @views ny[rng_up]   .= 1.0
+
+    cell = tree.cells[i]
+    isbdrycell = false
+    if !has_neighbor(cell, Cart2d.Xmin)
+      isbdrycell = true
+      push!(bxmin, cell.index)
+    end
+    if !has_neighbor(cell, Cart2d.Xmax)
+      isbdrycell = true
+      push!(bxmax, cell.index)
+    end
+    if !has_neighbor(cell, Cart2d.Ymin)
+      isbdrycell = true
+      push!(bymin, cell.index)
+    end
+    if !has_neighbor(cell, Cart2d.Ymax)
+      isbdrycell = true
+      push!(bymax, cell.index)
+    end
+    if !isbdrycell
+      push!(bulkidxs, cell.index)
+    end
   end
 
-  return Mesh2d(tree, boxes, extends, element_x, cache, offsets, CFL)
+  return Mesh2d(tree, boxes, extends, element_x, cache, offsets,
+                bulkidxs, (bxmin, bxmax, bymin, bymax), CFL)
 end
 
 

src/with_signature.jl

@@ -585,8 +585,9 @@ function _broadcast_face_2!(f, mesh::Mesh2d)
   Kx, Ky = mesh.tree.dims
   Nx, Ny = Npts, Npts
 
-  for (icell,cell_in) = enumerate(mesh.tree.cells)
+  @inbounds for icell in mesh.bulkidxs
 
+    cell_in = mesh.tree.cells[icell]
     bulk_in = cellindices(mesh, icell)
     rng_lhs, rng_rhs, rng_down, rng_up = faceindices(mesh, icell)