diff --git a/Project-1/example/demo.jl b/Project-1/example/demo.jl index b7ca259..60e283b 100644 --- a/Project-1/example/demo.jl +++ b/Project-1/example/demo.jl @@ -41,24 +41,21 @@ julia> yt0012(0.13) """ Base.@pure function yt0012(ξ::T where {T<:Real}) ε = 0.0050544107511712 - return 5 * 0.12 * ( - if 0 <= ξ <= 1 - 0.2969 * sqrt(ξ) - - 0.1260 * ξ - - 0.3516 * ξ^2 + - 0.2843 * ξ^3 - - 0.1015 * ξ^4 - elseif 1 < ξ <= 1+ε # 2-order smooth joining of the tail edge. - 0.03423616658680458 * sqrt(1+ε - ξ) + - 0.12491972188290168 * (1+ε - ξ) + - 11.579398976610927 * (1+ε - ξ)^2 + - 12.21202728162769 * (1+ε - ξ)^3 - else - 0 - end) + return 5 * + 0.12 * + (if 0 <= ξ <= 1 + 0.2969 * sqrt(ξ) - 0.1260 * ξ - 0.3516 * ξ^2 + 0.2843 * ξ^3 - 0.1015 * ξ^4 + elseif 1 < ξ <= 1 + ε # 2-order smooth joining of the tail edge. + 0.03423616658680458 * sqrt(1 + ε - ξ) + + 0.12491972188290168 * (1 + ε - ξ) + + 11.579398976610927 * (1 + ε - ξ)^2 + + 12.21202728162769 * (1 + ε - ξ)^3 + else + 0 + end) end -const Mx, My = 100+1, 20+1 +const Mx, My = 100 + 1, 20 + 1 include(normpath(joinpath(@__DIR__, "../../src/misc-util.jl"))) using .__CFD2021__misc_util__: tuplejoin @@ -127,8 +124,12 @@ interpolated_x = Array{Float64}(undef, (length(v) for v in x_lo)...); interpolated_y = Array{Float64}(undef, (length(v) for v in y_lo)...); @btime transfinite_interpolate_2d!(interpolated_x, x_lo, x_hi) @btime transfinite_interpolate_2d!(interpolated_y, y_lo, y_hi) -for u in axes(interpolated_x)[begin] plot!(p,interpolated_x[u,:], interpolated_y[u,:], label = nothing, color = myblue); end -for u in axes(interpolated_x)[ end ] plot!(p,interpolated_x[:,u], interpolated_y[:,u], label = nothing, color = myblue); end +for u in axes(interpolated_x)[begin] + plot!(p, interpolated_x[u, :], interpolated_y[u, :]; label=nothing, color=myblue) +end +for u in axes(interpolated_x)[end] + plot!(p, interpolated_x[:, u], interpolated_y[:, u]; label=nothing, color=myblue) +end # for u in axes(interpolated_x)[begin] plot!(p,interpolated_x[u,:],-interpolated_y[u,:], label = nothing, color = myblue); end # for u in axes(interpolated_x)[ end ] plot!(p,interpolated_x[:,u],-interpolated_y[:,u], label = nothing, color = myblue); end plot!(p, yt0012, 0, 1; label=nothing, color=:red) diff --git a/Project-1/src/transfinite-interpolate.jl b/Project-1/src/transfinite-interpolate.jl index 7d32ffa..79a03d9 100644 --- a/Project-1/src/transfinite-interpolate.jl +++ b/Project-1/src/transfinite-interpolate.jl @@ -70,10 +70,10 @@ function transfinite_interpolate_2d(f_lo::Tuple{Function,Function}, function interpolating_function(u::NTuple{2,T} where {T<:Number})::Number #= linear-part from edges bilinear-part from vertices =# - (1-u[2])*f_lo[1](u[1]) - (1-u[1])*(1-u[2])*v_cr[1,1] + - ( u[2])*f_hi[1](u[1]) - ( u[1])*(1-u[2])*v_cr[1,2] + - (1-u[1])*f_lo[2](u[2]) - (1-u[1])*( u[2])*v_cr[2,1] + - ( u[1])*f_hi[2](u[2]) - ( u[1])*( u[2])*v_cr[2,2] + return (1 - u[2]) * f_lo[1](u[1]) - (1 - u[1]) * (1 - u[2]) * v_cr[1, 1] + + (u[2]) * f_hi[1](u[1]) - (u[1]) * (1 - u[2]) * v_cr[1, 2] + + (1 - u[1]) * f_lo[2](u[2]) - (1 - u[1]) * (u[2]) * v_cr[2, 1] + + (u[1]) * f_hi[2](u[2]) - (u[1]) * (u[2]) * v_cr[2, 2] end return interpolating_function end @@ -107,60 +107,76 @@ julia> transfinite_interpolate_2d(([0,0.5,1], [0,0.5,1]), ([1,1.5,2], [1,1.5,2]) 1.0 1.5 2.0 ``` """ -function transfinite_interpolate_2d( # for regular Arrays - v_lo::Tuple{Vector{T}, Vector{S}} where {T <: Real, S <: Real}, - v_hi::Tuple{Vector{T}, Vector{S}} where {T <: Real, S <: Real}, - v_cr::Union{Nothing, SMatrix{2, 2, T} where T <: Real} = nothing)::Array{Real, 2} +function transfinite_interpolate_2d(v_lo::Tuple{Vector{T}, + Vector{S}} where {T<:Real,S<:Real}, + v_hi::Tuple{Vector{T}, + Vector{S}} where {T<:Real,S<:Real}; + v_cr::Union{Nothing,SMatrix{2,2,T} where T<:Real}=nothing)::Array{Real, + 2} #= check dimension-compatibility before allocating anything =# - length(v_lo) == length(v_hi) || throw( - DimensionMismatch(" incompatible grid specification ($(length(v_lo))-d at one end but $(length(v_hi))-d at the other).")) - N = length(v_lo); errmsg = "" #= In this function, dimensions are always indexed from 1, but not so for each axis of them =# - for dim in eachindex(v_lo) if ! (axes(v_lo[dim]) == axes(v_hi[dim])) - errmsg *= "grid mismatch in dimension $(dim) ($(axes(v_lo[dim])) at one end but $(axes(v_hi[dim])) at the other)." - end end; errmsg == "" || throw(DimensionMismatch(errmsg)) + length(v_lo) == length(v_hi) || + throw(DimensionMismatch(" incompatible grid specification ($(length(v_lo))-d at one end but $(length(v_hi))-d at the other).")) + N = length(v_lo) + errmsg = "" #= In this function, dimensions are always indexed from 1, but not so for each axis of them =# + for dim in eachindex(v_lo) + if !(axes(v_lo[dim]) == axes(v_hi[dim])) + errmsg *= "grid mismatch in dimension $(dim) ($(axes(v_lo[dim])) at one end but $(axes(v_hi[dim])) at the other)." + end + end + errmsg == "" || throw(DimensionMismatch(errmsg)) interpolated_array = Array{Float64}(undef, (length(v) for v in v_lo)...) transfinite_interpolate_2d!(interpolated_array, v_lo, v_hi, v_cr) return interpolated_array end -function transfinite_interpolate_2d( # for OffsetArrays - v_lo::Tuple{OffsetArray, OffsetArray}, - v_hi::Tuple{OffsetArray, OffsetArray}, - v_cr::Union{Nothing, SMatrix{2, 2, T} where T <: Real} = nothing)::OffsetArray +function transfinite_interpolate_2d(v_lo::Tuple{OffsetArray,OffsetArray}, + v_hi::Tuple{OffsetArray,OffsetArray}; + v_cr::Union{Nothing,SMatrix{2,2,T} where T<:Real}=nothing)::OffsetArray #= check dimension-compatibility before allocating anything =# - length(v_lo) == length(v_hi) || throw( - DimensionMismatch(" incompatible grid specification ($(length(v_lo))-d at one end but $(length(v_hi))-d at the other).")) - N = length(v_lo); errmsg = "" #= In this function, dimensions are always indexed from 1, but not so for each axis of them =# - for dim in eachindex(v_lo) if ! (axes(v_lo[dim]) == axes(v_hi[dim])) - errmsg *= "grid mismatch in dimension $(dim) ($(axes(v_lo[dim])) at one end but $(axes(v_hi[dim])) at the other)." - end end; errmsg == "" || throw(DimensionMismatch(errmsg)) + length(v_lo) == length(v_hi) || + throw(DimensionMismatch(" incompatible grid specification ($(length(v_lo))-d at one end but $(length(v_hi))-d at the other).")) + N = length(v_lo) + errmsg = "" #= In this function, dimensions are always indexed from 1, but not so for each axis of them =# + for dim in eachindex(v_lo) + if !(axes(v_lo[dim]) == axes(v_hi[dim])) + errmsg *= "grid mismatch in dimension $(dim) ($(axes(v_lo[dim])) at one end but $(axes(v_hi[dim])) at the other)." + end + end + errmsg == "" || throw(DimensionMismatch(errmsg)) interpolated_array = Array{Float64}(undef, (length(v) for v in v_lo)...) - interpolated_array = OffsetArray(interpolated_array, tuplejoin((axes(v) for v in v_lo)...)) + interpolated_array = OffsetArray(interpolated_array, + tuplejoin((axes(v) for v in v_lo)...)) transfinite_interpolate_2d!(interpolated_array, v_lo, v_hi, v_cr) return interpolated_array end -function transfinite_interpolate_2d!( # for regular Arrays - interpolated_array::Array{T, 2} where T <: Real, - v_lo::Tuple{Vector{T}, Vector{S}} where {T <: Real, S <: Real}, - v_hi::Tuple{Vector{T}, Vector{S}} where {T <: Real, S <: Real}, - v_cr::Union{Nothing, SMatrix{2, 2, T} where T <: Real} = nothing) - +function transfinite_interpolate_2d!(interpolated_array::Array{T,2} where {T<:Real}, + v_lo::Tuple{Vector{T}, + Vector{S}} where {T<:Real,S<:Real}, + v_hi::Tuple{Vector{T}, + Vector{S}} where {T<:Real,S<:Real}, + v_cr::Union{Nothing,SMatrix{2,2,T} where T<:Real}=nothing) if typeof(v_cr) <: Nothing # e.g, the following lines still assumes range ``1:2``. - v_cr = SA[(v_lo[1][begin]+v_lo[2][begin])/2 (v_lo[1][ end ]+v_hi[2][begin])/2 - (v_hi[1][begin]+v_lo[2][ end ])/2 (v_hi[1][ end ]+v_hi[2][ end ])/2] - end; M = [length(v)-1 for v in v_lo] # this is only for scaling, and never for iterating + v_cr = SA[(v_lo[1][begin] + v_lo[2][begin])/2 (v_lo[1][end] + v_hi[2][begin])/2 + (v_hi[1][begin] + v_lo[2][end])/2 (v_hi[1][end] + v_hi[2][end])/2] + end + M = [length(v) - 1 for v in v_lo] # this is only for scaling, and never for iterating # loop Version for u in CartesianIndices(interpolated_array) #= #= Deprecated =# Iterators.product(IndexCartesian(), [eachindex(v) for v in v_lo]...) =# interpolated_array[u] = #= linear-part from edges bilinear-part from vertices =# - (1-(u[2]-1)/M[2])*v_lo[1][u[1]] - (1-(u[1]-1)/M[1])*(1-(u[2]-1)/M[2])*v_cr[1,1] + - ( (u[2]-1)/M[2])*v_hi[1][u[1]] - ( (u[1]-1)/M[1])*(1-(u[2]-1)/M[2])*v_cr[1,2] + - (1-(u[1]-1)/M[1])*v_lo[2][u[2]] - (1-(u[1]-1)/M[1])*( (u[2]-1)/M[2])*v_cr[2,1] + - ( (u[1]-1)/M[1])*v_hi[2][u[2]] - ( (u[1]-1)/M[1])*( (u[2]-1)/M[2])*v_cr[2,2] + (1 - (u[2] - 1) / M[2]) * v_lo[1][u[1]] - + (1 - (u[1] - 1) / M[1]) * + (1 - (u[2] - 1) / M[2]) * + v_cr[1, 1] + ((u[2] - 1) / M[2]) * v_hi[1][u[1]] - + ((u[1] - 1) / M[1]) * (1 - (u[2] - 1) / M[2]) * v_cr[1, 2] + + (1 - (u[1] - 1) / M[1]) * v_lo[2][u[2]] - + (1 - (u[1] - 1) / M[1]) * ((u[2] - 1) / M[2]) * v_cr[2, 1] + + ((u[1] - 1) / M[1]) * v_hi[2][u[2]] - + ((u[1] - 1) / M[1]) * ((u[2] - 1) / M[2]) * v_cr[2, 2] end return interpolated_array @@ -170,18 +186,23 @@ function transfinite_interpolate_2d!(interpolated_array::OffsetArray, v_hi::Tuple{OffsetArray,OffsetArray}, v_cr::Union{Nothing,SMatrix{2,2,T} where T<:Real}=nothing) if typeof(v_cr) <: Nothing # e.g, the following lines still assumes range ``1:2``. - v_cr = SA[(v_lo[1][begin]+v_lo[2][begin])/2 (v_lo[1][ end ]+v_hi[2][begin])/2 - (v_hi[1][begin]+v_lo[2][ end ])/2 (v_hi[1][ end ]+v_hi[2][ end ])/2] - end; M = [length(v)-1 for v in v_lo] # this is only for scaling, and never for iterating + v_cr = SA[(v_lo[1][begin] + v_lo[2][begin])/2 (v_lo[1][end] + v_hi[2][begin])/2 + (v_hi[1][begin] + v_lo[2][end])/2 (v_hi[1][end] + v_hi[2][end])/2] + end + M = [length(v) - 1 for v in v_lo] # this is only for scaling, and never for iterating # loop Version for u in CartesianIndices(interpolated_array) #= #= Deprecated =# Iterators.product(IndexCartesian(), [eachindex(v) for v in v_lo]...) =# interpolated_array[u] = #= linear-part from edges bilinear-part from vertices =# - (1-u[2]/M[2])*v_lo[1][u[1]] - (1-u[1]/M[1])*(1-u[2]/M[2])*v_cr[1,1] + - ( u[2]/M[2])*v_hi[1][u[1]] - ( u[1]/M[1])*(1-u[2]/M[2])*v_cr[1,2] + - (1-u[1]/M[1])*v_lo[2][u[2]] - (1-u[1]/M[1])*( u[2]/M[2])*v_cr[2,1] + - ( u[1]/M[1])*v_hi[2][u[2]] - ( u[1]/M[1])*( u[2]/M[2])*v_cr[2,2] + (1 - u[2] / M[2]) * v_lo[1][u[1]] - + (1 - u[1] / M[1]) * (1 - u[2] / M[2]) * v_cr[1, 1] + + (u[2] / M[2]) * v_hi[1][u[1]] - + (u[1] / M[1]) * (1 - u[2] / M[2]) * v_cr[1, 2] + + (1 - u[1] / M[1]) * v_lo[2][u[2]] - + (1 - u[1] / M[1]) * (u[2] / M[2]) * v_cr[2, 1] + + (u[1] / M[1]) * v_hi[2][u[2]] - + (u[1] / M[1]) * (u[2] / M[2]) * v_cr[2, 2] end return interpolated_array end diff --git a/Project-1/test/runtests.jl b/Project-1/test/runtests.jl index 780a42d..958d09b 100644 --- a/Project-1/test/runtests.jl +++ b/Project-1/test/runtests.jl @@ -16,7 +16,6 @@ using Test using BenchmarkTools using OffsetArrays - @testset "transfinite_interpolate" begin include(normpath(joinpath(@__DIR__, "../../src/misc-util.jl"))) using .__CFD2021__misc_util__: tuplejoin @@ -26,57 +25,75 @@ using OffsetArrays @testset "transfinite_interpolate_2d" begin Mx, My = 3, 3 - @testset "on functions" begin - interp_func = transfinite_interpolate_2d( (x->2x, y->3y), (x->4x+3, y->2+5y) ) - @test interp_func((0.5,0.5)) == 3.0 + interp_func = transfinite_interpolate_2d((x -> 2x, y -> 3y), + (x -> 4x + 3, y -> 2 + 5y)) + @test interp_func((0.5, 0.5)) == 3.0 end @testset "on arrays" begin - x_lo = ([LinRange(0, 1, Mx)...], [LinRange(0, 0, My)...]); x_hi = ([LinRange(0, 1, Mx)...], [LinRange(1, 1, My)...]); - y_lo = ([LinRange(0, 0, Mx)...], [LinRange(0, 1, My)...]); y_hi = ([LinRange(1, 1, Mx)...], [LinRange(0, 1, My)...]); - in_situ_interp_x = Array{Float64}(undef, (length(v) for v in x_lo)...); - in_situ_interp_y = Array{Float64}(undef, (length(v) for v in y_lo)...); + x_lo = ([LinRange(0, 1, Mx)...], [LinRange(0, 0, My)...]) + x_hi = ([LinRange(0, 1, Mx)...], [LinRange(1, 1, My)...]) + y_lo = ([LinRange(0, 0, Mx)...], [LinRange(0, 1, My)...]) + y_hi = ([LinRange(1, 1, Mx)...], [LinRange(0, 1, My)...]) + in_situ_interp_x = Array{Float64}(undef, (length(v) for v in x_lo)...) + in_situ_interp_y = Array{Float64}(undef, (length(v) for v in y_lo)...) transfinite_interpolate_2d!(in_situ_interp_x, x_lo, x_hi) transfinite_interpolate_2d!(in_situ_interp_y, y_lo, y_hi) ab_situ_interp_x = transfinite_interpolate_2d(x_lo, x_hi) ab_situ_interp_y = transfinite_interpolate_2d(y_lo, y_hi) @testset "on ordinary arrays" begin - @test in_situ_interp_x == ab_situ_interp_x == [0.0 0.0 0.0; 0.5 0.5 0.5; 1.0 1.0 1.0] - @test in_situ_interp_y == ab_situ_interp_y == [0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0] + @test in_situ_interp_x == + ab_situ_interp_x == + [0.0 0.0 0.0; 0.5 0.5 0.5; 1.0 1.0 1.0] + @test in_situ_interp_y == + ab_situ_interp_y == + [0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0] end - x_lo = map(x->OffsetArray(x, OffsetArrays.Origin(0)), x_lo) - y_lo = map(x->OffsetArray(x, OffsetArrays.Origin(0)), y_lo) - x_hi = map(x->OffsetArray(x, OffsetArrays.Origin(0)), x_hi) - y_hi = map(x->OffsetArray(x, OffsetArrays.Origin(0)), y_hi) - in_situ_interp_x = OffsetArray(in_situ_interp_x, tuplejoin((axes(v) for v in x_lo)...)) - in_situ_interp_y = OffsetArray(in_situ_interp_y, tuplejoin((axes(v) for v in y_lo)...)) + x_lo = map(x -> OffsetArray(x, OffsetArrays.Origin(0)), x_lo) + y_lo = map(x -> OffsetArray(x, OffsetArrays.Origin(0)), y_lo) + x_hi = map(x -> OffsetArray(x, OffsetArrays.Origin(0)), x_hi) + y_hi = map(x -> OffsetArray(x, OffsetArrays.Origin(0)), y_hi) + in_situ_interp_x = OffsetArray(in_situ_interp_x, + tuplejoin((axes(v) for v in x_lo)...)) + in_situ_interp_y = OffsetArray(in_situ_interp_y, + tuplejoin((axes(v) for v in y_lo)...)) transfinite_interpolate_2d!(in_situ_interp_x, x_lo, x_hi) transfinite_interpolate_2d!(in_situ_interp_y, y_lo, y_hi) ab_situ_interp_x = transfinite_interpolate_2d(x_lo, x_hi) ab_situ_interp_y = transfinite_interpolate_2d(y_lo, y_hi) @testset "on OffsetArrays" begin - @test in_situ_interp_x == ab_situ_interp_x != [0.0 0.0 0.0; 0.5 0.5 0.5; 1.0 1.0 1.0] - @test in_situ_interp_y == ab_situ_interp_y != [0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0] - @test in_situ_interp_x == ab_situ_interp_x == OffsetArray([0.0 0.0 0.0; 0.5 0.5 0.5; 1.0 1.0 1.0], OffsetArrays.Origin(0)) - @test in_situ_interp_y == ab_situ_interp_y == OffsetArray([0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0], OffsetArrays.Origin(0)) + @test in_situ_interp_x == + ab_situ_interp_x != + [0.0 0.0 0.0; 0.5 0.5 0.5; 1.0 1.0 1.0] + @test in_situ_interp_y == + ab_situ_interp_y != + [0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0] + @test in_situ_interp_x == + ab_situ_interp_x == + OffsetArray([0.0 0.0 0.0; 0.5 0.5 0.5; 1.0 1.0 1.0], + OffsetArrays.Origin(0)) + @test in_situ_interp_y == + ab_situ_interp_y == + OffsetArray([0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0], + OffsetArrays.Origin(0)) end @testset "exception handling" begin - x_lo = ([LinRange(0, 1, Mx)...], [LinRange(0, 0, My)...]); x_hi = ([LinRange(0, 1, 2Mx+1)...], [LinRange(1, 1, 2My+1)...]); + x_lo = ([LinRange(0, 1, Mx)...], [LinRange(0, 0, My)...]) + x_hi = ([LinRange(0, 1, 2Mx + 1)...], [LinRange(1, 1, 2My + 1)...]) @test_throws DimensionMismatch transfinite_interpolate_2d(x_lo, x_hi) end - end # @testset "on arrays" end # @testset "transfinite_interpolate_2d" end # @testset "transfinite_interpolate" @testset "elliptic_grid_gen" begin include("../src/elliptic-grid-gen.jl") - import .elliptic_grid_gen + using .elliptic_grid_gen: elliptic_grid_gen @testset "inverted_poisson_jacobi" begin @test_skip elliptic_grid_gen.inverted_poisson_jacobi_step diff --git a/deps/build.jl b/deps/build.jl index acd4cf9..dc3695b 100644 --- a/deps/build.jl +++ b/deps/build.jl @@ -12,4 +12,4 @@ # See the License for the specific language governing permissions and # limitations under the License. using Dates -println(now(), "\tEmpty build of the project") \ No newline at end of file +println(now(), "\tEmpty build of the project") diff --git a/test/runtests.jl b/test/runtests.jl index d3e104a..4b37b03 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,9 +21,10 @@ include(normpath(joinpath(@__DIR__, "../src/CFD2021Projects.jl"))) @testset "tuplejoin" begin using .__CFD2021__misc_util__: tuplejoin @test tuplejoin((1, 2)) == (1, 2) - temp_v = ([0,0.5,1], [0,0.5,1], [0,0.5,1], [0,0.5,1]) + temp_v = ([0, 0.5, 1], [0, 0.5, 1], [0, 0.5, 1], [0, 0.5, 1]) temp_t = (axes(v) for v in temp_v) - @test tuplejoin(temp_t...) == (Base.OneTo(3), Base.OneTo(3), Base.OneTo(3), Base.OneTo(3)) + @test tuplejoin(temp_t...) == + (Base.OneTo(3), Base.OneTo(3), Base.OneTo(3), Base.OneTo(3)) end # @testset "tuplejoin" end # @testset "misc_utils"