New Base.axes for AbstractField (#3801)

* Use same Base.axes for all AbstractField

* Update abstract_field.jl

* Do better for FieldtimeSeries

src/AbstractOperations/AbstractOperations.jl

@@ -34,15 +34,6 @@ abstract type AbstractOperation{LX, LY, LZ, G, T} <: AbstractField{LX, LY, LZ, G
 
 const AF = AbstractField # used in unary_operations.jl, binary_operations.jl, etc
 
-function Base.axes(f::AbstractOperation)
-    idx = indices(f)
-    if idx === (:, : ,:)
-        return Base.OneTo.(size(f))
-    else
-        return Tuple(idx[i] isa Colon ? Base.OneTo(size(f, i)) : idx[i] for i = 1:3)
-    end
-end
-
 # We have no halos to fill
 @inline fill_halo_regions!(::AbstractOperation, args...; kwargs...) = nothing
 

src/Fields/abstract_field.jl

@@ -58,6 +58,39 @@ Base.size(f::AbstractField) = size(f.grid, location(f))
 Base.length(f::AbstractField) = prod(size(f))
 Base.parent(f::AbstractField) = f
 
+const Abstract3DField = AbstractField{<:Any, <:Any, <:Any, <:Any, <:Any, 3}
+const Abstract4DField = AbstractField{<:Any, <:Any, <:Any, <:Any, <:Any, 4}
+
+# TODO: to omit boundaries on Face fields, we have to return 2:N
+# when topo=Bounded, and loc=Face
+@inline axis(::Colon, N) = Base.OneTo(N)
+@inline axis(index::UnitRange, N) = index
+
+@inline function Base.axes(f::Abstract3DField)
+    Nx, Ny, Nz = size(f)
+    ix, iy, iz = indices(f)
+
+    ax = axis(ix, Nx)
+    ay = axis(iy, Ny)
+    az = axis(iz, Nz)
+
+    return (ax, ay, az)
+end
+
+@inline function Base.axes(f::Abstract4DField)
+    Nx, Ny, Nz, Nt = size(f)
+    ix, iy, iz = indices(f)
+
+    ax = axis(ix, Nx)
+    ay = axis(iy, Ny)
+    az = axis(iz, Nz)
+    at = Base.OneTo(Nt)
+
+    return (ax, ay, az, at)
+end
+
+
+
 """
     total_size(field::AbstractField)
 

src/Fields/field.jl

@@ -398,20 +398,6 @@ interior(f::Field, I...) = view(interior(f), I...)
 # Don't use axes(f) to checkbounds; use axes(f.data)
 Base.checkbounds(f::Field, I...) = Base.checkbounds(f.data, I...)
 
-@inline axis(::Colon, N) = Base.OneTo(N)
-@inline axis(index::UnitRange, N) = index
-
-@inline function Base.axes(f::Field)
-    Nx, Ny, Nz = size(f)
-    ix, iy, iz = f.indices
-
-    ax = axis(ix, Nx)
-    ay = axis(iy, Ny)
-    az = axis(iz, Nz)
-
-    return (ax, ay, az)
-end
-
 @propagate_inbounds Base.getindex(f::Field, inds...) = getindex(f.data, inds...)
 @propagate_inbounds Base.getindex(f::Field, i::Int)  = parent(f)[i]
 @propagate_inbounds Base.setindex!(f::Field, val, i, j, k) = setindex!(f.data, val, i, j, k)