Mercurial > repos > public > sbplib_julia

comparison grid.jl @ 38:2dce28c59429

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix 1d version of plotgridfunction
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 10 Jan 2019 14:55:57 +0100
parents 32a53cbee6c5
children 5ec57ec148ef
comparison
equal deleted inserted replaced
28:32a53cbee6c5 38:2dce28c59429
1 module grid 1 module grid
2 using Plots 2 using PyPlot
3 3
4 abstract type Grid end 4 abstract type Grid end
5 5
6 function numberOfDimensions(grid::Grid) 6 function numberOfDimensions(grid::Grid)
7 error("Not implemented for abstact type Grid") 7 error("Not implemented for abstact type Grid")
39 # TODO: Assert that the same values are not passed in both lims[1] and lims[2] 39 # TODO: Assert that the same values are not passed in both lims[1] and lims[2]
40 # i.e the domain length is positive for all dimensions 40 # i.e the domain length is positive for all dimensions
41 return new(nPointsPerDim, lims) 41 return new(nPointsPerDim, lims)
42 end 42 end
43 # 1D constructor which can be called as EquidistantGrid(m, (xl,xr)) 43 # 1D constructor which can be called as EquidistantGrid(m, (xl,xr))
44 function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Integer}) 44 function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real})
45 return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],))) 45 return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
46 end 46 end
47 47
48 end 48 end
49 49
92 for i ∈ eachindex(dx̄) 92 for i ∈ eachindex(dx̄)
93 dx̄[i] = (grid.limits[2][i]-grid.limits[1][i])/(grid.numberOfPointsPerDim[i]-1) 93 dx̄[i] = (grid.limits[2][i]-grid.limits[1][i])/(grid.numberOfPointsPerDim[i]-1)
94 end 94 end
95 dx̄ = Tuple(dx̄) 95 dx̄ = Tuple(dx̄)
96 96
97 nPoints = numberOfPoints(grid) 97 points = Vector{NTuple{numberOfDimensions(grid),Real}}(undef, numberOfPoints(grid))
98 points = Vector{NTuple{numberOfDimensions(grid),Real}}(undef, nPoints)
99 # Compute the points based on their Cartesian indices and the signed 98 # Compute the points based on their Cartesian indices and the signed
100 # grid spacings 99 # grid spacings
101 cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim) 100 cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
102 for i ∈ 1:nPoints 101 for i ∈ 1:numberOfPoints(grid)
103 ci = Tuple(cartesianIndices[i]) .-1 102 ci = Tuple(cartesianIndices[i]) .-1
104 points[i] = grid.limits[1] .+ dx̄.*ci 103 points[i] = grid.limits[1] .+ dx̄.*ci
104 end
105 # TBD: Keep? this? How do we want to represent points in 1D?
106 if numberOfDimensions(grid) == 1
107 points = broadcast(x -> x[1], points)
105 end 108 end
106 return points 109 return points
107 end 110 end
108 111
109 function plotOnGrid(grid::EquidistantGrid,v::Vector) 112 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
110 dim = numberOfDimensions(grid) 113 @assert dim<=numberOfDimensions(grid)
111 x = points(grid) 114 @assert dim>0
115 points = range(grid.limits[1][dim],stop=grid.limits[2][dim],length=grid.numberOfPointsPerDim[dim])
116 end
112 117
113 if dim ==1 118 function plotgridfunction(grid::EquidistantGrid, gridfunction)
114 plot(x,v) 119 if numberOfDimensions(grid) == 1
120 plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
121 elseif numberOfDimensions(grid) == 2
122 x = pointsalongdim(grid,1)
123 X = repeat(x,1,grid.numberOfPointsPerDim[2])
124 y = pointsalongdim(grid,2)
125 Y = repeat(y,1,grid.numberOfPointsPerDim[1])'
126 elseif numberOfDimensions(grid) == 3
127 x = pointsalongdim(grid,1)
128 X = repeat(x,1,grid.numberOfPointsPerDim[2])
129 y = pointsalongdim(grid,2)
130 Y = repeat(y,1,grid.numberOfPointsPerDim[1])'
115 else 131 else
116 error(string("Plot not implemented for dim =", string(dim))) 132 error(string("Plot not implemented for dimension =", string(dim)))
117 end 133 end
118 end 134 end
119 135
120 end 136 end