Mercurial > repos > public > sbplib_julia

comparison grid.jl @ 21:2dbdd00eaea0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement function returing the points of a EquidistantGrid Additional: - Implement function calculating the grid spacing of a equidistant grid - Change members of EquidistantGrid from vectors to tuples - Add a 1D constructor for convenience
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 18 Dec 2018 13:27:25 +0100
parents af8469bc1cb3
children f2dc3e09fffc
comparison
equal deleted inserted replaced
20:d1ac68092138 21:2dbdd00eaea0
14 error("Not yet implemented") 14 error("Not yet implemented")
15 end 15 end
16 16
17 abstract type BoundaryId end 17 abstract type BoundaryId end
18 18
19 # Move to seperate file. 19 # TODO: Move to seperate file.
20 # Prefer to use UInt here, but printing UInt returns hex.
20 struct EquidistantGrid <: Grid 21 struct EquidistantGrid <: Grid
21 nPointsPerDim::Vector{Int} 22 numberOfPointsPerDim::Tuple
22 limits::Vector{Pair{Real, Real}} 23 limits::NTuple{2,Tuple} # Stores the points at the lower and upper corner of the domain.
23 function EquidistantGrid(nPointsPerDim, lims) 24 # e.g (-1,0) and (1,2) for a domain of size (-1,1)x(0,2)
24 @assert length(lims) == length(nPointsPerDim) 25
26 # General constructor
27 function EquidistantGrid(nPointsPerDim::Tuple, lims::NTuple{2,Tuple})
28 @assert length(nPointsPerDim) > 0
29 @assert count(x -> x > 0, nPointsPerDim) == length(nPointsPerDim)
30 @assert length(lims[1]) == length(nPointsPerDim)
31 @assert length(lims[2]) == length(nPointsPerDim)
32 # TODO: Assert that the same values are not passed in both lims[1] and lims[2]
33 # i.e the domain length is positive for all dimensions
25 return new(nPointsPerDim, lims) 34 return new(nPointsPerDim, lims)
26 end 35 end
36 # 1D constructor which can be called as EquidistantGrid(m, (x_l,x_r))
37 function EquidistantGrid(nPointsPerDim::Int, lims::NTuple{2,Int})
38 return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
39 end
40
27 end 41 end
28 42
29 function numberOfDimensions(grid::EquidistantGrid) 43 function numberOfDimensions(grid::EquidistantGrid)
30 return length(grid.nPointsPerDim) 44 return length(grid.numberOfPointsPerDim)
31 end 45 end
32 46
33 function numberOfPoints(grid::EquidistantGrid) 47 function numberOfPoints(grid::EquidistantGrid)
34 numberOfPoints = grid.nPointsPerDim[1]; 48 numberOfPoints = grid.numberOfPointsPerDim[1];
35 for i = 2:length(grid.nPointsPerDim); 49 for i = 2:length(grid.numberOfPointsPerDim);
36 numberOfPoints = numberOfPoints*grid.nPointsPerDim[i] 50 numberOfPoints = numberOfPoints*grid.numberOfPointsPerDim[i]
37 end 51 end
38 return numberOfPoints 52 return numberOfPoints
39 end 53 end
40 54
55 # TODO: Decide if spacings should be positive or if it is allowed to be negative
56 # If defined as positive, then need to do something extra when calculating the
57 # points. The current implementation works for
58 function spacings(grid::EquidistantGrid)
59 h = Vector{Real}(undef, numberOfDimensions(grid))
60 for i ∈ eachindex(h)
61 h[i] = (grid.limits[2][i]-grid.limits[1][i])/(grid.numberOfPointsPerDim[i]-1)
62 end
63 return Tuple(h)
64 end
65
41 function points(grid::EquidistantGrid) 66 function points(grid::EquidistantGrid)
42 points = Vector{Real}(undef, numberOfPoints(grid)) 67 nPoints = numberOfPoints(grid)
43 for i = 1:numberOfDimensions(grid) 68 points = Vector{NTuple{numberOfDimensions(grid),Real}}(undef, nPoints)
44 lims = limitsForDimension(grid,i) 69 cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
45 points = range(lims.first, stop=lims.second, length=grid.nPointsPerDim[i]) 70 for i ∈ 1:nPoints
71 ci = Tuple(cartesianIndices[i]) .-1
72 points[i] = grid.limits[1] .+ spacings(grid).*ci
46 end 73 end
47 return points 74 return points
48 end 75 end
49 76
50 function limitsForDimension(grid::EquidistantGrid, dim::Int)
51 return grid.limits[dim]
52 end 77 end
53
54 end