Mercurial > repos > public > sbplib_julia

comparison EquidistantGrid.jl @ 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 9d53ecca34f7
children 1aaeb46ba5f4
comparison
equal deleted inserted replaced
122:6c6979ff17f4 124:631eb9b35d72
3 # by the exterior product of the vectors obtained by projecting (x̄₂-x̄₁) onto 3 # by the exterior product of the vectors obtained by projecting (x̄₂-x̄₁) onto
4 # the coordinate directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2) 4 # the coordinate directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2)
5 # the domain is defined as (-1,1)x(0,2). 5 # the domain is defined as (-1,1)x(0,2).
6 6
7 struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid 7 struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
8 numberOfPointsPerDim::NTuple{Dim, Int} # First coordinate direction stored first, then 8 size::NTuple{Dim, Int} # First coordinate direction stored first
9
10 limit_lower::NTuple{Dim, T} 9 limit_lower::NTuple{Dim, T}
11 limit_upper::NTuple{Dim, T} 10 limit_upper::NTuple{Dim, T}
11 spacing::NTuple{Dim, T}
12 12
13 # General constructor 13 # General constructor
14 function EquidistantGrid(nPointsPerDim::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T 14 function EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
15 @assert all(nPointsPerDim.>0) 15 @assert all(size.>0)
16 @assert all(limit_upper.-limit_lower .!= 0) 16 @assert all(limit_upper.-limit_lower .!= 0)
17 return new{Dim,T}(nPointsPerDim, limit_lower, limit_upper) 17 spacing = abs.(limit_upper.-limit_lower)./(size.-1)
18 return new{Dim,T}(size, limit_lower, limit_upper, spacing)
18 end 19 end
19
20 # # 1D constructor which can be called as EquidistantGrid(m, (xl,xr))
21 # function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real})
22 # return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
23 # end
24
25 end 20 end
26 21
27 # Returns the number of dimensions of an EquidistantGrid. 22 # Returns the number of dimensions of an EquidistantGrid.
28 # 23 #
29 # @Input: grid - an EquidistantGrid 24 # @Input: grid - an EquidistantGrid
30 # @Return: numberOfPoints - The number of dimensions 25 # @Return: dimension - The dimension of the grid
31 function numberOfDimensions(grid::EquidistantGrid) 26 function dimension(grid::EquidistantGrid)
32 return length(grid.numberOfPointsPerDim) 27 return length(grid.size)
33 end
34
35 # Computes the total number of points of an EquidistantGrid.
36 #
37 # @Input: grid - an EquidistantGrid
38 # @Return: numberOfPoints - The total number of points
39 function numberOfPoints(grid::EquidistantGrid)
40 return prod(grid.numberOfPointsPerDim)
41 end
42
43 # Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance
44 # between two points for each coordinate direction.
45 #
46 # @Input: grid - an EquidistantGrid
47 # @Return: h̄ - Grid spacing for each coordinate direction stored in a tuple.
48 function spacings(grid::EquidistantGrid)
49 return abs.(grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
50 end 28 end
51 29
52 function Base.eachindex(grid::EquidistantGrid) 30 function Base.eachindex(grid::EquidistantGrid)
53 CartesianIndices(grid.numberOfPointsPerDim) 31 CartesianIndices(grid.size)
54 end 32 end
55 33
56 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered 34 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered
57 # such that points in the first coordinate direction varies first, then the second 35 # such that points in the first coordinate direction varies first, then the second
58 # and lastely the third (if applicable) 36 # and lastely the third (if applicable)
60 # @Input: grid - an EquidistantGrid 38 # @Input: grid - an EquidistantGrid
61 # @Return: points - the points of the grid. 39 # @Return: points - the points of the grid.
62 function points(grid::EquidistantGrid) 40 function points(grid::EquidistantGrid)
63 # TODO: Make this return an abstract array? 41 # TODO: Make this return an abstract array?
64 physical_domain_size = (grid.limit_upper .- grid.limit_lower) 42 physical_domain_size = (grid.limit_upper .- grid.limit_lower)
65 indices = Tuple.(CartesianIndices(grid.numberOfPointsPerDim)) 43 indices = Tuple.(CartesianIndices(grid.size))
66 return broadcast(I -> grid.limit_lower .+ physical_domain_size.*(I.-1), indices) 44 return broadcast(I -> grid.limit_lower .+ physical_domain_size.*(I.-1), indices)
67 end 45 end
68 46
69 function pointsalongdim(grid::EquidistantGrid, dim::Integer) 47 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
70 @assert dim<=numberOfDimensions(grid) 48 @assert dim<=dimension(grid)
71 @assert dim>0 49 @assert dim>0
72 points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.numberOfPointsPerDim[dim]) 50 points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.size[dim])
73 end 51 end
74
75 using PyPlot, PyCall
76
77 function plotgridfunction(grid::EquidistantGrid, gridfunction)
78 if numberOfDimensions(grid) == 1
79 plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
80 elseif numberOfDimensions(grid) == 2
81 mx = grid.numberOfPointsPerDim[1]
82 my = grid.numberOfPointsPerDim[2]
83 X = repeat(pointsalongdim(grid,1),1,my)
84 Y = permutedims(repeat(pointsalongdim(grid,2),1,mx))
85 plot_surface(X,Y,reshape(gridfunction,mx,my));
86 else
87 error(string("Plot not implemented for dimension ", string(numberOfDimensions(grid))))
88 end
89 end