Line data Source code
1 : # https://stackoverflow.com/questions/51466537/how-to-plot-a-vector-field-in-julia
2 :
3 : """
4 : X, Y = meshgrid(t)
5 : X, Y = meshgrid(x, y)
6 :
7 : Replacement for Matlab's `meshgrid`.
8 : """
9 6 : meshgrid(x, y = x) =
10 : (repeat(x', inner=(length(y), 1)), repeat(y, inner=(1, length(x))))
11 :
12 : struct MeshGrid{T}
13 : X::Matrix{T}
14 : Y::Matrix{T}
15 :
16 7 : function MeshGrid{T}(X::Matrix{T}, Y::Matrix{T}) where {T}
17 7 : @assert size(X) == size(Y)
18 7 : new(X, Y)
19 : end
20 : end
21 :
22 7 : MeshGrid(X::Matrix{T}, Y::Matrix{T}) where {T} = MeshGrid{T}(X, Y)
23 :
24 6 : MeshGrid(g::Tuple{X, Y}) where {X, Y} = MeshGrid(collect(g[1]), collect(g[2]))
25 :
26 6 : MeshGrid(x, y = x) = MeshGrid(meshgrid(x, y))
27 :
28 : """
29 : grid = MeshGrid(X, Y)
30 : grid = MeshGrid(x, y = x)
31 : grid = MeshGrid(meshgrid(x, y))
32 :
33 : Generate data structure that captures contents generated by
34 : [`meshgrid`](@ref).
35 :
36 : The constructor calls [`meshgrid`](@ref) for non-matrix arguments
37 : `x, y`. A `MeshGrid` object can be deconstructed by
38 :
39 : X, Y = grid[]
40 :
41 : You can "enumerate" all vertex positions as `SVector{2, T}` with
42 :
43 : collect(grid[:])
44 : map(grid[:]) do x
45 : @SVector[x..., 0]
46 : end
47 :
48 : You can pass `grid::MeshGrid` to `Plots.plot`
49 :
50 : plot(grid, scatter=false, ...)
51 :
52 : See also [`meshgrid`](@ref), [`perturb!`](@ref), [`triangles`](@ref)
53 : """ MeshGrid
54 :
55 : """
56 : X, Y = MeshGrid(X, Y)
57 :
58 : Deconstruct [`MeshGrid`](@ref) to tuple `(X, Y)`.
59 : """
60 16 : Base.getindex(g::MeshGrid) = (g.X, g.Y)
61 :
62 : """
63 : collect(grid[:])
64 :
65 : Enumerate vertices of [`MeshGrid`](@ref).
66 : """
67 1 : function Base.getindex(g::MeshGrid, ::Colon)
68 1 : (SA[xy...] for xy in zip(vec(g.X), vec(g.Y)))
69 : end
70 :
71 : """
72 : m, n = size(grid)
73 :
74 : Get dimensions of [`MesgGrid`](@ref), same as
75 :
76 : m, n = size(first(grid[]))
77 : """
78 4 : function Base.size(g::MeshGrid)
79 4 : size(first(g[]))
80 : end
81 :
82 : """
83 : hx, hy = spacing(grid)
84 :
85 : Get distances between knots in x- and y-directions, i.e.,
86 : the cell size of the [`dual`](@ref) grid.
87 :
88 : See also [`MeshGrid`](@ref)
89 : """
90 5 : function spacing(g::MeshGrid)
91 5 : X, Y = g[]
92 5 : hx = X[1, 2] - X[1, 1]
93 5 : hy = Y[2, 1] - Y[1, 1]
94 5 : hx, hy
95 : end
96 :
97 : """
98 : d = dual(grid)
99 :
100 : Get dual [`MeshGrid`](@ref) such that each knot in `grid`
101 : refers to a cell in the dual `d`.
102 :
103 : See also [`primal`](@ref)
104 : """
105 1 : function dual(grid::MeshGrid)
106 1 : hx, hy = spacing(grid)
107 1 : X, Y = grid[]
108 1 : m, n = size(X)
109 :
110 1 : x = range(X[1, 1] - hx/2, X[1, end] + hx/2, length=n+1)
111 1 : y = range(Y[1, 1] - hy/2, Y[end, 1] + hy/2, length=m+1)
112 :
113 1 : MeshGrid(x, y)
114 : end
115 :
116 : """
117 : p = primal(grid)
118 :
119 : Get primal [`MeshGrid`](@ref) from [`dual`](@ref) `grid` such that
120 : every cell in `grid` refers to a knot in `p`.
121 : """
122 1 : function primal(grid::MeshGrid)
123 1 : hx, hy = spacing(grid)
124 1 : X, Y = grid[]
125 1 : m, n = size(X)
126 :
127 1 : x = range(X[1, 1] + hx/2, X[1, end] - hx/2, length=n-1)
128 1 : y = range(Y[1, 1] + hy/2, Y[end, 1] - hy/2, length=m-1)
129 :
130 1 : MeshGrid(x, y)
131 : end
132 :
133 : """
134 : grid = MeshGrid(x, y)
135 : x, y = ranges(grid)
136 :
137 : Get ranges `x` and `y` that generate [`MeshGrid`](@ref) `grid`.
138 : """
139 1 : function ranges(grid::MeshGrid)
140 1 : X, Y = grid[]
141 1 : m, n = size(X)
142 :
143 1 : x = range(X[1, 1], X[1, end], length=n)
144 1 : y = range(Y[1, 1], Y[end, 1], length=m)
145 :
146 1 : x, y
147 : end
148 :
149 : """
150 : perturb!(grid, p=0.05; boundary::Bool = true)
151 :
152 : Perturb vertex positions of grid with uniform random noise.
153 : `p` is the maximum amplitude relative to the cell size. If
154 : `boundary = false`, boundary vertices will not be perturbed.
155 :
156 : !!! warning
157 : This method does not check whether cells "flip orientation"!
158 :
159 : See also [`MeshGrid`](@ref)
160 : """
161 3 : function perturb!(g::MeshGrid, p = 0.05; boundary = true)
162 1 : R = CartesianIndices(g.X)
163 1 : Ifirst, Ilast = first(R), last(R)
164 1 : I1 = oneunit(Ifirst)
165 :
166 1 : if !boundary
167 1 : @assert all(size(g.X) .> 2)
168 1 : Ifirst += I1
169 1 : Ilast -= I1
170 : end
171 :
172 1 : hx = abs(g.X[1, 2] - g.X[1, 1]) * 2p
173 1 : hy = abs(g.Y[2, 1] - g.Y[1, 1]) * 2p
174 :
175 1 : for i in Ifirst:Ilast
176 1 : g.X[i] += (rand() - 0.5) * hx
177 1 : g.Y[i] += (rand() - 0.5) * hy
178 1 : end
179 : end
180 :
181 : """
182 : triangles(grid) :: Vector
183 :
184 : Collect triangles in `MeshGrid` with each cell split into
185 : two triangles. Each triangle is represented as a triple of
186 : linear vertex indices.
187 :
188 : See also [`MeshGrid`](@ref)
189 : """
190 1 : triangles(g::MeshGrid) = _triangles(size(g.X)...)
191 :
192 1 : function _triangles(m, n)
193 1 : sz = max(0, (m-1) * (n-1) * 2)
194 2 : t = Vector{Tuple{Int, Int, Int}}(undef, sz)
195 :
196 1 : L = LinearIndices((m, n))
197 :
198 1 : k = 1
199 1 : for i in 1:m-1, j in 1:n-1
200 2 : t[k ] = (L[i+1, j], L[i+1, j+1], L[i, j])
201 2 : t[k+1] = (L[i, j], L[i+1, j+1], L[i, j+1])
202 2 : k += 2
203 2 : end
204 :
205 1 : t
206 : end
|
