Mercurial > repos > public > sbplib_julia

comparison EquidistantGrid.jl @ 51:614b56a017b9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 11 Jan 2019 11:55:13 +0100
parents
children 8c6db1f6d8e0
comparison
equal deleted inserted replaced
46:50c6c252d954 51:614b56a017b9
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̄₂
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)
5 # the domain is defined as (-1,1)x(0,2).
6 struct EquidistantGrid <: AbstractGrid
7 numberOfPointsPerDim::Tuple # First coordinate direction stored first, then
8 # second, then third.
9 limits::NTuple{2,Tuple} # Stores the two points which defines the range of
10 # the e.g (-1,0) and (1,2) for a domain of size
11 # (-1,1)x(0,2)
12
13 # General constructor
14 function EquidistantGrid(nPointsPerDim::Tuple, lims::NTuple{2,Tuple})
15 @assert length(nPointsPerDim) > 0
16 @assert count(x -> x > 0, nPointsPerDim) == length(nPointsPerDim)
17 @assert length(lims[1]) == length(nPointsPerDim)
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
23 # 1D constructor which can be called as EquidistantGrid(m, (xl,xr))
24 function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real})
25 return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
26 end
27
28 end
29
30 # Returns the number of dimensions of an EquidistantGrid.
31 #
32 # @Input: grid - an EquidistantGrid
33 # @Return: numberOfPoints - The number of dimensions
34 function numberOfDimensions(grid::EquidistantGrid)
35 return length(grid.numberOfPointsPerDim)
36 end
37
38 # Computes the total number of points of an EquidistantGrid.
39 #
40 # @Input: grid - an EquidistantGrid
41 # @Return: numberOfPoints - The total number of points
42 function numberOfPoints(grid::EquidistantGrid)
43 numberOfPoints = grid.numberOfPointsPerDim[1];
44 for i = 2:length(grid.numberOfPointsPerDim);
45 numberOfPoints = numberOfPoints*grid.numberOfPointsPerDim[i]
46 end
47 return numberOfPoints
48 end
49
50 # Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance
51 # between two points for each coordinate direction.
52 #
53 # @Input: grid - an EquidistantGrid
54 # @Return: h̄ - Grid spacing for each coordinate direction stored in a tuple.
55 function spacings(grid::EquidistantGrid)
56 h̄ = Vector{Real}(undef, numberOfDimensions(grid))
57 for i ∈ eachindex(h̄)
58 h̄[i] = abs(grid.limits[2][i]-grid.limits[1][i])/(grid.numberOfPointsPerDim[i]-1)
59 end
60 return Tuple(h̄)
61 end
62
63 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered
64 # such that points in the first coordinate direction varies first, then the second
65 # and lastely the third (if applicable)
66 #
67 # @Input: grid - an EquidistantGrid
68 # @Return: points - the points of the grid.
69 function points(grid::EquidistantGrid)
70 # Compute signed grid spacings
71 dx̄ = Vector{Real}(undef, numberOfDimensions(grid))
72 for i ∈ eachindex(dx̄)
73 dx̄[i] = (grid.limits[2][i]-grid.limits[1][i])/(grid.numberOfPointsPerDim[i]-1)
74 end
75 dx̄ = Tuple(dx̄)
76
77 points = Vector{NTuple{numberOfDimensions(grid),Real}}(undef, numberOfPoints(grid))
78 # Compute the points based on their Cartesian indices and the signed
79 # grid spacings
80 cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
81 for i ∈ 1:numberOfPoints(grid)
82 ci = Tuple(cartesianIndices[i]) .-1
83 points[i] = grid.limits[1] .+ dx̄.*ci
84 end
85 # TBD: Keep? this? How do we want to represent points in 1D?
86 if numberOfDimensions(grid) == 1
87 points = broadcast(x -> x[1], points)
88 end
89 return points
90 end
91
92 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
93 @assert dim<=numberOfDimensions(grid)
94 @assert dim>0
95 points = range(grid.limits[1][dim],stop=grid.limits[2][dim],length=grid.numberOfPointsPerDim[dim])
96 end
97
98 using PyPlot, PyCall
99 # using Plots; pyplot()
100
101 function plotgridfunction(grid::EquidistantGrid, gridfunction)
102 if numberOfDimensions(grid) == 1
103 plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
104 elseif numberOfDimensions(grid) == 2
105 mx = grid.numberOfPointsPerDim[1];
106 my = grid.numberOfPointsPerDim[2];
107 x = pointsalongdim(grid,1)
108 X = repeat(x,1,my)
109 y = pointsalongdim(grid,2)
110 Y = repeat(y,1,mx)
111 # plot_surface(X,Y,reshape(gridfunction,mx,my))
112 fig = figure()
113 ax = fig[:add_subplot](1,1,1, projection = "3d")
114 ax[:plot_surface](X,Y,reshape(gridfunction,mx,my))
115 plt[:show]()
116 else
117 error(string("Plot not implemented for dimension ", string(numberOfDimensions(grid))))
118 end
119 end