|
| 1 | +module MultivariateTriangle |
| 2 | + |
| 3 | +using ApproxFun, Makie, MultivariateOrthogonalPolynomials |
| 4 | +import AbstractPlotting: mesh, mesh! |
| 5 | +import ApproxFunOrthogonalPolynomials: _padua_length |
| 6 | + |
| 7 | +export contourf, contourf |
| 8 | + |
| 9 | +contourf(f::Fun; kwds...) = _mesh(meshdata(f)...; shading=false, kwds...) |
| 10 | +contourf!(s, f::Fun; kwds...) = _mesh!(s, meshdata(f)...; shading=false, kwds...) |
| 11 | + |
| 12 | + |
| 13 | +function _mesh(p, T, v; resolution=(400,400), kwds...) |
| 14 | + T_mat = Array{Int}(undef, length(T), 3) |
| 15 | + for k = 1:length(T) |
| 16 | + T_mat[k,:] .= T[k] |
| 17 | + end |
| 18 | + s = Scene(resolution=resolution) |
| 19 | + mesh!(s, [first.(p) last.(p)], T_mat; color=v, kwds...) |
| 20 | +end |
| 21 | + |
| 22 | + |
| 23 | +function _surface(p, T, v; resolution=(400,400), kwds...) |
| 24 | + T_mat = Array{Int}(undef, length(T), 3) |
| 25 | + for k = 1:length(T) |
| 26 | + T_mat[k,:] .= T[k] |
| 27 | + end |
| 28 | + # s = Scene(resolution=resolution) |
| 29 | + mesh(first.(p), last.(p), vec(v), T_mat; kwds...) |
| 30 | +end |
| 31 | + |
| 32 | + |
| 33 | + |
| 34 | +function _mesh!(s, p, T, v; kwds...) |
| 35 | + T_mat = Array{Int}(undef, length(T), 3) |
| 36 | + for k = 1:length(T) |
| 37 | + T_mat[k,:] .= T[k] |
| 38 | + end |
| 39 | + mesh!(s, [first.(p) last.(p)], T_mat; color=v, kwds...) |
| 40 | +end |
| 41 | + |
| 42 | +function meshdata(f::Fun{<:PiecewiseSpace}) |
| 43 | + pTv = MultivariateTriangle.meshdata.(components(f)) |
| 44 | + p = vcat(first.(pTv)...) |
| 45 | + T = pTv[1][2] |
| 46 | + cs = length(pTv[1][1]) |
| 47 | + for k = 2:length(pTv) |
| 48 | + append!(T, (t -> (cs.+t)).(pTv[k][2])) |
| 49 | + cs += length(pTv[k][1]) |
| 50 | + end |
| 51 | + |
| 52 | + v = vcat(last.(pTv)...) |
| 53 | + |
| 54 | + p, T, v |
| 55 | +end |
| 56 | + |
| 57 | +function meshdata(f::Fun{<:TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}}}) |
| 58 | + p = points(f) |
| 59 | + v = values(f) |
| 60 | + n = length(p) |
| 61 | + T = Vector{NTuple{3,Int}}() |
| 62 | + d_x,d_y = factors(domain(f)) |
| 63 | + a_x,b_x = endpoints(d_x) |
| 64 | + a_y,b_y = endpoints(d_y) |
| 65 | + if iseven(_padua_length(n)) |
| 66 | + l = floor(Int, (1+sqrt(1+8n))/4) |
| 67 | + |
| 68 | + push!(p, Vec(b_x,b_y)) |
| 69 | + push!(p, Vec(a_x,b_y)) |
| 70 | + |
| 71 | + push!(v, f(b_x,b_y)) |
| 72 | + push!(v, f(a_x,b_y)) |
| 73 | + |
| 74 | + for p = 0:l-2 |
| 75 | + for k = (2p*l)+1:(2p*l)+l-1 |
| 76 | + push!(T, (k+1, k, l+k+1)) |
| 77 | + end |
| 78 | + for k = (2p*l)+1:(2p*l)+l-1 |
| 79 | + push!(T, (k, l+k, l+k+1)) |
| 80 | + end |
| 81 | + for k = (2p*l)+l+1:(2p*l)+2l-1 |
| 82 | + push!(T, (k+1, k, l+k)) |
| 83 | + end |
| 84 | + for k = (2p*l)+l+2:(2p*l)+2l |
| 85 | + push!(T, (k, k+l-1, l+k)) |
| 86 | + end |
| 87 | + end |
| 88 | + for p=0:l-3 |
| 89 | + push!(T, ((2p+1)*l+1, (2p+2)*l+1, (2p+3)*l+1)) |
| 90 | + end |
| 91 | + for p =0:l-2 |
| 92 | + push!(T, ((2p+1)*l, (2p+2)*l, (2p+3)*l)) |
| 93 | + end |
| 94 | + push!(T, (1, n+1, l+1)) |
| 95 | + push!(T, (n-2l+1, n+2, n-l+1)) |
| 96 | + else |
| 97 | + l = floor(Int, (3+sqrt(1+8n))/4) |
| 98 | + |
| 99 | + push!(p, Vec(a_x,b_y)) |
| 100 | + push!(p, Vec(a_x,a_y)) |
| 101 | + |
| 102 | + push!(v, f(a_x,b_y)) |
| 103 | + push!(v, f(a_x,a_y)) |
| 104 | + |
| 105 | + for p = 0:l-2 |
| 106 | + for k = p*(2l-1)+1:p*(2l-1)+l-1 |
| 107 | + push!(T, (k+1, k, l+k)) |
| 108 | + end |
| 109 | + for k = p*(2l-1)+1:p*(2l-1)+l-2 |
| 110 | + push!(T, (k+1, l+k, l+k+1)) |
| 111 | + end |
| 112 | + end |
| 113 | + for p = 0:l-3 |
| 114 | + for k = p*(2l-1)+l+1:p*(2l-1)+2l-2 |
| 115 | + push!(T, (k+1, k, l+k)) |
| 116 | + end |
| 117 | + for k = p*(2l-1)+l+1:p*(2l-1)+2l-1 |
| 118 | + push!(T, (k, k+l-1, l+k)) |
| 119 | + end |
| 120 | + end |
| 121 | + |
| 122 | + for p=0:l-3 |
| 123 | + push!(T, (p*(2l-1) + 1, p*(2l-1) + l+1, p*(2l-1) + 2l)) |
| 124 | + end |
| 125 | + |
| 126 | + for p=0:l-3 |
| 127 | + push!(T, (p*(2l-1) + l, p*(2l-1) + 2l-1, p*(2l-1) + 3l-1)) |
| 128 | + end |
| 129 | + |
| 130 | + push!(T, (n-2l+2, n+1, n-l+2)) |
| 131 | + push!(T, (n-l+1, n+2, n)) |
| 132 | + end |
| 133 | + |
| 134 | + p, T, v |
| 135 | +end |
| 136 | +meshdata(f::Fun{<:Space{<:Triangle}}) = |
| 137 | + triangle_meshdata(points(f), values(f), (domain(f).a, domain(f).b, domain(f).c), |
| 138 | + f.((domain(f).a, domain(f).b, domain(f).c))) |
| 139 | + |
| 140 | +function triangle_meshdata(p, v, (a, b, c), (fa, fb, fc)) |
| 141 | + n = length(p) |
| 142 | + T = Vector{NTuple{3,Int}}() |
| 143 | + |
| 144 | + |
| 145 | + if iseven(_padua_length(n)) |
| 146 | + l = floor(Int, (1+sqrt(1+8n))/4) |
| 147 | + |
| 148 | + push!(p, b) |
| 149 | + push!(p, c) |
| 150 | + |
| 151 | + push!(v, fb) |
| 152 | + push!(v, fc) |
| 153 | + |
| 154 | + for p = 0:l-2 |
| 155 | + for k = (2p*l)+1:(2p*l)+l-1 |
| 156 | + push!(T, (k+1, k, l+k+1)) |
| 157 | + end |
| 158 | + for k = (2p*l)+1:(2p*l)+l-1 |
| 159 | + push!(T, (k, l+k, l+k+1)) |
| 160 | + end |
| 161 | + for k = (2p*l)+l+1:(2p*l)+2l-1 |
| 162 | + push!(T, (k+1, k, l+k)) |
| 163 | + end |
| 164 | + for k = (2p*l)+l+2:(2p*l)+2l |
| 165 | + push!(T, (k, k+l-1, l+k)) |
| 166 | + end |
| 167 | + end |
| 168 | + for p=0:l-3 |
| 169 | + push!(T, ((2p+1)*l+1, (2p+2)*l+1, (2p+3)*l+1)) |
| 170 | + end |
| 171 | + for p =0:l-2 |
| 172 | + push!(T, ((2p+1)*l, (2p+2)*l, (2p+3)*l)) |
| 173 | + end |
| 174 | + push!(T, (1, n+1, l+1)) |
| 175 | + push!(T, (n-2l+1, n+2, n-l+1)) |
| 176 | + else |
| 177 | + l = floor(Int, (3+sqrt(1+8n))/4) |
| 178 | + |
| 179 | + push!(p, Vec(c)) |
| 180 | + push!(p, Vec(a)) |
| 181 | + |
| 182 | + push!(v, fc) |
| 183 | + push!(v, fa) |
| 184 | + |
| 185 | + for p = 0:l-2 |
| 186 | + for k = p*(2l-1)+1:p*(2l-1)+l-1 |
| 187 | + push!(T, (k+1, k, l+k)) |
| 188 | + end |
| 189 | + for k = p*(2l-1)+1:p*(2l-1)+l-2 |
| 190 | + push!(T, (k+1, l+k, l+k+1)) |
| 191 | + end |
| 192 | + end |
| 193 | + for p = 0:l-3 |
| 194 | + for k = p*(2l-1)+l+1:p*(2l-1)+2l-2 |
| 195 | + push!(T, (k+1, k, l+k)) |
| 196 | + end |
| 197 | + for k = p*(2l-1)+l+1:p*(2l-1)+2l-1 |
| 198 | + push!(T, (k, k+l-1, l+k)) |
| 199 | + end |
| 200 | + end |
| 201 | + |
| 202 | + for p=0:l-3 |
| 203 | + push!(T, (p*(2l-1) + 1, p*(2l-1) + l+1, p*(2l-1) + 2l)) |
| 204 | + end |
| 205 | + |
| 206 | + for p=0:l-3 |
| 207 | + push!(T, (p*(2l-1) + l, p*(2l-1) + 2l-1, p*(2l-1) + 3l-1)) |
| 208 | + end |
| 209 | + |
| 210 | + push!(T, (n-2l+2, n+1, n-l+2)) |
| 211 | + push!(T, (n-l+1, n+2, n)) |
| 212 | + end |
| 213 | + |
| 214 | + p, T, v |
| 215 | +end |
| 216 | + |
| 217 | + |
| 218 | + |
| 219 | +end # module |
0 commit comments