Mercurial > repos > public > sbplib_julia

comparison grid.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 50c6c252d954
children cb9a789338a1
comparison
equal deleted inserted replaced
46:50c6c252d954 51:614b56a017b9
1 module grid 1 module Grid
2 using Plots
3 pyplot()
4 2
5 abstract type Grid end 3 # TODO: Where is this used?
6
7 function numberOfDimensions(grid::Grid)
8 error("Not implemented for abstact type Grid")
9 end
10
11 function numberOfPoints(grid::Grid)
12 error("Not implemented for abstact type Grid")
13 end
14
15 function points(grid::Grid)
16 error("Not implemented for abstact type Grid")
17 end
18
19 # TODO: Should this be here?
20 abstract type BoundaryId end 4 abstract type BoundaryId end
21 5
22 # EquidistantGrid is a grid with equidisant grid spacing per coordinat 6 include("AbstractGrid.jl")
23 # direction. The domain is defined through the two points P1 = x̄₁, P2 = x̄₂ 7 include("EquidistantGrid.jl")
24 # by the exterior product of the vectors obtained by projecting (x̄₂-x̄₁) onto
25 # the coordinate directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2)
26 # the domain is defined as (-1,1)x(0,2).
27 struct EquidistantGrid <: Grid
28 numberOfPointsPerDim::Tuple # First coordinate direction stored first, then
29 # second, then third.
30 limits::NTuple{2,Tuple} # Stores the two points which defines the range of
31 # the e.g (-1,0) and (1,2) for a domain of size
32 # (-1,1)x(0,2)
33
34 # General constructor
35 function EquidistantGrid(nPointsPerDim::Tuple, lims::NTuple{2,Tuple})
36 @assert length(nPointsPerDim) > 0
37 @assert count(x -> x > 0, nPointsPerDim) == length(nPointsPerDim)
38 @assert length(lims[1]) == length(nPointsPerDim)
39 @assert length(lims[2]) == length(nPointsPerDim)
40 # TODO: Assert that the same values are not passed in both lims[1] and lims[2]
41 # i.e the domain length is positive for all dimensions
42 return new(nPointsPerDim, lims)
43 end
44 # 1D constructor which can be called as EquidistantGrid(m, (xl,xr))
45 function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real})
46 return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
47 end
48 8
49 end 9 end
50
51 # Returns the number of dimensions of an EquidistantGrid.
52 #
53 # @Input: grid - an EquidistantGrid
54 # @Return: numberOfPoints - The number of dimensions
55 function numberOfDimensions(grid::EquidistantGrid)
56 return length(grid.numberOfPointsPerDim)
57 end
58
59 # Computes the total number of points of an EquidistantGrid.
60 #
61 # @Input: grid - an EquidistantGrid
62 # @Return: numberOfPoints - The total number of points
63 function numberOfPoints(grid::EquidistantGrid)
64 numberOfPoints = grid.numberOfPointsPerDim[1];
65 for i = 2:length(grid.numberOfPointsPerDim);
66 numberOfPoints = numberOfPoints*grid.numberOfPointsPerDim[i]
67 end
68 return numberOfPoints
69 end
70
71 # Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance
72 # between two points for each coordinate direction.
73 #
74 # @Input: grid - an EquidistantGrid
75 # @Return: h̄ - Grid spacing for each coordinate direction stored in a tuple.
76 function spacings(grid::EquidistantGrid)
77 h̄ = Vector{Real}(undef, numberOfDimensions(grid))
78 for i ∈ eachindex(h̄)
79 h̄[i] = abs(grid.limits[2][i]-grid.limits[1][i])/(grid.numberOfPointsPerDim[i]-1)
80 end
81 return Tuple(h̄)
82 end
83
84 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered
85 # such that points in the first coordinate direction varies first, then the second
86 # and lastely the third (if applicable)
87 #
88 # @Input: grid - an EquidistantGrid
89 # @Return: points - the points of the grid.
90 function points(grid::EquidistantGrid)
91 # Compute signed grid spacings
92 dx̄ = Vector{Real}(undef, numberOfDimensions(grid))
93 for i ∈ eachindex(dx̄)
94 dx̄[i] = (grid.limits[2][i]-grid.limits[1][i])/(grid.numberOfPointsPerDim[i]-1)
95 end
96 dx̄ = Tuple(dx̄)
97
98 points = Vector{NTuple{numberOfDimensions(grid),Real}}(undef, numberOfPoints(grid))
99 # Compute the points based on their Cartesian indices and the signed
100 # grid spacings
101 cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
102 for i ∈ 1:numberOfPoints(grid)
103 ci = Tuple(cartesianIndices[i]) .-1
104 points[i] = grid.limits[1] .+ dx̄.*ci
105 end
106 # TBD: Keep? this? How do we want to represent points in 1D?
107 if numberOfDimensions(grid) == 1
108 points = broadcast(x -> x[1], points)
109 end
110 return points
111 end
112
113 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
114 @assert dim<=numberOfDimensions(grid)
115 @assert dim>0
116 points = range(grid.limits[1][dim],stop=grid.limits[2][dim],length=grid.numberOfPointsPerDim[dim])
117 end
118
119 function plotgridfunction(grid::EquidistantGrid, gridfunction)
120 if numberOfDimensions(grid) == 1
121 plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
122 elseif numberOfDimensions(grid) == 2
123 mx = grid.numberOfPointsPerDim[1];
124 my = grid.numberOfPointsPerDim[2];
125 x = pointsalongdim(grid,1)
126 X = repeat(x,1,my)
127 y = pointsalongdim(grid,2)
128 Y = repeat(y,1,mx)'
129 surface(X,Y,reshape(gridfunction,mx,my))
130 else
131 error(string("Plot not implemented for dimension ", string(numberOfDimensions(grid))))
132 end
133 end
134
135 # Evaluate function f on the grid g
136 function evalOn(g::Grid, f::Function)
137 F(x) = f(x...)
138 return F.(points(g))
139 end
140
141 end