Mercurial > repos > public > sbplib_julia

changeset 124:631eb9b35d72 cell_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make grid spacing a property of EquidistantGrid. Create plotting module for sbplib.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 13 Feb 2019 10:37:52 +0100
parents 6c6979ff17f4
children 22642722a8ec
files AbstractGrid.jl EquidistantGrid.jl diffOp.jl sbpPlot.jl
diffstat 4 files changed, 36 insertions(+), 68 deletions(-) [+]
line wrap: on
line diff
--- a/AbstractGrid.jl	Tue Feb 12 15:18:18 2019 +0100
+++ b/AbstractGrid.jl	Wed Feb 13 10:37:52 2019 +0100
@@ -1,10 +1,6 @@
 abstract type AbstractGrid end
 
-function numberOfDimensions(grid::AbstractGrid)
-    error("Not implemented for abstact type AbstractGrid")
-end
-
-function numberOfPoints(grid::AbstractGrid)
+function dimension(grid::AbstractGrid)
     error("Not implemented for abstact type AbstractGrid")
 end
 
--- a/EquidistantGrid.jl	Tue Feb 12 15:18:18 2019 +0100
+++ b/EquidistantGrid.jl	Wed Feb 13 10:37:52 2019 +0100
@@ -5,52 +5,30 @@
 # the domain is defined as (-1,1)x(0,2).
 
 struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
-    numberOfPointsPerDim::NTuple{Dim, Int} # First coordinate direction stored first, then
-
+    size::NTuple{Dim, Int} # First coordinate direction stored first
     limit_lower::NTuple{Dim, T}
     limit_upper::NTuple{Dim, T}
+    spacing::NTuple{Dim, T}
 
     # General constructor
-    function EquidistantGrid(nPointsPerDim::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
-        @assert all(nPointsPerDim.>0)
+    function EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
+        @assert all(size.>0)
         @assert all(limit_upper.-limit_lower .!= 0)
-        return new{Dim,T}(nPointsPerDim, limit_lower, limit_upper)
+        spacing = abs.(limit_upper.-limit_lower)./(size.-1)
+        return new{Dim,T}(size, limit_lower, limit_upper, spacing)
     end
-
-    # # 1D constructor which can be called as EquidistantGrid(m, (xl,xr))
-    # function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real})
-    #     return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
-    # end
-
 end
 
 # Returns the number of dimensions of an EquidistantGrid.
 #
 # @Input: grid - an EquidistantGrid
-# @Return: numberOfPoints - The number of dimensions
-function numberOfDimensions(grid::EquidistantGrid)
-    return length(grid.numberOfPointsPerDim)
-end
-
-# Computes the total number of points of an EquidistantGrid.
-#
-# @Input: grid - an EquidistantGrid
-# @Return: numberOfPoints - The total number of points
-function numberOfPoints(grid::EquidistantGrid)
-    return prod(grid.numberOfPointsPerDim)
-end
-
-# Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance
-# between two points for each coordinate direction.
-#
-# @Input: grid - an EquidistantGrid
-# @Return: h̄ - Grid spacing for each coordinate direction stored in a tuple.
-function spacings(grid::EquidistantGrid)
-    return abs.(grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
+# @Return: dimension - The dimension of the grid
+function dimension(grid::EquidistantGrid)
+    return length(grid.size)
 end
 
 function Base.eachindex(grid::EquidistantGrid)
-    CartesianIndices(grid.numberOfPointsPerDim)
+    CartesianIndices(grid.size)
 end
 
 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered
@@ -62,28 +40,12 @@
 function points(grid::EquidistantGrid)
     # TODO: Make this return an abstract array?
     physical_domain_size = (grid.limit_upper .- grid.limit_lower)
-    indices = Tuple.(CartesianIndices(grid.numberOfPointsPerDim))
+    indices = Tuple.(CartesianIndices(grid.size))
     return broadcast(I -> grid.limit_lower .+ physical_domain_size.*(I.-1), indices)
 end
 
 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
-    @assert dim<=numberOfDimensions(grid)
+    @assert dim<=dimension(grid)
     @assert dim>0
-    points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.numberOfPointsPerDim[dim])
+    points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.size[dim])
 end
-
-using PyPlot, PyCall
-
-function plotgridfunction(grid::EquidistantGrid, gridfunction)
-    if numberOfDimensions(grid) == 1
-        plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
-    elseif numberOfDimensions(grid) == 2
-        mx = grid.numberOfPointsPerDim[1]
-        my = grid.numberOfPointsPerDim[2]
-        X = repeat(pointsalongdim(grid,1),1,my)
-        Y = permutedims(repeat(pointsalongdim(grid,2),1,mx))
-        plot_surface(X,Y,reshape(gridfunction,mx,my));
-    else
-        error(string("Plot not implemented for dimension ", string(numberOfDimensions(grid))))
-    end
-end
--- a/diffOp.jl	Tue Feb 12 15:18:18 2019 +0100
+++ b/diffOp.jl	Wed Feb 13 10:37:52 2019 +0100
@@ -59,9 +59,7 @@
 
 # Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable
 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
-    N = D.grid.numberOfPointsPerDim
-    closuresize = closureSize(D.op)
-    for I ∈ regionindices(N, closuresize, (r1,r2))
+    for I ∈ regionindices(D.grid.size, closureSize(D.op), (r1,r2))
         @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2]))
         @inbounds u[I] = apply(D, v, indextuple)
     end
@@ -83,10 +81,7 @@
 
 using TiledIteration
 function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
-    N = D.grid.numberOfPointsPerDim
-    closuresize = closureSize(D.op)
-    ri = regionindices(N, closuresize, (r1,r2))
-
+    ri = regionindices(D.grid.size, closureSize(D.op), (r1,r2))
     for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2)) # TBD: Is this the right way, the right size?
         for j ∈ tileaxs[2], i ∈ tileaxs[1]
             I = ri[i,j]
@@ -114,19 +109,17 @@
 
 # u = L*v
 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
-    h = Grid.spacings(L.grid)[1]
-    uᵢ = L.a * apply(L.op, h, v, i)
+    uᵢ = L.a * apply(L.op, L.grid.spacing[1], v, i)
     return uᵢ
 end
 
 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
-    h = Grid.spacings(L.grid)
     # 2nd x-derivative
     @inbounds vx = view(v, :, Int(I[2]))
-    @inbounds uᵢ = L.a*apply(L.op, h[1], vx , I[1])
+    @inbounds uᵢ = L.a*apply(L.op, L.grid.spacing[1], vx , I[1])
     # 2nd y-derivative
     @inbounds vy = view(v, Int(I[1]), :)
-    @inbounds uᵢ += L.a*apply(L.op, h[2], vy, I[2])
+    @inbounds uᵢ += L.a*apply(L.op, L.grid.spacing[2], vy, I[2])
     return uᵢ
 end
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sbpPlot.jl	Wed Feb 13 10:37:52 2019 +0100
@@ -0,0 +1,17 @@
+module sbpPlot
+using PyPlot, PyCall
+
+function plotgridfunction(grid::EquidistantGrid, gridfunction)
+    if dimension(grid) == 1
+        plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
+    elseif dimension(grid) == 2
+        mx = grid.size[1]
+        my = grid.size[2]
+        X = repeat(pointsalongdim(grid,1),1,my)
+        Y = permutedims(repeat(pointsalongdim(grid,2),1,mx))
+        plot_surface(X,Y,reshape(gridfunction,mx,my));
+    else
+        error(string("Plot not implemented for dimension ", string(dimension(grid))))
+    end
+end
+end