Mercurial > repos > public > sbplib_julia

comparison EquidistantGrid.jl @ 75:93c833019857

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make EquidistantGrid more concrete
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 17 Jan 2019 15:19:33 +0100
parents c47aaedcf71e
children 2be36b38389d
comparison
equal deleted inserted replaced
74:c47aaedcf71e 75:93c833019857
1 # EquidistantGrid is a grid with equidistant grid spacing per coordinat 1 # EquidistantGrid is a grid with equidistant grid spacing per coordinat
2 # direction. The domain is defined through the two points P1 = x̄₁, P2 = x̄₂ 2 # direction. The domain is defined through the two points P1 = x̄₁, P2 = x̄₂
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 struct EquidistantGrid <: AbstractGrid 6
7 numberOfPointsPerDim::Tuple # First coordinate direction stored first, then 7 struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
8 # second, then third. 8 numberOfPointsPerDim::NTuple{Dim, Int} # First coordinate direction stored first, then
9 limits::NTuple{2,Tuple} # Stores the two points which defines the range of 9
10 # the e.g (-1,0) and (1,2) for a domain of size 10 limit_lower::NTuple{Dim, T}
11 # (-1,1)x(0,2) 11 limit_upper::NTuple{Dim, T}
12 12
13 # General constructor 13 # General constructor
14 function EquidistantGrid(nPointsPerDim::Tuple, lims::NTuple{2,Tuple}) 14 function EquidistantGrid(nPointsPerDim::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
15 @assert length(nPointsPerDim) > 0 15 @assert all(nPointsPerDim.>0)
16 @assert count(x -> x > 0, nPointsPerDim) == length(nPointsPerDim) 16 @assert all(limit_upper.-limit_lower .!= 0)
17 @assert length(lims[1]) == length(nPointsPerDim) 17 return new{Dim,T}(nPointsPerDim, limit_lower, limit_upper)
18 @assert length(lims[2]) == length(nPointsPerDim)
19 # TODO: Assert that the same values are not passed in both lims[1] and lims[2]
20 # i.e the domain length is positive for all dimensions
21 return new(nPointsPerDim, lims)
22 end 18 end
23 # 1D constructor which can be called as EquidistantGrid(m, (xl,xr)) 19
24 function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real}) 20 # # 1D constructor which can be called as EquidistantGrid(m, (xl,xr))
25 return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],))) 21 # function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real})
26 end 22 # return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
23 # end
27 24
28 end 25 end
29 26
30 # Returns the number of dimensions of an EquidistantGrid. 27 # Returns the number of dimensions of an EquidistantGrid.
31 # 28 #
38 # Computes the total number of points of an EquidistantGrid. 35 # Computes the total number of points of an EquidistantGrid.
39 # 36 #
40 # @Input: grid - an EquidistantGrid 37 # @Input: grid - an EquidistantGrid
41 # @Return: numberOfPoints - The total number of points 38 # @Return: numberOfPoints - The total number of points
42 function numberOfPoints(grid::EquidistantGrid) 39 function numberOfPoints(grid::EquidistantGrid)
43 numberOfPoints = grid.numberOfPointsPerDim[1]; 40 return prod(grid.numberOfPointsPerDim)
44 for i = 2:length(grid.numberOfPointsPerDim);
45 numberOfPoints = numberOfPoints*grid.numberOfPointsPerDim[i]
46 end
47 return numberOfPoints
48 end 41 end
49 42
50 # Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance 43 # Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance
51 # between two points for each coordinate direction. 44 # between two points for each coordinate direction.
52 # 45 #
53 # @Input: grid - an EquidistantGrid 46 # @Input: grid - an EquidistantGrid
54 # @Return: h̄ - Grid spacing for each coordinate direction stored in a tuple. 47 # @Return: h̄ - Grid spacing for each coordinate direction stored in a tuple.
55 function spacings(grid::EquidistantGrid) 48 function spacings(grid::EquidistantGrid)
56 a = grid.limits[1] 49 return abs.(grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
57 b = grid.limits[2]
58
59 return abs.(b.-a)./(grid.numberOfPointsPerDim.-1)
60 end 50 end
61 51
62 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered 52 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered
63 # such that points in the first coordinate direction varies first, then the second 53 # such that points in the first coordinate direction varies first, then the second
64 # and lastely the third (if applicable) 54 # and lastely the third (if applicable)
65 # 55 #
66 # @Input: grid - an EquidistantGrid 56 # @Input: grid - an EquidistantGrid
67 # @Return: points - the points of the grid. 57 # @Return: points - the points of the grid.
68 function points(grid::EquidistantGrid) 58 function points(grid::EquidistantGrid)
69 # Compute signed grid spacings 59 dx̄ = (grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
70 dx̄ = Vector{Real}(undef, numberOfDimensions(grid))
71 for i ∈ eachindex(dx̄)
72 dx̄[i] = (grid.limits[2][i]-grid.limits[1][i])/(grid.numberOfPointsPerDim[i]-1)
73 end
74 dx̄ = Tuple(dx̄)
75 60
76 points = Vector{NTuple{numberOfDimensions(grid),Real}}(undef, numberOfPoints(grid)) 61 points = Vector{typeof(dx̄)}(undef, numberOfPoints(grid))
77 # Compute the points based on their Cartesian indices and the signed 62 # Compute the points based on their Cartesian indices and the signed
78 # grid spacings 63 # grid spacings
79 cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim) 64 cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
80 for i ∈ 1:numberOfPoints(grid) 65 for i ∈ 1:numberOfPoints(grid)
81 ci = Tuple(cartesianIndices[i]) .-1 66 ci = Tuple(cartesianIndices[i]) .-1
82 points[i] = grid.limits[1] .+ dx̄.*ci 67 points[i] = grid.limit_lower .+ dx̄.*ci
83 end 68 end
69
84 # TBD: Keep? this? How do we want to represent points in 1D? 70 # TBD: Keep? this? How do we want to represent points in 1D?
85 if numberOfDimensions(grid) == 1 71 if numberOfDimensions(grid) == 1
86 points = broadcast(x -> x[1], points) 72 points = broadcast(x -> x[1], points)
87 end 73 end
88 return points 74 return points