Mercurial > repos > public > sbplib_julia

annotate EquidistantGrid.jl @ 77:2be36b38389d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use limit_lower and limit_upper in pointsalongdim
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 24 Jan 2019 13:33:07 +0100
parents 93c833019857
children 8d505e9bc715
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
1 # EquidistantGrid is a grid with equidistant grid spacing per coordinat
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
2 # direction. The domain is defined through the two points P1 = x̄₁, P2 = x̄₂
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
3 # by the exterior product of the vectors obtained by projecting (x̄₂-x̄₁) onto
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
4 # the coordinate directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
5 # the domain is defined as (-1,1)x(0,2).
75
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
6
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
7 struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
8 numberOfPointsPerDim::NTuple{Dim, Int} # First coordinate direction stored first, then
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
9
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
10 limit_lower::NTuple{Dim, T}
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
11 limit_upper::NTuple{Dim, T}
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
12
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
13 # General constructor
75
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
14 function EquidistantGrid(nPointsPerDim::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
15 @assert all(nPointsPerDim.>0)
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
16 @assert all(limit_upper.-limit_lower .!= 0)
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
17 return new{Dim,T}(nPointsPerDim, limit_lower, limit_upper)
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
18 end
75
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
19
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
20 # # 1D constructor which can be called as EquidistantGrid(m, (xl,xr))
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
21 # function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real})
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
22 # return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
23 # end
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
24
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
25 end
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
26
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
27 # Returns the number of dimensions of an EquidistantGrid.
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
28 #
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
29 # @Input: grid - an EquidistantGrid
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
30 # @Return: numberOfPoints - The number of dimensions
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
31 function numberOfDimensions(grid::EquidistantGrid)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
32 return length(grid.numberOfPointsPerDim)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
33 end
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
34
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
35 # Computes the total number of points of an EquidistantGrid.
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
36 #
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
37 # @Input: grid - an EquidistantGrid
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
38 # @Return: numberOfPoints - The total number of points
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
39 function numberOfPoints(grid::EquidistantGrid)
75
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
40 return prod(grid.numberOfPointsPerDim)
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
41 end
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
42
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
43 # Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
44 # between two points for each coordinate direction.
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
45 #
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
46 # @Input: grid - an EquidistantGrid
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
47 # @Return: h̄ - Grid spacing for each coordinate direction stored in a tuple.
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
48 function spacings(grid::EquidistantGrid)
75
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
49 return abs.(grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
50 end
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
52 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
53 # such that points in the first coordinate direction varies first, then the second
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
54 # and lastely the third (if applicable)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
55 #
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
56 # @Input: grid - an EquidistantGrid
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
57 # @Return: points - the points of the grid.
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
58 function points(grid::EquidistantGrid)
75
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
59 dx̄ = (grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
60
75
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
61 points = Vector{typeof(dx̄)}(undef, numberOfPoints(grid))
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
62 # Compute the points based on their Cartesian indices and the signed
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
63 # grid spacings
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
64 cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
65 for i ∈ 1:numberOfPoints(grid)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
66 ci = Tuple(cartesianIndices[i]) .-1
75
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
67 points[i] = grid.limit_lower .+ dx̄.*ci
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
68 end
75
93c833019857 Make EquidistantGrid more concrete
Jonatan Werpers <jonatan@werpers.com>
parents: 74
diff changeset
69
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
70 # TBD: Keep? this? How do we want to represent points in 1D?
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
71 if numberOfDimensions(grid) == 1
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
72 points = broadcast(x -> x[1], points)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
73 end
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
74 return points
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
75 end
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
76
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
77 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
78 @assert dim<=numberOfDimensions(grid)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
79 @assert dim>0
77
2be36b38389d Use limit_lower and limit_upper in pointsalongdim
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 75
diff changeset
80 points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.numberOfPointsPerDim[dim])
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
81 end
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
82
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
83 using PyPlot, PyCall
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
84
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
85 function plotgridfunction(grid::EquidistantGrid, gridfunction)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
86 if numberOfDimensions(grid) == 1
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
87 plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
88 elseif numberOfDimensions(grid) == 2
58
f3ea15c8ebb6 Clean up plotgridfunction now that the matplotlib backend works. The issue was not in the julia code, but rather the installation of matplotlib.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 53
diff changeset
89 mx = grid.numberOfPointsPerDim[1]
f3ea15c8ebb6 Clean up plotgridfunction now that the matplotlib backend works. The issue was not in the julia code, but rather the installation of matplotlib.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 53
diff changeset
90 my = grid.numberOfPointsPerDim[2]
f3ea15c8ebb6 Clean up plotgridfunction now that the matplotlib backend works. The issue was not in the julia code, but rather the installation of matplotlib.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 53
diff changeset
91 X = repeat(pointsalongdim(grid,1),1,my)
f3ea15c8ebb6 Clean up plotgridfunction now that the matplotlib backend works. The issue was not in the julia code, but rather the installation of matplotlib.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 53
diff changeset
92 Y = permutedims(repeat(pointsalongdim(grid,2),1,mx))
f3ea15c8ebb6 Clean up plotgridfunction now that the matplotlib backend works. The issue was not in the julia code, but rather the installation of matplotlib.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 53
diff changeset
93 plot_surface(X,Y,reshape(gridfunction,mx,my));
51
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
94 else
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
95 error(string("Plot not implemented for dimension ", string(numberOfDimensions(grid))))
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
96 end
614b56a017b9 Split grid.jl into AbstractGrid.jl and EquidistantGrid.jl
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
97 end