Mercurial > repos > public > sbplib_julia

comparison EquidistantGrid.jl @ 134:79699dda29be

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in cell_based_test
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 21 Feb 2019 16:27:28 +0100
parents 155bbecf18bb
children 99308f68e548
comparison
equal deleted inserted replaced
84:48079bd39969 134:79699dda29be
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 inverse_spacing::NTuple{Dim, T} # The reciprocal of the grid spacing
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 inverse_spacing = (size.-1)./abs.(limit_upper.-limit_lower)
18 return new{Dim,T}(size, limit_lower, limit_upper, inverse_spacing)
18 end 19 end
20 end
19 21
20 # # 1D constructor which can be called as EquidistantGrid(m, (xl,xr)) 22 function Base.eachindex(grid::EquidistantGrid)
21 # function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real}) 23 CartesianIndices(grid.size)
22 # return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
23 # end
24
25 end 24 end
26 25
27 # Returns the number of dimensions of an EquidistantGrid. 26 # Returns the number of dimensions of an EquidistantGrid.
28 # 27 #
29 # @Input: grid - an EquidistantGrid 28 # @Input: grid - an EquidistantGrid
30 # @Return: numberOfPoints - The number of dimensions 29 # @Return: dimension - The dimension of the grid
31 function numberOfDimensions(grid::EquidistantGrid) 30 function dimension(grid::EquidistantGrid)
32 return length(grid.numberOfPointsPerDim) 31 return length(grid.size)
33 end 32 end
34 33
35 # Computes the total number of points of an EquidistantGrid. 34 # Returns the spacing of the grid
36 # 35 #
37 # @Input: grid - an EquidistantGrid 36 function spacing(grid::EquidistantGrid)
38 # @Return: numberOfPoints - The total number of points 37 return 1.0./grid.inverse_spacing
39 function numberOfPoints(grid::EquidistantGrid)
40 return prod(grid.numberOfPointsPerDim)
41 end 38 end
42 39
43 # Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance 40 # Computes the points of an EquidistantGrid as an array of tuples with
44 # between two points for each coordinate direction. 41 # the same dimension as the grid.
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
51
52 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered
53 # such that points in the first coordinate direction varies first, then the second
54 # and lastely the third (if applicable)
55 # 42 #
56 # @Input: grid - an EquidistantGrid 43 # @Input: grid - an EquidistantGrid
57 # @Return: points - the points of the grid. 44 # @Return: points - the points of the grid.
58 function points(grid::EquidistantGrid) 45 function points(grid::EquidistantGrid)
59 dx̄ = (grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1) 46 # TODO: Make this return an abstract array?
60 47 indices = Tuple.(CartesianIndices(grid.size))
61 points = Vector{typeof(dx̄)}(undef, numberOfPoints(grid)) 48 h = spacing(grid)
62 # Compute the points based on their Cartesian indices and the signed 49 return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices)
63 # grid spacings
64 cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
65 for i ∈ 1:numberOfPoints(grid)
66 ci = Tuple(cartesianIndices[i]) .-1
67 points[i] = grid.limit_lower .+ dx̄.*ci
68 end
69
70 # TBD: Keep? this? How do we want to represent points in 1D?
71 if numberOfDimensions(grid) == 1
72 points = broadcast(x -> x[1], points)
73 end
74 return points
75 end 50 end
76 51
77 function pointsalongdim(grid::EquidistantGrid, dim::Integer) 52 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
78 @assert dim<=numberOfDimensions(grid) 53 @assert dim<=dimension(grid)
79 @assert dim>0 54 @assert dim>0
80 points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.numberOfPointsPerDim[dim]) 55 points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.size[dim])
81 end 56 end
82
83 using PyPlot, PyCall
84
85 function plotgridfunction(grid::EquidistantGrid, gridfunction)
86 if numberOfDimensions(grid) == 1
87 plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
88 elseif numberOfDimensions(grid) == 2
89 mx = grid.numberOfPointsPerDim[1]
90 my = grid.numberOfPointsPerDim[2]
91 X = repeat(pointsalongdim(grid,1),1,my)
92 Y = permutedims(repeat(pointsalongdim(grid,2),1,mx))
93 plot_surface(X,Y,reshape(gridfunction,mx,my));
94 else
95 error(string("Plot not implemented for dimension ", string(numberOfDimensions(grid))))
96 end
97 end